---
title: Subspace-Orbit Randomized SVD
url: https://www.emergentmind.com/topics/subspace-orbit-randomized-svd-sor-svd
type: topic
---

# Subspace-Orbit Randomized SVD

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 [1804.00462].

## 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 [1804.00462].

## 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 $T_1 = A \Omega$.
- Compute $T_2 = A^\top T_1$.

For enhanced accuracy, $q$ steps of a power iteration may be applied. With this optional refinement, initialize $T_2 \leftarrow \Omega$ and, for $i=0,\ldots,q$, alternate:
- $T_1 \leftarrow A T_2$
- $T_2 \leftarrow A^\top T_1$

The parameter $q$ 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 [1804.00462].

## 3. Orthonormal Basis Construction

For both $T_1$ and $T_2$, thin QR decompositions are computed:
- $T_1 = Q_1 R_1$ with $Q_1 \in \mathbb{R}^{m \times \ell}$,
- $T_2 = Q_2 R_2$ with $Q_2 \in \mathbb{R}^{n \times \ell}$.

Here, $Q_1$ provides an approximate orthonormal basis for $\mathrm{Range}(A)$, while $Q_2$ approximates $\mathrm{Range}(A^\top)$. 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 [1804.00462].

## 4. Low-Rank Approximation and Error Bounds

The core matrix is constructed as $M = Q_1^\top A Q_2 \in \mathbb{R}^{\ell \times \ell}$, then subjected to a rank-$k$ truncated singular value decomposition: $M \approx \tilde{U}_k \Sigma_k \tilde{V}_k^\top$, with $\Sigma_k = \mathrm{diag}(\hat{\sigma}_1,\ldots,\hat{\sigma}_k)$. The low-rank approximation of $A$ is lifted from the subspace:
\[
\hat{A} = (Q_1 \tilde{U}_k) \Sigma_k (Q_2 \tilde{V}_k)^\top.
\]
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 $A = U_k \Sigma_k V_k^\top + U_0 \Sigma_0 V_0^\top = A_k + A_0$, the following holds for $p=2$ or $F$:
\[
\|A - \hat{A}\|_p \leq \|A_0\|_p + \|A_k - Q_1 Q_1^\top A_k\|_F + \|A_k - A_k Q_2 Q_2^\top\|_F,
\]
where projection errors diminish rapidly as $p$ and $q$ 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 [1804.00462].

## 5. Algorithmic Workflow and Pseudocode

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

| Step | Description                                | Output/Notation          |
|------|--------------------------------------------|-------------------------|
| 1    | Draw $\Omega \in \mathbb{R}^{n \times \ell}$ (i.i.d. Gaussian) | $\Omega$                |
| 2    | Compute $T_1 = A \Omega$                   | $T_1$                   |
| 3    | Compute $T_2 = A^\top T_1$                 | $T_2$                   |
| 4    | Thin QR: $[Q_1, R_1] = \text{qr}(T_1)$; $[Q_2, R_2] = \text{qr}(T_2)$ | $Q_1$, $Q_2$            |
| 5    | Form $M = Q_1^\top A Q_2$                  | $M \in \mathbb{R}^{\ell \times \ell}$ |
| 6    | Truncated SVD: $[\tilde{U}_k, \Sigma_k, \tilde{V}_k] = \text{svds}(M, k)$ | $\tilde{U}_k$, $\Sigma_k$, $\tilde{V}_k$ |
| 7    | Reconstruct: $\hat{A} = (Q_1 \tilde{U}_k) \Sigma_k (Q_2 \tilde{V}_k)^\top$ | $\hat{A}$               |

With $q$ power steps, steps 2–3 are replaced by the power iteration loop as detailed above [1804.00462].

## 6. Computational Complexity and Scalability

The algorithm is optimized for high performance on modern hardware. For $\ell = k + p \ll n$ and a matrix–vector multiply cost $C_{mult} \approx 2mn$:
- Basic 3-pass: cost $\approx 3\ell C_{mult} + O(\ell^2(m+n) + \ell^3)$,
- 2-pass: cost $\approx 2\ell C_{mult} + O(\ell^2(m+n))$,
- Additional $2q$ passes for power iterations.

Overall complexity is $O(mnk)$ flops, with dominant operations naturally suited to level-3 BLAS routines. Only $2$–$3$ 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 [1804.00462].

## 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 $Q_1$, whereas SOR-SVD builds both $Q_1$ and $Q_2$—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 $Q_2$; SOR-SVD resolves this by using the data-driven $T_1 = A\Omega$, thereby enhancing accuracy and stability.
- Compared to subspace iteration strategies that rely on repeated application of $(A A^\top)^q A\Omega$, SOR-SVD achieves comparable accuracy with fewer, less expensive orthogonalizations and improved communication efficiency [1804.00462].

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.

Source: https://www.emergentmind.com/topics/subspace-orbit-randomized-svd-sor-svd