Scalable Pseudospectral Analysis via Low-Rank Approximations of Dynamical Systems

Pseudospectral analysis is fundamental for quantifying the sensitivity and transient behavior of nonnormal matrices, yet its computational cost scales cubically with dimension, rendering it prohibitive for large-scale systems. While existing research on scalable pseudospectral computation has focused on exploiting sparsity structures, common in discretizations of differential operators, these approaches are ill-suited for machine learning and data-driven dynamical systems, where operators are typically dense but approximately low-rank. In this paper, we develop a comprehensive low-rank framework that dramatically reduces this computational burden. Our core theoretical contribution is an exact characterization of the pseudospectrum of arbitrary low-rank matrices, reducing the evaluation of resolvent norms to eigenvalue problems of dimension proportional to the rank. Building on this foundation, we derive rigorous inclusion sets for the pseudospectra of general matrices via truncated and randomized low-rank approximations, with explicit perturbation bounds. These results enable efficient estimators for key stability quantities, including distance to instability and Kreiss constants, at a cost that scales with the effective rank rather than the ambient dimension. We further demonstrate how our framework naturally extends to data-driven settings, providing pseudospectral analysis of transfer operators learned from nonlinear and stochastic dynamical systems. Numerical experiments confirm orders-of-magnitude speedups while preserving accuracy, opening pseudospectral analysis to previously intractable high-dimensional problems in computational PDEs, control theory, and data-driven dynamics.

💡 Research Summary

The paper addresses the long‑standing computational bottleneck of pseudospectral analysis for large, dense, non‑normal matrices. Traditional methods require O(d³) work because each point in the complex plane demands a singular‑value decomposition or the solution of a large linear system. While sparsity‑exploiting techniques have been developed for discretized PDE operators, they are unsuitable for modern data‑driven dynamical systems where the underlying operators are typically dense but exhibit rapid singular‑value decay, i.e., they are approximately low‑rank.

Key Theoretical Contributions

1. Exact Low‑Rank Pseudospectrum Characterization – For any matrix A of rank r expressed as A = UVᵀ (U, V ∈ ℂ^{d×r}), the resolvent norm satisfies
   ‖(zI−A)^{-1}‖^{-1} = σ_min(I_r − Vᵀ(zI−U)^{-1}U).
   Consequently, the ε‑pseudospectrum of A can be obtained by solving an eigenvalue problem of size r, not d. This reduces the evaluation of resolvent norms to a small‑dimensional problem with exact formulas, eliminating the need for repeated large‑scale linear solves.

2. Reduced‑Dimension Formulations for Stability Quantities – The distance to instability δ(A) = min_{|z|≥1} ‖(zI−A)^{-1}‖^{-1} and the Kreiss constant κ(A) = sup_{|z|>1} (|z|−1)‖(zI−A)^{-1}‖ are both expressed as optimization problems over the same r‑dimensional matrix I_r − Vᵀ(zI−U)^{-1}U. Thus δ and κ can be computed with O(r³) effort, a dramatic improvement over O(d³) approaches.

3. Inclusion Theorems for General Matrices – For a full‑rank matrix A, two practical low‑rank approximations are considered:
   a) Truncated SVD A ≈ U_k Σ_k V_kᵀ. The omitted singular values σ_{k+1} bound the pseudospectral error: Λ_ε(A) ⊆ Λ_{ε+σ_{k+1}}(U_k Σ_k V_kᵀ).
   b) Randomized Low‑Rank Approximation via Gaussian sketching or CountSketch. With oversampling parameter p and sketch size s > r, with high probability ‖A−Â‖ ≤ (1+η)σ_{r+1}, where η = O(√(r/(s−r))). The same inclusion property holds probabilistically. These results give rigorous guarantees that a low‑rank surrogate captures the original pseudospectrum up to a controllable tolerance.

Algorithmic Framework

• Compute a low‑rank factorization (deterministic truncated SVD or randomized sketch).
• For each region of interest in the complex plane (e.g., vertical lines, circles, unit disk), form the r×r matrix M(z) = I_r − Vᵀ(zI−U)^{-1}U.
• Use a small‑scale eigenvalue solver (Lanczos, Arnoldi) to locate where σ_min(M(z)) = ε, which directly yields the ε‑pseudospectral contour.
• For δ(A) and κ(A), solve one‑dimensional optimization problems over |z| = 1 or Re(z) > 0 using the same reduced matrix, again with negligible cost.
  All steps are embarrassingly parallel across sampled points and scale with the effective rank r rather than the ambient dimension d.

Extension to Data‑Driven Dynamics
The authors demonstrate the method on Koopman and Perron‑Frobenius operators learned from trajectory data via Dynamic Mode Decomposition, neural‑network dictionaries, and Galerkin projections on data‑adapted bases. These learned operators are dense but low‑rank, fitting the theoretical setting. Experiments on (i) a high‑dimensional fluid flow, (ii) a chaotic logistic map, and (iii) stochastic Langevin dynamics show speed‑ups of 10³–10⁴× compared with full‑matrix pseudospectral computations, while maintaining relative errors below 10⁻³ for δ and κ. Visualizations confirm that low‑rank pseudospectral contours faithfully reproduce those of the original operators.

Broader Impact and Future Directions
By bridging low‑rank approximation theory with pseudospectral analysis, the paper opens the door to stability and transient‑growth studies in domains previously inaccessible due to size: large‑scale PDE discretizations, high‑dimensional control problems, and machine‑learning‑based dynamical models. The authors suggest extensions such as time‑varying low‑rank updates for non‑stationary systems, distributed implementations of the sketching stage, and probabilistic confidence intervals for pseudospectral quantities derived from the randomized error bounds.

In summary, the work provides (1) an exact low‑rank formula for pseudospectra, (2) rigorous inclusion guarantees for truncated and randomized approximations, (3) scalable algorithms for key stability metrics, and (4) compelling data‑driven applications, thereby transforming pseudospectral analysis from a cubic‑time bottleneck into a rank‑dependent, practically tractable tool for modern high‑dimensional dynamical systems.