An Improved Incremental Singular Value Decomposition and New Error Bounds

Technical Summary: An Improved Incremental Singular Value Decomposition and New Error Bounds

Problem Statement
The incremental Singular Value Decomposition (SVD), originally proposed by Brand [2, 3] and extended by Fareed et al. [10] to weighted inner products, is a standard procedure for low-rank approximation in streaming data scenarios where the full matrix cannot be stored. The algorithm updates a truncated SVD ($U \approx Q\Sigma R^\top$) as new columns arrive. However, in floating-point arithmetic, every update multiplies the running singular basis $Q$ by a small orthogonal factor. While individual multiplications preserve orthogonality up to rounding errors, these errors accumulate over the stream length $n$. Eventually, $Q$ loses its numerical orthogonality, causing the singular values to become inaccurate and the stopping criteria to become unreliable.

The standard remedy is periodic reorthogonalization (e.g., via modified Gram–Schmidt), triggered when orthogonality degrades. However, the frequency required to ensure a prescribed precision level remained an open question. As noted by Brand [3, p. 28], it was unclear how often reorthogonalization is necessary to guarantee a specific overall level of numerical precision, particularly in weighted settings (e.g., finite-element Proper Orthogonal Decomposition), where reorthogonalization is computationally expensive.

Methodology
The authors propose a restructuring of the incremental SVD algorithm that delays updates to the running singular factors. The core methodology relies on two observations:

1. Rank-non-increasing dominance: In low-rank applications, most incoming columns lie within the current subspace (residual < $\text{tol}$). These "rank-non-increasing" columns do not expand the subspace.
2. Batching: The cumulative effect of $s$ consecutive rank-non-increasing columns can be reproduced by a single thin SVD of a short and wide buffer matrix containing the projection coefficients, rather than performing $s$ separate SVDs of bordered matrices.

The proposed algorithm (Algorithm 5) maintains a buffer for projection coefficients of rank-non-increasing columns. The running factor $Q$ remains unchanged during these steps. When a "rank-increasing" column arrives (residual $\ge \text{tol}$), the algorithm performs the following:

• It empties the buffer by computing a single thin SVD of the accumulated coefficients and updating $Q$ and $R$ in a batch.
• It processes the new column and may apply a second stopping test to the smallest singular value of the local bordered matrix to discard negligible rank.
• It applies classical Gram–Schmidt with two passes (CGS-2) only at these rank-increasing events to restore the $W$-orthogonality of the appended column.

This restructuring reduces the number of large orthogonal multiplications applied to $Q$ from $O(n)$ (the stream length) to $O(r + J_{E2})$, where $r$ is the numerical rank and $J_{E2}$ is the number of rank-increasing events with truncation.

Main Contributions
The work provides two primary theoretical advances and one practical algorithmic improvement:

1. Stream-length-independent orthogonality bound:
   The authors prove (Theorem 4.1) that under a standard probabilistic model for rounding errors, the loss of $W$-orthogonality of the computed left factor $Q$ is bounded by:
   $$\|I - Q^\top W Q\|_2 \leq c_1 (r + J_{E2}) r \sqrt{m} \log(2/\delta) \varepsilon_{\text{mach}} + O(\varepsilon_{\text{mach}}^2)$$
   Crucially, this bound is independent of the stream length $n$. It depends on $m$ (column dimension) only via a factor of $\sqrt{m}$. This resolves Brand's open question: reorthogonalization is not required at every step; the proposed batching strategy ensures that orthogonality is maintained uniformly in $n$ at machine precision.

2. Sharp forward error bound:
   The work establishes a new forward error bound for the factorization in the $W$-weighted operator norm (Theorem 4.5):
   $$\|U - Q\Sigma R^\top\|_W \leq \sqrt{n} \, \text{tol} + (\text{rounding error terms, uniform in } n)$$
   This sharpens the existing bound by Fareed and Singler [9], which was $n \, \text{tol}$. The authors demonstrate via a constructive example that the $\sqrt{n}$ rate is achievable, proving that the previous bound was too loose by a factor of $\sqrt{n}$.

3. Algorithmic efficiency:
   By reducing the number of expensive reorthogonalization steps and batching updates, the proposed algorithm significantly reduces computational cost.

Results
Numerical experiments were conducted on a 2D parabolic test problem (both unweighted and weighted finite-element settings) and synthetic data.

• Orthogonality: While existing algorithms (Brand's direct update and Fareed's reorthogonalized update) showed orthogonality degradation growing with $n$ (drifting up to $10^{-3}$ in weighted cases), the proposed algorithm maintained orthogonality at machine precision ($\sim 10^{-14}$) for stream lengths up to $n=10,000$ and column dimensions up to $m=100,000$.
• Error bounds: Experiments confirmed the $\sqrt{n} \, \text{tol}$ stopping rate. In a constructive example, the measured error matched the $\sqrt{n} \, \text{tol}$ prediction, whereas the Fareed–Singler bound ($n \, \text{tol}$) was too loose by a factor growing with $\sqrt{n}$.
• Performance: The proposed algorithm was 4.5$\times$ to 34$\times$ faster than its closest competitors. In the weighted case, the existing algorithm spent 94% of its time on reorthogonalization, while the proposed algorithm reduced these costs by three orders of magnitude.

Significance
The work claims to clarify the open question regarding the frequency of reorthogonalization in incremental SVD. By restructuring the algorithm to batch rank-non-increasing updates, the following is achieved:

1. Uniform orthogonality: The running basis remains numerically orthogonal regardless of stream length, rendering heuristic, periodic reorthogonalization checks—which lack a priori guarantees—unnecessary.
2. Tighter error analysis: It provides a mathematically rigorous, sharp error bound ($\sqrt{n} \, \text{tol}$), improving the state of the art by a factor of $\sqrt{n}$.
3. Practical acceleration: The method offers significant speed advantages in applications dominated by reorthogonalization costs, such as finite-element POD, without compromising accuracy.

The authors note that the algorithm has already been applied in preliminary versions to time-fractional PDEs, integro-differential equations, and PDE-constrained optimization, relying on the now-established orthogonality and error guarantees. Future work includes tightening the analysis of orthogonality drift (currently conservative by $\sim 10^3$) and optimizing the projection step, which currently dominates runtime.