dymad.numerics¶

class dymad.numerics.DM(n_components, n_neighbors, alpha=1, epsilon=None)¶

Bases: object

Diffusion map

Original MATLAB implementation by Tyrus Berry & John Harlim

This implementation

is class-based

accelerates epsilon estimation by a minimizer

adds an approximate interpolation method

Parameters:

n_components – Number of diffusion map components
n_neighbors – Number of nearest neighbors for graph construction
alpha – Normalization parameter
epsilon – Kernel bandwidth; if None, it will be estimated

fit(x)¶

fit_transform(x)¶

load_state_dict(d)¶

Return type:: None

state_dict()¶

Return type:: dict[str, Any]

transform(x, ifsym=False)¶: Nystrom extension for diffusion maps.

class dymad.numerics.DMF(n_components, n_neighbors=None, alpha=1, epsilon=None)¶

Bases: object

Diffusion map with dense matrix implementation

A Knn option is added to emulate the DM class, that is a sparsified version of DM.

Also includes a KRR implementation for sanity check of other classes.

Rule of thumb comparing DM and DMF:

For eigenpairs, full DM gives accurate eigenvectors, but less accurate eigenvalues.

For KRR, full DM, i.e., without Knn, performs the best.

Whenever memory allows, use full DM.

fit(x)¶

fit_krr(X, y, ridge)¶

fit_transform(x)¶

load_state_dict(d)¶

Return type:: None

predict_krr(X)¶

state_dict()¶

Return type:: dict[str, Any]

transform(x)¶: Nystrom extension for diffusion maps.

class dymad.numerics.DimensionEstimator(data=None, Knn=None, tree=None, bracket=None, tol=0.2)¶

Bases: object

Estimate the intrinsic dimension of a point cloud using kernel method.

Based on https://doi.org/10.1016/j.acha.2015.01.001 See Fig. 6

In implementation, we analytically evaluate

d(log(S(e)))/d(log(e))

There are three operation modes:

Given data only, Knn=None: use all pairwise distances

Given data and Knn=int: use kNN distances, a KDTree will be built

Given tree and Knn=int: use kNN distances from the tree, data will not be used

Parameters:

data – Input data, shape (N, d).
Knn – Number of nearest neighbors to use. If None, use all pairwise distances.
tree – A precomputed KDTree instance. If given, data is not used.
bracket – Bracket for the scalar minimization
tol – Tolerance for biased rounding of fractional dimension

plot(N=20, fig=None, sty='b-')¶

sanity_check(K=None, ifref=True, ifnrm=True)¶

class dymad.numerics.Manifold(data, d, K=None, g=None, T=None, iforit=False, extT=None)¶

Bases: object

classmethod from_tensors(t)¶

gmls(x, Y, ret_der=False)¶

plot2d(N, scl=1)¶

plot3d(N, scl=1)¶

precompute()¶

to_tensors()¶

class dymad.numerics.ManifoldAltTree(data, d, K=None, g=None, T=None, iforit=False, extT=None, tree_data=None, tree_transform=None)¶

Bases: Manifold

The only difference from Manifold is that a KDTree on a different set of points is used to find the kNN points for GMLS.

Specifically, instead of finding kNN points of x from the original data X, we first transform x to z using a given transform function, and then find the kNN points of z from a given set of points Z (with a KDTree built on Z).

The rest is still the same as Manifold.

class dymad.numerics.ManifoldAnalytical(data, d, K=None, g=None, fT=None)¶

Bases: Manifold

precompute()¶

class dymad.numerics.VBDM(n_components, n_neighbors, Kb=None, operator=None)¶

Bases: object

Variable bandwidth diffusion map

https://www.sciencedirect.com/science/article/pii/S1063520315000020

Original MATLAB implementation by Tyrus Berry & John Harlim

This implementation

is class-based

accelerates epsilon estimation by a minimizer

adds an approximate interpolation method

It works well when the sampling is dense

Parameters:

n_components – Number of diffusion map components
n_neighbors – Number of nearest neighbors for graph construction
Kb – Number of nearest neighbors for density estimation
operator – ‘lb’ for Laplace-Beltrami, ‘kb’ for Kolmogorov backward

fit(x)¶

fit_transform(x)¶

load_state_dict(d)¶

Return type:: None

state_dict()¶

Return type:: dict[str, Any]

transform(x, ret_den=False)¶: Interpolate eigenfunctions to new x

transform_naive(x)¶

dymad.numerics.central_diff(f, x, h=1e-06, v=None)¶

Central difference.

Parameters:

f – function handle that accepts and returns numpy arrays.
x – input array.
h – step size.
v – directions for directional derivative. If None, return full Jacobian.

Returns:

derivative of f at x, possibly directional.

Return type:

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

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

dymad.numerics.complex_grid(grid)¶

Args:

grid: If a real array of 2xN, a meshgrid of complex values is formed.: If a list of two floats (a,b) and an int (N), a uniform grid is created over [-a, a]x[-b, b], each side N points. Otherwise use as is.

dymad.numerics.complex_plot(grid, sv, levels, fig=None, mode='line', lwid=2, lsty=None)¶

dymad.numerics.complex_step(f, x, h=1e-20, v=None)¶

Complex step differentiation.

Parameters:

f – function handle that accepts and returns numpy arrays.
x – input array.
h – step size.
v – directions for directional derivative. If None, return full Jacobian.

Returns:

derivative of f at x, possibly directional.

Return type:

dymad.numerics.denoise(data, *, method='savgol', axis=0, **kwargs)¶

Apply a model-independent denoising method while preserving the input array type.

Return type:: ndarray | Tensor

dymad.numerics.denoising_metrics(*, original, denoised)¶

Aggregate denoising deltas and roughness statistics across one or more signals.

Return type:: dict[str, float]

dymad.numerics.disc2cont(z, dt)¶

dymad.numerics.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.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.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.fe_step(f, t, x, dt, **kwargs)¶

Perform a single Forward Euler integration step.

Parameters:

f – Function that computes the derivative dx/dt = f(t, x, **kwargs).
t – Current time.
x – Current state vector.
dt – Time step for the integration.
**kwargs – Additional arguments to pass to the function f.

dymad.numerics.generate_coef(order, eps)¶

dymad.numerics.generate_weak_weights(dt, n_integration_points, poly_order=4, int_rule_order=4)¶

Generate weights for a weak formulation approach, returning arrays (C, D) and the number of “windows” K along the time dimension.

Steps:

Compute a length scale:

\[L = \frac{t_{N-1} - t_0}{2}\]

Generate Jacobi polynomial basis (P0) and its derivative (P1), each with alpha=1, beta=1, on a grid of size n_integration_points in [-1,1].
Construct weighting arrays for the integrals:

\[w_0 = 1 - h^2,\quad w_1 = -\frac{2h}{L}\]

Compute Newton-Cotes integration weights on that grid.
Combine everything to form:

\[C = -(P_1 w_0 + P_0 w_1) w,\quad D = P_0 w_0 w\]

Parameters:

dt (float) – Time step size.
n_integration_points (int) – Number of integration points (N) in [-1, 1].
poly_order (int) – Order of the Jacobi polynomial basis.
int_rule_order (int) – Order of Newton-Cotes integration rule (1..4).

Returns:

Weights C and D:

C: shape (poly_order, n_integration_points),
D: shape (poly_order, n_integration_points).

Return type:

Tuple[np.ndarray, np.ndarray]

dymad.numerics.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.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.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.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.rational_kernel(theta, order, eps)¶

dymad.numerics.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.rk4_step(f, t, x, dt, **kwargs)¶

Perform a single Runge-Kutta 4th order (RK4) integration step.

Parameters:

f – Function that computes the derivative dx/dt = f(t, x, **kwargs).
t – Current time.
x – Current state vector.
dt – Time step for the integration.
**kwargs – Additional arguments to pass to the function f.

dymad.numerics.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.tangent_1circle(x)¶

dymad.numerics.tangent_2torus(x, R)¶

dymad.numerics.torch_jacobian(f, x, v=None, dtype=torch.float64)¶

Jacobian using torch.func.jacobian.

Parameters:

f – function handle that accepts and returns torch tensors.
x – input array.
v – directions for directional derivative. If None, return full Jacobian.

Returns:

derivative of f at x, possibly directional.

Return type:

dymad.numerics.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.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

Modules

`complex`
`denoising`
`dm`
`gradients`
`linalg`
`manifold`
`spectrum`
`time_int`
`weak`