Linalg

Eigendecomposition

resolvent4py.linalg.eigendecomposition.arnoldi_iteration(L: LinearOperator, action: Callable[[Vec, Vec | None], Vec], krylov_dim: int, verbose: int | None = 0) → Tuple[BV, ndarray][source]

Perform the Arnoldi iteration algorithm to compute an orthonormal basis and the corresponding Hessenberg matrix for the range of the linear operator specified by L and action.

Parameters:
Returns:

tuple with an orthonormal basis for the Krylov subspace and the Hessenberg matrix

Return type:

(BV with krylov_dim columns, numpy.ndarray of size krylov_dim x krylov_dim)

resolvent4py.linalg.eigendecomposition.check_eig_convergence(action: Callable[[Vec, Vec], Vec], D: ndarray, V: BV, monitor: bool | None = False) → array[source]

Check convergence of the eigenpairs by measuring \(\lVert L v - \lambda v\rVert\) for each pair \((\lambda, v)\).

Parameters:
Returns:

Error vector (each entry is the error of the corresponding eigen pair)

Return type:

np.array

resolvent4py.linalg.eigendecomposition.eig(L: LinearOperator, action: Callable[[Vec, Vec], Vec], krylov_dim: int, n_evals: int, process_evals: Callable[[ndarray], ndarray] | None = None, verbose: int | None = 0) → Tuple[ndarray, BV][source]

Compute the eigendecomposition of the linear operator specified by L and action. For example, to compute the eigenvalues of \(L\) closest to the origin, set action = L.solve and process_evals = lambda x: 1./x.

Parameters:
Returns:

tuple with the desired eigenvalues and corresponding eigenvectors

Return type:

(numpy.ndarray of size n_evals x n_evals, SLEPc.BV with n_evals columns)

resolvent4py.linalg.eigendecomposition.match_right_and_left_eigenvectors(V: BV, W: BV, Dv: ndarray, Dw: ndarray) → Tuple[BV, BV, ndarray, ndarray][source]

Scale and sort the right and left eigenvectors and corresponding eigenvalues of an underlying operator \(L\), so that

\[W^* L V = D_v = D_w,\quad W^* V = I \in\mathbb{R}^{m\times m}.\]
Parameters:

  • V (SLEPc.BV) – m right eigenvectors

  • W (SLEPc.BV) – m left eigenvectors

  • Dv (numpy.ndarray of size m x m) – right eigenvalues

  • Dv – eigenvalues computed from the right eigendecomposition of \(L\)

  • Dw (numpy.ndarray of size m x m) – eigenvalues computed from the left eigendecomposition of \(L\). (Attention: these have already been complex conjugated.)

Returns:

tuple \((V, W, D_v, D_w)\) with the biorthogonalized eigenvectors and corresponding eigenvalues

Return type:

(SLEPc.BV with m columns, SLEPc.BV with m columns, numpy.ndarray of size m x m, numpy.ndarray of size m x m)

Randomized SVD

resolvent4py.linalg.randomized_svd.check_randomized_svd_convergence(action: Callable[[Vec, Vec], Vec], U: BV, S: ndarray, V: BV, monitor: bool | None = False) → array[source]

Check the convergence of the singular value triplets by measuring \(\lVert Av/\sigma - u\rVert\) for every triplet \((u, \sigma, v)\).

Parameters:

  • action (Callable[[PETSc.Vec, PETSc.Vec], PETSc.Vec]) – one of LinearOperator.apply() or LinearOperator.solve()

  • U (SLEPc.BV) – left singular vectors

  • D (numpy.ndarray) – diagonal 2D numpy array with the singular values

  • V (SLEPc.BV) – right singular vectors

Returns:

Error vector (each entry is the error of the corresponding singular triplet)

Return type:

np.array

resolvent4py.linalg.randomized_svd.randomized_svd(L: LinearOperator, action: Callable[[BV, BV], BV], n_rand: int, n_loops: int, n_svals: int, verbose: int | None = 0) → Tuple[BV, ndarray, BV][source]

Compute the singular value decomposition (SVD) of the linear operator specified by L and action using a randomized SVD algorithm. (See [Halko2011].) For example, with L.solve_mat we compute

\[L^{-1} = U \Sigma V^*.\]
Parameters:

  • L (LinearOperator) – instance of the LinearOperator class

  • action (Callable[[SLEPc.BV, SLEPc.BV], SLEPc.BV]) – one of LinearOperator.apply_mat() or LinearOperator.solve_mat()

  • n_rand (int) – number of random vectors

  • n_loops (int) – number of randomized svd power iterations (see [Ribeiro2020] for additional details on this parameter)

  • n_svals (int) – number of singular triplets to return

  • verbose (Optional[int], default is 0) – defines verbosity of output to terminal (useful to monitor progress during time stepping). = 0 no printout to terminal, = 1 monitor randomized SVD iterations.

Returns:

leading n_svals left singular vectors, singular values and right singular vectors

Return type:

Tuple[SLEPc.BV, np.ndarray, SLEPc.BV]

References

[Halko2011]

Halko et al., Finding Structure With Randomness: Probabilistic Algorithms For Constructing Approximate Matrix Decompositions, SIAM Review, 2011

[Ribeiro2020]

Ribeiro et al., Randomized resolvent analysis, Physical Review Fluids, 2020

Resolvent Analysis via Time Stepping

resolvent4py.linalg.resolvent_analysis_time_stepping.resolvent_analysis_rsvd_dt(L: LinearOperator, dt: float, omega: float, n_omegas: int, n_periods: int, n_rand: int, n_loops: int, n_svals: int, tol: float | None = 0.001, verbose: int | None = 0) → Tuple[BV, ndarray, BV][source]

Perform resolvent analysis using randomized linear algebra and time stepping. In particular, it can be shown that

\[x_\omega e^{i\omega t} = \left(i\omega I - A\right)^{-1}f_\omega e^{i\omega t} \to \int_0^t e^{A(t-\tau)}f_\omega e^{i\omega\tau} d\tau\]

for sufficiently long time \(t \gg 1\) and for any complex-valued forcing \(f(t) = f_\omega e^{i\omega t}\) (assuming that \(A\) is stable). Computing the integral on the right-hand side can be done by integrating

\[\frac{d}{dt}x(t) = Ax(t) + f(t)\]

forward in time with initial condition \(x(0) = 0\). Thus, the action of the resolvent operator \(R(i\omega) = \left(i\omega I - A\right)^{-1}\) on a vector \(f_\omega\) can be computed by time-stepping the ODE above. For now, time integration is performed explicitly via the Adams-Bashforth scheme. (See [Martini2021] and [Farghadan2025] for more details.)

Note

This function is meant for systems whose dimension is so large that linear systems of the form \((i\omega I - A)x_\omega = f_\omega\) cannot be solved easily. Typically, this happens when \(\text{dim}(x_\omega) \sim O(10^7)\) or larger. If you have a “small enough” system, then we highly recommend using randomized_svd.randomized_svd() instead for the singular value decomposition of the resolvent.

Parameters:

  • L (LinearOperator) – linear operator representing \(A\)

  • dt (float) – target time step \(\Delta t\) to integrate the ODE

  • omega (float) – fundamental frequency \(\omega\) of the forcing

  • n_omegas (int) – number of integer harmonics of \(\omega\) to resolve

  • n_periods (int) – number of periods \(T = 2\pi/\omega\) to integrate the ODE through. (This number should be large enough that we can expect transients to have decayed.)

  • n_rand (int) – number of random forcing vectors for randomized SVD

  • n_loops (int) – number of power iterations for randomized SVD

  • n_svals (int) – number of singular values/vectors to output

  • tol (Optional[float], default is \(10^{-3}\)) – integrate the ODE forward for n_periods or until \(\lVert x(kT) - x((k-1)T) \rVert < \mathrm{tol}\).

  • verbose (Optional[int], default is 0) – defines verbosity of output to terminal (useful to monitor progress during time stepping). = 0 no printout to terminal, = 1 monitor randomized SVD iterations, = 2 monitor randomized SVD iterations and time-stepping progress.

Returns:

left singular vectors, singular values, and right singular vectors of the resolvent operators \(R(i\omega)\) evaluated at frequencies \(\Omega = \omega\{0, 1, 2, \ldots, n_{\omega}\}\) if the linear operator \(A\) is real-valued; otherwise at frequencies \(\Omega = \omega\{0, 1, 2, \ldots, n_{\omega}, -n_{\omega},-(n_{\omega}-1) \ldots, -1\}\)

Return type:

Tuple[List[SLEPc.BV], List[np.ndarray], List[SLEPc.BV]]

References

[Martini2021]

Martini et al., Efficient computation of global resolvent modes, Journal of Fluid Mechanics, 2021

[Farghadan2025]

Farghadan et al., Scalable resolvent analysis for three-dimensional flows, Journal of Computational Physics, 2025