dymad.numerics.linalg¶

Functions

`check_direction`(v1, v2)	The cosine values between v1 and v2.
`check_orthogonality`(U, V[, M])	Expecting U.H * M * V = I
`eig_low_rank`(U, V)	Approach like Exact DMD.
`expm_full_rank`(W, s, b)	Compute B_i = b @ exp(s_i * W) for i=1..m, in batch.
`expm_low_rank`(U, V, s, b)	Compute B_i = b @ exp(s_i * U V^T) for i=1..m, in batch.
`logm_low_rank`(V, U[, dt])	Given A = V U^T dt (n x n, rank r) with V, U (n x r), compute logm(A).
`make_random_matrix`(Ndim, Nrnk, zrng, wrng[, dt])	Random (Ndim x Ndim) matrix of rank Nrnk, with randomized eigenvalues ranged in zrng and wrng.
`mode_split`(l, U[, comp])	Split eigenvalues and modes according to the requested components.
`randomized_svd`(loader, N, k[, oversample, ...])	Two-pass randomized SVD for a matrix stored as row-blocks:
`real_lowrank_from_eigpairs`(lam, U, V[, tol])	Construct a real low-rank factorization V_real @ B @ U_real.T from (possibly complex) eigenpairs of a real matrix A = sum_i lam_i * v_i * u_i^H.
`scaled_eig`(A[, B])	Suppose A U = B U L, V^H A = L V^H B Ideally one should have double diagonalization (for non-degenerate case): V^H B U = I and V^H A U = L but by default each column of U and V is normalized by the length, and the double diagonalization is not satisfied.
`truncate_sequence`(seq, order)	Truncation of scalar sequence.
`truncated_lstsq`(A, B[, tsvd])	Solve the linear system AX = B by least squares.
`truncated_svd`(X, order)	A vanilla interface for different types of truncation order.

dymad.numerics.linalg.check_direction(v1, v2)¶: The cosine values between v1 and v2.

dymad.numerics.linalg.check_orthogonality(U, V, M=None)¶: Expecting U.H * M * V = I

dymad.numerics.linalg.eig_low_rank(U, V)¶

Approach like Exact DMD.

Suppose the full matrix is A = U V^T, with U, V of shape (n, r).

We compute the eigendecomposition of the small matrix A_tilde = V^T U of shape (r, r)

A_tilde = W @ L @ W_inv

Then A = (U W) @ L @ (V W_inv)^H

Furthermore, we scale the left and right eigenvectors so that they are orthonormal.

dymad.numerics.linalg.expm_full_rank(W, s, b)¶

Compute B_i = b @ exp(s_i * W) for i=1..m, in batch.

Parameters:

W (Tensor) – (n, n) full rank matrix
s (Tensor) – (m,) list/1D tensor of scalars
b (Tensor) – (batch, n) rows are left-multipliers

Returns:

(m, batch, n) where out[i] = b @ exp(s[i] * W)

Return type:

out

dymad.numerics.linalg.expm_low_rank(U, V, s, b)¶

Compute B_i = b @ exp(s_i * U V^T) for i=1..m, in batch.

Uses the identity: exp(sA) = I + U [s * phi_1(s S)] V^T, S = V^T U. So: b @ exp(sA) = b + (bU) [s * phi_1(s S)] V^T.

Parameters:

U (Tensor) – (n, r)
V (Tensor) – (n, r)
s (Tensor) – (m,) list/1D tensor of scalars
b (Tensor) – (batch, n) rows are left-multipliers

Returns:

(m, batch, n) where out[i] = b @ exp(s[i] * U V^T)

Return type:

out

dymad.numerics.linalg.logm_low_rank(V, U, dt=1.0)¶

Given A = V U^T dt (n x n, rank r) with V, U (n x r), compute logm(A).

Technically, this logarithm is ill-defined if A is not full rank, but here we compute the following:

Suppose A has the eigendecomposition A = W @ L @ R^H and we let logm(A) = W @ (log(L)/dt) @ R^H

For computational purpose, we return two real factors of logm(A), with the help of real_lowrank_from_eigpairs.

Notes

This approach does not work when A has negative real eigenvalues.

Parameters:

V (ndarray) – (n, r) real arrays Tall factors of A = V U^T. Columns need not be orthonormal.
U (ndarray) – (n, r) real arrays Tall factors of A = V U^T. Columns need not be orthonormal.

Returns:

(n, r) real arrays: Tall factors of logm(A).

Return type:

V_out, U_out

dymad.numerics.linalg.make_random_matrix(Ndim, Nrnk, zrng, wrng, dt=-1)¶: Random (Ndim x Ndim) matrix of rank Nrnk, with randomized eigenvalues ranged in zrng and wrng. If dt>0 is given, the eigenvalues will be mapped to discrete-time. The eigenpairs are always assumed to be conjugate.

dymad.numerics.linalg.mode_split(l, U, comp='ri')¶

Split eigenvalues and modes according to the requested components.

Consider r pairs of eigenvalues and modes, modes in shape (r, n), this function splits them as (r, k) and (r, k, n), where k is the number of requested components.

For example, if comp=’ri’, then k=2, and the two components are the real and imaginary parts of the eigenvalues and modes.

Parameters:

l (ndarray) – (r,) eigenvalues
U (ndarray) – (r, n) modes
comp (str) – ‘r’ - real, ‘i’ - imag, ‘a’ - amplitude, ‘p’ - phase, can be composed like ‘ri’, ‘ap’ etc.; default is ‘ri’.

Returns:

(r, k) U_split: (r, k, n)

Return type:

l_split

dymad.numerics.linalg.randomized_svd(loader, N, k, oversample=10, n_iter=0, return_u=False, dtype=<class 'numpy.float64'>, seed=0)¶

Two-pass randomized SVD for a matrix stored as row-blocks:

A^T = [A_0^T, A_1^T, …, A_{N-1}^T]

where A_i shape (n_i, d), A shape (N=sum_i n_i, d)

Two passes estimate k singular values and right singular vectors. A third pass is needed if left singular vectors are requested.

Parameters:

loader (callable) – A function that takes a block index and returns the corresponding block.
N (int) – Number of blocks
k (int) – Target rank (k > 0).
oversample (int, default 10) – Extra dimensions to improve spectral accuracy (total sketch dim l = k + oversample).
n_iter (int, default 0) – Number of power iterations (each adds two more passes).
return_u (bool, default False) – Whether to return left singular vectors.
dtype (numpy dtype, default np.float64) – Internal precision for computation.
seed (int or None, default 0) – RNG seed for the Gaussian sketch.

Returns:

(N, k) ndarray, Left singular vectors (columns). S : (k,) ndarray, Singular values in descending order. Vt : (d, k) ndarray, Right singular vectors (columns).

Return type:

dymad.numerics.linalg.real_lowrank_from_eigpairs(lam, U, V, tol=1e-12)¶

Construct a real low-rank factorization V_real @ B @ U_real.T from (possibly complex) eigenpairs of a real matrix A = sum_i lam_i * v_i * u_i^H.

Parameters:

lam (ndarray) – (r,) eigenvalues
U (ndarray) – (n, r) left eigenvectors
V (ndarray) – (n, r) right eigenvectors
tol (float) – tolerance to identify real vs complex eigenvalues

Returns:

(r_real, r_real) block-diagonal (1x1 for real, 2x2 for conjugate pairs) U_real: (m, r_real) V_real: (n, r_real)

Return type:

dymad.numerics.linalg.scaled_eig(A, B=None)¶

Suppose A U = B U L, V^H A = L V^H B Ideally one should have double diagonalization (for non-degenerate case): V^H B U = I and V^H A U = L but by default each column of U and V is normalized by the length, and the double diagonalization is not satisfied. Here we scale both U and V so that they are approximately orthonormal to each other (w.r.t. B); also the scaling is such that the norms of u_i and v_i are equal.

However, if one needs to project quantities to, e.g., U, use pseudo-inverse of U instead of V for numerical robustness.

dymad.numerics.linalg.truncate_sequence(seq, order)¶

Truncation of scalar sequence.

Possible order parameters

Float: Max value to retain
Integer: Keep first N values
‘full’: Retain all pairs

dymad.numerics.linalg.truncated_lstsq(A, B, tsvd=None)¶

Solve the linear system AX = B by least squares.

If truncated SVD is used, the function returns the two factors of X.

Parameters:

A (np.ndarray) – Coefficient matrix.
B (np.ndarray) – Right-hand side matrix.
tsvd (int or float, optional) – If provided, use truncated SVD with this order.

Returns:

Solution matrix X, or its two factors.

Return type:

np.ndarray

dymad.numerics.linalg.truncated_svd(X, order)¶

A vanilla interface for different types of truncation order.

Possible order parameters

Float, positive: Energy percentage
Float, negative: Optimal truncation by Gavish&Donoho.
Integer, positive: Keep first N pairs
Integer, negative: Remove last N pairs
‘full’: Retain all pairs