Subspace-Orbit Randomized SVD

• Subspace-Orbit Randomized SVD is a matrix decomposition method that produces efficient low-rank approximations by leveraging two-sided subspace projections and random sampling.
• It achieves tight spectral and Frobenius error bounds with minimal data passes, ensuring both precision and scalability for handling large, high-dimensional matrices.
• The algorithm utilizes iterative power steps and orthonormal basis construction, making it highly compatible with modern multicore and GPU architectures for fast computations.

The Subspace-Orbit Randomized Singular Value Decomposition (SOR-SVD) is a matrix decomposition technique designed to produce efficient and accurate low-rank approximations of large, dense data matrices. Utilizing random sampling and subspace iteration principles, SOR-SVD operates with minimal passes on the data and achieves computational complexity suitable for large-scale modern computing platforms, including multicore and GPU architectures. SOR-SVD distinguishes itself from prior randomized SVD methods through a two-sided subspace projection approach and provably tight spectral and Frobenius error bounds, rendering it particularly effective for high-dimensional data applications where rapid, reliable matrix approximation is required (Kaloorazi et al., 2018).

1. Mathematical Formulation and Setup

Given a real-valued matrix $A \in \mathbb{R}^{m \times n}$ with $m \geq n$ and numerical rank $k \ll \min(m, n)$, the SOR-SVD algorithm constructs a rank-$k$ approximation $\hat{A}$ through randomized linear algebraic techniques. The procedure begins by selecting an oversampling parameter $p \geq 2$—typically in the range of 5–10—and sets $\ell = k + p$. A Gaussian test matrix $\Omega \in \mathbb{R}^{n \times \ell}$ with entries drawn i.i.d. from $\mathcal{N}(0, 1)$ is used to probe both the range of $A$ and its adjoint. This double sampling strategy is a key distinguishing feature (Kaloorazi et al., 2018).

2. Subspace-Orbit Sampling and Power Iteration

The core of SOR-SVD is its two-sided subspace-orbit sampling procedure. In the basic two-pass algorithm:

• Compute $m \geq n$0.
• Compute $m \geq n$1.

For enhanced accuracy, $m \geq n$2 steps of a power iteration may be applied. With this optional refinement, initialize $m \geq n$3 and, for $m \geq n$4, alternate:

• $m \geq n$5
• $m \geq n$6

The parameter $m \geq n$7 determines the extent to which small singular values are suppressed, directly improving the subspace capture at the expense of additional passes through the data. This iterative enhancement is crucial for matrices with slowly-decaying singular value spectra, enabling the algorithm to adaptively optimize the low-rank subspace approximation (Kaloorazi et al., 2018).

3. Orthonormal Basis Construction

For both $m \geq n$8 and $m \geq n$9, thin QR decompositions are computed:

• $k \ll \min(m, n)$0 with $k \ll \min(m, n)$1,
• $k \ll \min(m, n)$2 with $k \ll \min(m, n)$3.

Here, $k \ll \min(m, n)$4 provides an approximate orthonormal basis for $k \ll \min(m, n)$5, while $k \ll \min(m, n)$6 approximates $k \ll \min(m, n)$7. This symmetric treatment of row and column subspaces—enabled by the two-sided sketch—provides empirical and theoretical advantages over one-sided randomized SVDs, as the basis quality directly impacts the accuracy of the final low-rank reconstruction (Kaloorazi et al., 2018).

4. Low-Rank Approximation and Error Bounds

The core matrix is constructed as $k \ll \min(m, n)$8, then subjected to a rank-$k \ll \min(m, n)$9 truncated singular value decomposition: $k$0, with $k$1. The low-rank approximation of $k$2 is lifted from the subspace: $k$3 This technique ensures that most computational effort is directed towards operations (matrix–matrix multiplies, QR decompositions) that scale favorably on parallel platforms. Error analysis establishes deterministic and average-case spectral/Frobenius-norm bounds. Specifically, for $k$4, the following holds for $k$5 or $k$6: $k$7 where projection errors diminish rapidly as $k$8 and $k$9 increase. The expected Frobenius norm error and lower bounds for approximate singular values are likewise specified, demonstrating near-optimality as the oversampling and power iteration parameters grow (Kaloorazi et al., 2018).

5. Algorithmic Workflow and Pseudocode

The basic two-pass SOR-SVD algorithm is summarized as:

Step	Description	Output/Notation
1	Draw $\hat{A}$0 (i.i.d. Gaussian)	$\hat{A}$1
2	Compute $\hat{A}$2	$\hat{A}$3
3	Compute $\hat{A}$4	$\hat{A}$5
4	Thin QR: $\hat{A}$6; $\hat{A}$7	$\hat{A}$8, $\hat{A}$9
5	Form $p \geq 2$0	$p \geq 2$1
6	Truncated SVD: $p \geq 2$2	$p \geq 2$3, $p \geq 2$4, $p \geq 2$5
7	Reconstruct: $p \geq 2$6	$p \geq 2$7

With $p \geq 2$8 power steps, steps 2–3 are replaced by the power iteration loop as detailed above (Kaloorazi et al., 2018).

6. Computational Complexity and Scalability

The algorithm is optimized for high performance on modern hardware. For $p \geq 2$9 and a matrix–vector multiply cost $\ell = k + p$0:

• Basic 3-pass: cost $\ell = k + p$1,
• 2-pass: cost $\ell = k + p$2,
• Additional $\ell = k + p$3 passes for power iterations.

Overall complexity is $\ell = k + p$4 flops, with dominant operations naturally suited to level-3 BLAS routines. Only $\ell = k + p$5–$\ell = k + p$6 data passes are required, with the option for even single-pass variants. All large-scale work—mat-mul and tall-skinny QR—are efficiently parallelizable. Communication-avoiding QR (CAQR) can further reduce data movement. The small SVD is performed in-core and is negligible in communication cost (Kaloorazi et al., 2018).

7. Relationship to Other Randomized SVD Algorithms

SOR-SVD generalizes and improves upon previous randomized SVD methods:

• The classical one-sided R-SVD (Halko–Martinsson–Tropp) computes a single subspace projection, forming only $\ell = k + p$7, whereas SOR-SVD builds both $\ell = k + p$8 and $\ell = k + p$9—enabling better control over both the column and row subspaces and yielding tighter error bounds.
• TSR-SVD (“two-sided randomized SVD”) sketches both subspaces with separate random matrices but is hampered by a poorly constructed $\Omega \in \mathbb{R}^{n \times \ell}$0; SOR-SVD resolves this by using the data-driven $\Omega \in \mathbb{R}^{n \times \ell}$1, thereby enhancing accuracy and stability.
• Compared to subspace iteration strategies that rely on repeated application of $\Omega \in \mathbb{R}^{n \times \ell}$2, SOR-SVD achieves comparable accuracy with fewer, less expensive orthogonalizations and improved communication efficiency (Kaloorazi et al., 2018).

SOR-SVD is thus positioned as a robust, practical, and highly parallelizable approach for low-rank approximation, supporting scalable deployment in large-scale applications such as robust PCA, where replacing dense SVD computations with SOR-SVD leads to dramatic performance improvements.