nitrogen.linalg

Linear algebra routines, including iterative eigensolvers. This module uses the SciPy LinearOperator class to represent matrix-vector product routines.

class nitrogen.linalg.RealOrthogonalProjector(*args, **kwargs)

Bases: scipy.sparse.linalg.interface.LinearOperator

A projection operator to orthogonalize against a given set of vectors

\[\hat{P} = 1 - U U^T\]
Parameters

U ((N,n) ndarray) – n column vectors to orthogonalize against. The columns of U are assumed to be real and orthonormal.

nitrogen.linalg.aslinearoperator(A)

Return A as a LinearOperator.

‘A’ may be any of the following types:
  • ndarray

  • matrix

  • sparse matrix (e.g. csr_matrix, lil_matrix, etc.)

  • LinearOperator

  • An object with .shape and .matvec attributes

See the LinearOperator documentation for additional information.

Notes

If ‘A’ has no .dtype attribute, the data type is determined by calling LinearOperator.matvec() - set the .dtype attribute to prevent this call upon the linear operator creation.

Examples

>>> from scipy.sparse.linalg import aslinearoperator
>>> M = np.array([[1,2,3],[4,5,6]], dtype=np.int32)
>>> aslinearoperator(M)
<2x3 MatrixLinearOperator with dtype=int32>
nitrogen.linalg.bounds(H, k=10, m=3)

Estimate lower and upper bounds to the spectral range of the Hermitian operator H.

Parameters
  • H (LinearOperator) – A Hermitian operator.

  • k (int, optional) – The number of Lanczos steps to take before estimating the bounds. The default is 10.

  • m (int, optional) – The safe-guard level for the bound estimate. This must be less than or equal to k. The default is 3.

Returns

  • low (float) – An estimate of the lower bound of the smallest eigenvalue.

  • high (float) – An estimate of the upper bound of the largest eigenvalue.

Notes

These spectral bound estimates are based the work of Y. Zhou and R.-C. Li [ZL11]. Their bounds (2.6), (2.7), and (2.8) correspond to m = 1, k, and 3, respectively.

In general, a larger value of m will provide a safer, but looser, bound on the spectral range.

References

ZL11

Y. Zhou and R.-C. Li, “Bound the spectrum of large Hermitian matrices.” Linear Algebra and its Applications, 435, 480-493 (2011). https://doi.org/10.1016/j.laa.2010.06.034

nitrogen.linalg.chebauto(H, K, v0=None)

Calculate the auto-correlation function of a wavepacket via the Chebyshev propagator of of a Hermitian operator H.

Parameters
  • H (LinearOperator) – A Hermitian operator.

  • K (int) – The number of propagation steps.

  • v0 (ndarray, optional) – The initial vector. If None, a uniform vector is used. The default is None.

Returns

  • C (ndarray) – A (2*K+1,) array containing the auto-correlation function of v0 in the Chebyshev order domain.

  • scale ((float, float)) – The mean energy and half-width, (ebar, de) used to scale the Hamiltonian. The unscaled energy is related to the Chebyshev angle parameter as E = de * cos(theta) + ebar.

See also

bounds

Used to scale the Hamiltonian

chebspec

Calculate a spectrum from the Chebyshev auto-correlation function.

Notes

See Ref. [CG99] for a description of Chebyshev propagators.

References

CG99(1,2)

R. Chen, H. Guo. “The Chebyshev propagator for quantum systems.” Computer Physics Communications, 119, 19-31 (1999). https://doi.org/10.1016/S0010-4655(98)00179-9

nitrogen.linalg.chebspec(H, K, v0=None, window='gaussian', window_scale=1.0, sample_factor=1, norm='area')

Calculate the auto-correlation function of a wavepacket via the Chebyshev propagator of of a Hermitian operator H.

Parameters
  • H (LinearOperator) – A Hermitian operator.

  • K (int) – The number of propagation steps.

  • v0 (ndarray, optional) – The initial vector. If None, a uniform vector is used. The default is None.

  • window ({'gaussian', None}) – The window function. If None, no window function is used. The default is ‘gaussian’. See Notes for details.

  • window_scale (float) – A scaling parameter for the window function. The default is 1.0. See Notes for details.

  • sample_factor (int) – The over-sampling factor. Must be an integer >= 1. The default is 1.

  • norm ({'area', 'peak'}) – The normalization convention for the energy spectrum. The default is ‘area’. See Notes for details.

Returns

  • E (ndarray) – The energy axis.

  • G (ndarray) – The energy spectrum, normalized according to norm.

  • fwhm (ndarray) – The energy-dependent FWHM line width.

  • C (ndarray) – The windowed auto-correlation function.

See also

bounds

Used to scale the Hamiltonian

chebauto

Calculates the Chebyshev auto-correlation function.

chebwindow

Standard windows functions.

Notes

This spectral method follows that described in Ref. [CG99].

If the norm parameter is ‘area’ (default), then the returned spectrum G is normalized such that its energy integral is equal to the norm-squared of the initial vector, v0.

The non-linear relationship between energy and the Chebyshev angle lead to an energy-dependent line width. Therefore, the peak amplitude of a resolved transition will not be proportional to the spectral intensity globally. The norm is ‘peak’, then the spectrum is rescaled such that the peak amplitude of individual, resolved transitions does correspond to the spectral intensity. The line widths are not affected, however, such that the integrated area now no longer has a meaningful norm.

nitrogen.linalg.chebwindow(N, window, window_scale=1.0)

Calculate the window function for a Chebyshev order-domain correlation function.

Parameters
  • N (int) – Length of correlation function.

  • window ({'gaussian', 'boxcar', None}) – The window function type.

  • window_scale (float) – A scaling parameter for the window function. Must be >= 1.0 The default is 1.0. See Notes for details.

Returns

  • f (ndarray) – The window function.

  • hw (float) – The approximate half-width-half-maximum value (in units of radians).

Notes

The window functions are always scaled to a dimensionless order-time parameter, tau, ranging from 0 to 1 over the total propagation time. The ‘gaussian’ window is equal to np.exp(-(3.0 * window_scale * tau)**2). The ‘boxcar’ window is 1 for tau < 1/window_scale and 0 for tau >- 1/window_scale.

nitrogen.linalg.eigstrp(H, k=5, pad=10, tol=1e-10, maxiter=None, v0=None, rper=20, P=None, pper=1, printlevel=0, eigval='smallest', projection_level=2, refine_per=20)

Thick-restart Lanczos iterative eigensolver with projection.

Parameters
  • H (LinearOperator) – A real symmetric or complex Hermitian operator.

  • k (int, optional) – The number of eigenvalues requested. The default is 5.

  • pad (int, optional) – The total number of Lanczos vectors kept after a restart is k + pad. If pad is too small, the convergence will be less than optimal. The default is 10.

  • tol (float, optional) – The absolute eigenvalue convergence tolerance. The default is 1e-10.

  • maxiter (int, optional) – The maximum number of iterations. If None, no maximum is set. The default is None.

  • v0 (ndarray, optional) – Initial vector. The default is random (None).

  • rper (int, optional) – The restart period. The default is 20.

  • P (LinearOperator, optional) – Projection operator applied to the Lanczos vectors. The default is None.

  • pper (int, optional) – The projection period. The default is 1.

  • printlevel (int, optional) – Print level. The default is 0.

  • eigval ({'smallest', 'largest'}, optional) – Calculat the smallest or largest eigenvalues. The default is ‘smallest’.

  • projection_level (integer, optional) – The projection level determines how aggressively the projection operator is applied at each iteration. If 2 (default), it will be applied immediately after the matrix-vector routine and again after orthogonalization. If 1, it will be applied only after orthogonalization. If 0, it will not be applied.

  • refine_per (integer, optional) – After refine_per restarts, the target Ritz vectors will be explicitly rediagonalized and reorthogonalized. The default is 20.

Returns

  • w (ndarray) – Array of k eigenvalues

  • v (ndarray) – Array of k eigenvectors. v[:,i] is the eigenvector corresponding to eigenvalue w[i].

Notes

The restarted Lanczos method is based on that described by Wu and Simon [WS00].

References

WS00

K. Wu and H. Simon, “Thick-Restart Lanczos Method for Large Symmetric Eigenvalue Problems.” SIAM Journal on Matrix Analysis and Applications, 22(2), 602-616. (2000) https://doi.org/10.1137/S0895479898334605

nitrogen.linalg.full(H)

Return the full matrix of a LinearOperator

Parameters

H (LinearOperator) – A matrix operator in LinearOperator form.

Returns

Hfull – The full, dense matrix.

Return type

ndarray

nitrogen.linalg.msqrth(A)

Calculate the square-root matrix of a hermitian matrix, A.

Parameters

A (array_like) – The matrix

Returns

rtA – The square-root matrix.

Return type

ndarray

Notes

The square root is calculated via eigendecomposition. The square root of each eigenvalue is determined via NumPy’s sqrt function, which defines the branch convention.

nitrogen.linalg.simple_corr(H, dt, n, v0)

Simple finite time-step propagation for an autocorrelation function.

Parameters
  • H (ndarray) – The Hamiltonian matrix.

  • dt (float) – The time step (in units of \(t/\hbar\))

  • n (integer) – The number of time steps.

  • v0 (ndarray) – The initial vector.

Returns

C – The autocorrelation function.

Return type

ndarray

nitrogen.linalg.packed

Packed-storage matrices and routines. The lower triangle is stored in row-major order. (For symmetric matrices, this is equivalent to upper triangle column-major order. For Hermitian matrices, this is the conjugate of the upper triangle in column-major order.)

nitrogen.linalg.packed.IJ2k(I, J)

Assuming a symmetric array, calculate the 1D packed storage index for the (I,J) = (J,I) element

nitrogen.linalg.packed.chol_sp(A, out=None)

Cholesky decomposition of a symmetric matrix in packed format. If real symmetric, then A should be positive definite. If complex symmetric (not Hermitian), then H should have non-zero pivots.

Parameters
  • A ((np, ...) ndarray) – A is stored in 1D packed format (see nitrogen.linalg.packed)

  • out ((np, ...) ndarray) – Output buffer. If None, this will be created. If out = H, then in-place decomposition is performed

Returns

out – The lower triangle Cholesky decomposition L in packed storage. H = L @ L.T

Return type

ndarray

nitrogen.linalg.packed.inv_sp(A, out=None)

Inverse of a real symmetric, positive definite matrix in packed format. (Complex symmetric matrices may not be defined.)

Parameters
  • A ((np, ...) ndarray) – The packed storage array.

  • out ((np, ...) ndarray) – Output buffer. If None, this will be created. If out = A, then in-place inversion is performed

Returns

out – The result.

Return type

ndarray

nitrogen.linalg.packed.inv_tp(L, out=None)

Invert a triangular matrix in lower row-packed (or upper column-packed) storage.

Parameters
  • L ((np, ...) ndarray) – L is stored in 1D packed format (see nitrogen.linalg.packed)

  • out ((np, ...) ndarray) – Output buffer. If None, this will be created. If out = L, then in-place inversion is performed

Returns

out – The inverse of the triangular matrix in packed storage.

Return type

ndarray

nitrogen.linalg.packed.k2IJ(k)

Calculate the full 2D index (I,`J`) for a given packed 1D index k. The returned index is always in the lower triangle.

Parameters

k (int) – The packed 1D index.

Returns

I,J – The unpacked 2D index.

Return type

np.uint64

nitrogen.linalg.packed.llt_tp(L, out=None)

L @ L.T of a lower triangular matrix.

Parameters
  • L ((np, ...) ndarray) – Lower triangular matrix in packed storage.

  • out ((np, ...) ndarray) – Output buffer. If None, this will be created. If out = L, then in-place multiplication is performed

Returns

out – The symmetric result in packed storage.

Return type

ndarray

nitrogen.linalg.packed.ltl_tp(L, out=None)

L.T @ L with a lower triangular matrix.

Parameters
  • L ((np, ...) ndarray) – Lower triangular matrix in packed storage.

  • out ((np, ...) ndarray) – Output buffer. If None, this will be created. If out = L, then in-place multiplication is performed

Returns

out – The symmetric result in packed storage.

Return type

ndarray

nitrogen.linalg.packed.n2N(n)

Calculate the square matrix rank N for a packed storage size n

Parameters

n (int) – The packed length.

Returns

N – The matrix rank.

Return type

np.uint64

nitrogen.linalg.packed.symfull(P)

Return the full array for a symmetric array in packed storage

Parameters

P (ndarray) – The lower triangle of a symmetric array in packed storage.

Returns

The full symmetric matrix

Return type

ndarray

nitrogen.linalg.packed.trAB_sp(A, B, out=None)

The trace of A @ B, in symmetric packed stroage.

Parameters
  • A ((np, ...) ndarray) – Symmetric matrix in packed storage.

  • B ((np, ...) ndarray) – Symmetric matrix in packed storage.

  • out ((...) ndarray, optional) – Output buffer. If None, this will be created. If scalar, this is ignored.

Returns

out – The result.

Return type

ndarray or scalar

nitrogen.linalg.packed.trilfull(L)

Return the full array for a lower triangle array. in packed storage

Parameters

L (ndarray) – The lower triangle in packed storage

Returns

The full lower triangle array.

Return type

ndarray