PaperlayerView exact revision on arXiv

Source-linked AI summary

Solving a sequence of eigenvalue problems using snapshot and row subsampling

Francesco Hrobat, Yuji Nakatsukasa

arXiv:2610.19326v1math.NAcs.NA

TL;DR

Sequences of largest- or smallest-eigenvalue problems arise in many applications, and the paper develops a snapshot-based approach with row subsampling. Its analysis gives exponential convergence results under analytic-extension assumptions, while the selected rows remain reusable when the snapshot matrix is unchanged and the initial selection is good.

  • Problem

    Solving sequences of largest- or smallest-eigenvalue problems is important in applications, and reduced-order methods seek solutions in a lower-dimensional subspace.

    Source: §1 · p. 1

  • Method

    The algorithm combines snapshots with sketched Rayleigh–Ritz and row subsampling, selecting rows from the snapshot matrix for the reduced computation.

    Source: §1 · pp. 2–3

  • Results

    With analytic A(p), snapshot singular values decay exponentially, and eigenvector-approximation convergence is asymptotically exponential for any distinct snapshot distribution; a well-chosen row subset remains valid while X is unchanged.

    Source: §1 · p. 3, §3.1 · p. 8, §3.2 · p. 14

  • Takeaways & Limitations

    Once a suitable row subset is chosen, the algorithm can compute SA(p) without accessing the full matrix A(p).

    Source: §1 · p. 3

    The theoretical error analysis focuses on a simple isolated eigenvalue, with extension to multiple eigenvalues under Assumption 2.4.

    Source: §3.1 · p. 6

Abstract

from arXiv · show

A sequence of eigenvalue problems requires solving $A_i\vec{x}_i= \lambda_i\vec{x}_i$ where $i \in \{1, \dots, N\}$, $N \gg 1$ and (usually) $\lambda_i$ is the largest or smallest eigenvalue of $A_i$. Often, the matrices $A_i$ come from a discretization of a smooth map depending on a parameter $p \rightarrow A(p)$, that is $A_i = A(p_i)$, so the matrices $A_i$ change gradually along the sequence. Under this latter assumption, we introduce an efficient algorithm for solving sequences of eigenvalue problems. Our algorithm combines the ideas of reduced order model and row subsampling, as done in [E. Jones and Y. Nakatsukasa, arXiv, 2025] for linear systems. In addition, in our approach, the same set of rows is guaranteed to work well independently of the value of the parameter $p$ or the matrix $A(p)$ as long as the snapshot matrix is unchanged and given that the row subset was chosen well at the beginning. Furthermore, we implement an update strategy of the snapshot matrix. Along these results, we present theoretical analysis of the snapshot matrix for eigenvalue problems, showing an exponential decay of the singular values, and we give bounds on the error on eigenvectors for unseen value of the parameter. This analysis generalizes that in [M. Guido, D. Kressner, and P. Ricci, J. Sci. Comput., 2024] and shows that, regardless of the distribution of the snapshot points, the convergence is always (asymptotically) exponential in the number of snapshot points if the function $A(p)$ admits an analytical extension. We further generalize our algorithm for computing largest or smallest singular values and relative singular vectors and to generalized eigenvalue problems. We conclude with numerical results showing the efficiency of our approach.

1. Introduction.

The introduction frames sequences of largest- or smallest-eigenvalue problems and combines snapshots with row-subsampled Rayleigh–Ritz to reduce their solution cost. It also previews theoretical guarantees for snapshot quality and eigenvector approximation.

  • For new parameter values, the eigenvector is sought in the span of snapshot eigenvectors, reducing the search dimension to r ≪ n,N.Source: §1 · p. 2
  • POD uses leading left singular vectors of the snapshot matrix as a reduced projection basis, optionally truncating according to singular values.Source: §1 · p. 2
  • The method combines snapshots with sketched Rayleigh–Ritz, replacing the projected matrix with one formed using a row-subsampling matrix S.Source: §1 · p. 2
  • After selecting rows, the method can compute SA(p) without accessing all of A(p), and reduces the least-squares problem from n×r to s×r.Source: §1 · p. 3
  • The introduction previews exponential singular-value decay for analytic A(p) and eigenvector-error bounds for unseen parameters, with asymptotically exponential convergence for distinct snapshot distributions.Source: §1 · p. 3
  • The paper distinguishes time-dependent problems, whose matrices are accessed sequentially, from parameter-dependent problems, whose matrices can be accessed independently.Source: §1 · p. 3

2. Snapshot matrix for eigenvalue problems.

Under analytic-extension and eigenvalue-separation assumptions, the snapshot matrix has exponentially decaying singular values. The analysis also treats multiple-eigenvalue eigenspaces and relates the decay to reduced-basis approximation.

  • The analysis assumes an analytic extension of A(p) and a simple, separated eigenvalue; it establishes a continuous eigenvector map that is analytic on the domain.Source: §2 · p. 3
  • For an analytic eigenvector map, the snapshot matrix singular values decay exponentially, supporting the use of a reduced basis.Source: §2 · p. 4
  • The singular-value result justifies POD, while a spectral gap is fundamental to eigenvector analyticity; when an eigenvalue is not separated, the analysis can instead use a larger eigenvalue set.Source: §2 · p. 4
  • For multiple eigenvalues, a uniformly enclosed, separated eigenvalue cluster admits an analytic basis, extending the snapshot result with mr replacing r.Source: §2 · pp. 4–5
  • In the 1000×1000 symmetric random-matrix example, the estimated decay rate ρ≈3.1 reflects observed behavior for Chebyshev and equispaced snapshots.Source: Figure 1 · p. 5

3. Theoretical results.

The analysis bounds how well snapshot spaces capture eigenvectors at unseen parameter values, then examines Rayleigh–Ritz with row subsampling. Under the stated assumptions, the bounds cover different snapshot distributions, and a well-chosen row subset can be reused while the snapshot matrix is unchanged.

3.1. Theoretical results on interpolation and extrapolation.

  • The projection-error analysis assumes a simple isolated eigenvalue; a separated group of eigenvalues can instead be treated as a larger eigenspace under the stated assumptions.A spectral gap is fundamental to the analyticity of the eigenvector.Source: §3.1 · p. 6
  • Theorems 3.1 and 3.2 bound projection error for equispaced and Chebyshev snapshots, respectively, with results covering interpolation and extrapolation.The paper relates equispaced snapshots to time-dependent problems and Chebyshev snapshots to parameter-dependent problems; it reports no observed Runge phenomenon for the projection error.Source: §3.1 · pp. 6–7
  • For snapshot points whose empirical measures converge weak-star to a unit positive measure, the eigenvector projection error has an asymptotically exponential rate under the stated analyticity assumptions.The rate is described through the limiting measure’s logarithmic potential; the analysis treats new parameter values inside or outside [−1, 1].Source: §3.1 · pp. 8–9

3.2. Theoretical results on subsampling for Rayleigh-Ritz.

  • The subsampling guarantee is conditional: the row subset must be chosen well, X must be suitably conditioned, and the reduced basis should be nearly invariant under A(p).Using leading left singular vectors of X through POD is given as an example for controlling its condition number.Source: §3.2 · p. 13
  • For Hermitian A(p), the reduced matrix formed with X need not be Hermitian, so its computed eigenvalues are not always real.The paper notes that orthogonalizing X and using an orthonormal basis Q changes this issue.Source: §3.2 · p. 14
  • When the snapshot matrix is unchanged, a row subset chosen well at the outset can be reused across parameter values, and the subsampling bound depends on X rather than A(p).The paper attributes error in that setting to snapshot-matrix performance as p varies, rather than to subsampling.Source: §3.2 · p. 14

4. Main algorithm.

SubApSnapEig builds a snapshot basis and row subset offline, then uses them to approximate eigenpairs at new parameter values. The online phase can reuse the same snapshot and rows, with periodic or residual-triggered updates.

4. Main algorithm.

  • The offline phase computes eigenpairs at snapshot parameters and collects their eigenvectors into a snapshot matrix X.X may be replaced by leading left singular vectors Q, with truncation based on singular values.Source: §4 · p. 14
  • The algorithm selects rows from X or Q to form S, or can use a more conservative strategy that samples from a matrix with additional columns.The conservative option can support a cheap residual estimate but requires access to A(p) and products A(p)X; nearest-snapshot row selection is another option that avoids accessing A(p) at new parameters.Source: §4 · p. 15
  • The same X and S can serve many new parameters, but the algorithm advises updating them periodically or when the estimated residual is large.When updating X, it refines the latest computed eigenvectors with Krylov iterations and adds them to X.Source: §4 · pp. 15–16
  • With s = O(mr), the stated online costs include O(n(mr)^2) for the sampled product and O((mr)^3) each for least squares and eigenpair computation.The passage also gives offline-phase cost components, including O(n^3r) for computing X.Source: §4 · p. 15
  • For a new parameter, the online phase solves a subsampled least-squares problem and computes the corresponding eigenpairs of the reduced matrix.SX can be precomputed, and SA(p*)X can be formed from sampled rows, potentially without accessing the full matrix.Source: §4 · pp. 15–16

5. Implementation and residual estimation.

This section introduces implementation details for the algorithm.

5. Implementation and residual estimation.

  • Section 5 discusses implementation details for the algorithm.Source: §5 · p. 16

Computation of SA(p∗)X.

The sampled product is formed by sampling rows of A(p*) before multiplying by X; orthogonalization can improve numerical stability when snapshots are redundant.

Computation of SA(p∗)X.

  • Forming SA(p*)X as (SA(p*))X can avoid accessing all entries of A(p) when single-entry computation is expensive.The passage identifies entry-wise access and costly entries as a setting where sampling first is more efficient.Source: Computation of SA(p∗)X · p. 16
  • Orthogonalizing X with a rank-revealing QR or SVD is particularly useful when mr is large and X is numerically rank deficient.The passage says orthogonalization improves numerical stability and conditioning of least-squares problems.Source: Computation of SA(p∗)X · p. 16

Update of X and S.

The algorithm updates its snapshots and row subset to incorporate new eigenvector information and maintain performance across parameter values.

Update of X and S.

  • Update X after a prescribed number of parameters or when the estimated residual is large.The passage highlights time-dependent problems as a setting where later parameters may move farther from initial snapshot points.Source: Update of X and S · p. 16
  • New information enters X by refining computed eigenvectors with a few Krylov iterations, initialized from the previously computed eigenvector.The passage notes that computed eigenvectors otherwise remain in the span of X.Source: Update of X and S · pp. 16–17
  • The row subset S is chosen using X or its left singular vectors Q, without computing A(p)X for the update.Source: Update of X and S · p. 17

Subsampling strategies.

The section surveys row-selection options and develops a leverage-score-based residual estimate to help decide when the snapshot matrix or selected rows should be updated.

Subsampling strategies.

  • Row selection can use deterministic methods such as LUPP, CPQR, and Osinsky’s approach, or randomized leverage scores, determinantal point processes, and adaptive randomized pivoting.Subsampling [X, A(p)X] requires accessing A(p) and selecting more rows, but gives stronger theoretical bounds.Source: Subsampling strategies · p. 17

5.1. Residual estimate and change of S.

  • With leverage-score sampling and oversampling, a subspace-embedding property supports a cheap residual estimate for the least-squares solution.The estimate is developed for rows scaled by 1/√(s p_i), where p_i is the probability of selecting row i.Source: §5.1 · p. 17
  • When subsampling from [X, AX], the residual estimate requires only an O(sr) computation using quantities already formed for the least-squares problem.The estimate uses (SA(p)X)y(p).Source: §5.1 · p. 17
  • Because computing ∥A(p)∥ requires accessing the entire matrix, the method reuses the norm saved at the nearest snapshot point for a new parameter.This norm reuse is part of the practical relative-error calculation.Source: §5.1 · p. 18
  • In practice, subsampling from X alone gives only a rough residual estimate, while non-oversampled LUPP or CPQR provides no residual information.A large estimate suggests updating X or S; without oversampling, the residual cannot be estimated this way.Source: §5.1 · p. 18

6. Related algorithms.

The section situates the approach among methods that reuse information from earlier problems and distinguishes its eigenvalue setting from DEIM’s vector approximation setting.

6. Related algorithms.

  • Prior approaches that exploit previously computed eigenpairs include Krylov recycling, rational Krylov methods, Taylor and Chebyshev expansions, and snapshot-based model reduction.DEIM is identified as a related method that also combines snapshots with subsampling.Source: §6 · p. 18
  • DEIM approximates a parameterized vector in a low-dimensional snapshot subspace and solves a subsampled least-squares problem.Its setup assumes individual entries of the target vector can be evaluated relatively quickly.Source: §6 · p. 18
  • The DEIM entry-access assumption is usually not satisfied when the target vector is an eigenvector of A(p).The paper contrasts this with eigenvalue problems, where the eigenvector is the quantity to approximate.Source: §6 · p. 18

7. SVD and GEP.

The paper extends its snapshot-and-subsampling approach to singular-value problems and generalized eigenvalue problems, with different costs and row-selection choices.

7. SVD and GEP.

  • For SVD, the method uses separate snapshot matrices for right and left singular vectors, then solves the resulting least-squares problem.Rows can be selected from the left-vector snapshot matrix without accessing the entire A(p).Source: §7 · p. 18
  • An alternative SVD strategy sketches the approximation of Xℓ†A(p), equivalent to Rayleigh-Ritz on A(p)(A(p))*, but its least-squares problem has dimension s × n.Computing left singular vectors costs O(nr²) for this strategy, versus O(r³) for the first method.Source: §7 · p. 19
  • For generalized eigenvalue problems, the paper sets up a sequence of matrix pairs (A(p), B(p)) and derives a corresponding least-squares problem.The supplied passages introduce the formulation but do not provide further details of that least-squares system.Source: §7 · p. 19
  • For generalized problems, rows may be selected from B(p)X at snapshot points and reused from the nearest snapshot for a new parameter.In time-dependent problems, updating these rows as snapshot points move farther away requires accessing the entire B(p).Source: §7 · p. 19

8. Numerical Experiments.

The experiments compare SubApSnapEig with MOR across differential operators, PageRank, a Gaussian kernel, graph centrality, and Transformer attention. Results show substantial online speedups in several examples, with accuracy trade-offs and snapshot updates affecting performance.

8. Numerical Experiments.

  • The experiments compare SubApSnapEig using LU or leverage-score row extraction against MOR, which projects A(p) online without row subsampling.Initial snapshot eigenpairs are computed with MATLAB eigs; the snapshot matrix is updated only in the time-dependent graph example.Source: §8 · p. 19

8.1. Example 1: Exponential Differential Operator.

  • 20× faster: SubApSnapEig outpaces MOR on the exponential differential operator, while MOR’s eigenvalue and eigenvector errors are about four orders of magnitude smaller.The SubApSnapEig errors are of order 10^-4, and the error does not significantly dip around snapshot points.Source: §8.1 · p. 20, Figure 4 · p. 21

8.2. Example 2: Differential operator with multiple eigenvalues.

  • 25× faster: SubApSnapEig is faster than MOR for the multiple-eigenvalue differential operator, although its errors are two orders of magnitude larger.Here the errors dip around snapshot points, and the LU variant is particularly faster because each row has more nonzeros.Source: §8.2 · p. 21, Figure 5 · p. 22

8.3. Example 3:

  • On PageRank, the methods have nearly equal errors; LU SubApSnapEig has a two-orders-smaller eigenvalue error, while MOR is only 2× slower.The authors report that the accuracy behavior is unclear.Source: §8.3 · p. 23, Figure 6 · p. 23

8.4. Example 4: Gaussian Kernel.

  • 250× faster: SubApSnapEig reduces online time for the Gaussian kernel, while MOR’s eigenvalue error is several orders smaller over the first half of the parameter domain.The other algorithms’ errors are described as pretty similar.Source: §8.4 · p. 24, Table 4 · p. 24

8.5. Example 5: Graph Eigenvector Centrality.

  • For graph eigenvector centrality, snapshot updates reduce upcoming parameter errors, and not updating X produces worse accuracy; SubApSnapEig is approximately 6× faster online.More frequent snapshot updates make the algorithm less efficient relative to MOR.Source: §8.5 · p. 25, Figure 8 · p. 26, Table 5 · p. 26

8.6. Example 6: Transformer.

  • On Transformer attention, LU SubApSnapEig overestimates the leading singular value by a factor of 2, while the other methods achieve 10^-2 singular-value and 10^-1 singular-vector accuracy.SubApSnapEig is around 20× faster than MOR; the reported run uses 50 initial snapshots without updating the snapshot matrix.Source: Figure 8 · p. 26, §8.6 · p. 27, Table 6 · p. 27
Loading 2610.19326v1…