Frequently in online or data-driven applications, new information becomes available in the form of a stream (Syamantak et al. [KS24]). Because processing the data for analysis or inference often involves matrix factorizations like the SVD or an eigendecomposition, we develop new updating methods. As repeatedly refactoring a large matrix is expensive, we propose low-rank updates to a previous factorization. Thus we consider the model

\begin{equation*} \Abar = A + C W^T, \end{equation*}

where the previous data \(A\) is \(m \times n\), and the updates \(C\) and \(W\) are \(m \times t\) and \(n \times t\). Simply computing \(\Abar \) costs \(m n t \) multiplications, which we therefore regard as an optimal complexity. For some well known factorizations, efficient updating methods are known. For instance, the methods of Gill et al. [GGMS74] update the Cholesky or LDL\(\null ^\textrm {T}\) factorization in \(\mathcal {O}(mnt) \) flops, while Golub and Van Loan [GV13, Sec 12.5] describe a method for updating the QR factorization. Also, Bunch et al. [BNS78] describe an algorithm for updating the eigendecomposition, while Brand [Bra06] and Moonen et al. [MVDV92] develop methods for the SVD. The complexity of the latter methods also scales as \(\mathcal {O}(mnt)\), but updating an eigendecomposition or SVD typically involves iterative nonlinear equation solves. In SVD computation, the first step is to reduce the matrix to (upper) bidiagonal form before computing the singular values of the bidiagonal (e.g., implemented in LAPACK’s bdsqr and gebrd [ABB+99]). Since the bidiagonalization and SVD are closely related, it is not surprising that attempts have been made to replace the SVD with the bidiagonalization for low-rank matrix approximations (Simon and Zha [SZ00]). Even though the SVD guarantees the best low-rank approximation, the bidiagonalization can be computed with a predetermined number of orthogonal transformations, making it computationally attractive.

The most stable method for computing a bidiagonal factorization uses sequences of orthogonal Householder reflectors. Because this requires and overwrites the matrix elements in memory, it is best suited for dense systems. Its complexity scales as \(\mathcal {O}(mn^2)\) flops and \(mn\) memory and it is therefore limited to small or medium problems. A second approach accesses the data only via matrix-vector products within a short two-vector recursion. This Golub-Kahan bidiagonalization (GKB) produces a partial bidiagonalization after \(k\) iterations. For a general \(A\) the GKB complexity is \(\mathcal {O}(kmn)\). When the data is sparse or otherwise structured, GKB can exploit the structure with potentially much fewer flops (but without further modifications suffers from rapid loss of orthogonality). Our algorithms reuse an existing bidiagonal factorization \(A = Q B P^T\) to compute the next factorization

\begin{equation*} \Qbar \Bbar \Pbar ^T = Q B P^T + C W^T \end{equation*}

at reduced cost. To exploit previous information fully, we develop sparsity-exploiting bidiagonalization algorithms. One method is gk-bidiag, which we compare to LAPACK’s bidiagonalization functions (Table 1). We also propose methods such as a compact representation of products of Householder matrices combined with the GKB iteration.

References

• [ABB+99]  Edward Anderson, et al. LAPCK users’ guide. SIAM, 1999

• [BNS78]  James R Bunch, Christopher P Nielsen, and Danny C Sorensen. Rank-one modification of the symmetric eigenproblem. Numerische Mathematik, 31(1):31–48, 1978.

• [Bra06]  Matthew Brand. Fast low-rank modifications of the thin singular value decomposition. Linear Algebra and its Applications, 415(1):20–30, 2006.

• [DH11]  Timothy A Davis and Yifan Hu. The University of Florida sparse matrix collection. ACM Trans. Math. Softw. (TOMS), 38(1):1–25, 2011.

• [GGMS74]  Philip E Gill, Gene H Golub, Walter Murray, and Michael A Saunders. Methods for modifying matrix factorizations. Mathematics of Computation, 28(126):505–535, 1974.

• [GV13]  G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins Studies in the Mathematical Sciences. The Johns Hopkins University Press, Baltimore, fourth edition, 2013.

• [KS24]  Syamantak Kumar and Purnamrita Sarkar. Streaming PCA for Markovian data. Advances in Neural Information Processing Systems (NeurlIPS), 36, 2024.

• [MVDV92]  Marc Moonen, Paul Van Dooren, and Joos Vandewalle. A singular value decomposition updating algorithm for subspace tracking. SIAM Journal on Matrix Analysis and Applications, 13(4):1015–1038, 1992.

• [SZ00]  Horst D Simon and Hongyuan Zha. Low-rank matrix approximation using the Lanczos bidiagonalization process with applications. SIAM Journal on Scientific Computing, 21(6):2257–2274, 2000.