svd.cpp File Reference

SVD via one-sided Jacobi + randomized truncated SVD. More...

#include "linalg/svd/svd.hpp"
#include "linalg/factorization/qr.hpp"
#include <algorithm>
#include <cmath>
#include <stdexcept>

Go to the source code of this file.

Namespaces
namespace	num

Functions
SVDResult	num::svd (const Matrix &A, Backend backend=Backend::seq, real tol=1e-12, idx max_sweeps=100)
	Full SVD of an mxn matrix via one-sided Jacobi.
SVDResult	num::svd_truncated (const Matrix &A, idx k, Backend backend=default_backend, idx oversampling=10, Rng *rng=nullptr)
	Randomized truncated SVD – top-k singular triplets.

Detailed Description

SVD via one-sided Jacobi + randomized truncated SVD.

– One-sided Jacobi SVD -----------------------------------------------—

Apply Givens rotations G to the columns of A from the right until the columns are mutually orthogonal. Each rotation zeros (AᵀA)[p,q]:

Given col_p and col_q of the working matrix A: alpha = col_p * col_p, beta = col_q * col_q, gamma = col_p * col_q zeta = (beta - alpha) / (2gamma) t = sign(zeta) / (|zeta| + sqrt(1 + zeta^2)) <- same formula as Jacobi eig c = 1/sqrt(1+t^2), s = c*t

col_p' = c*col_p + s*col_q col_q' = -s*col_p + c*col_q V[:,p]' = c*V[:,p] + s*V[:,q] <- accumulate right sing. vecs V[:,q]' = -s*V[:,p] + c*V[:,q]

After convergence: S[i] = ||A[:,i]||, U[:,i] = A[:,i]/S[i]

– Randomized SVD -----------------------------------------------------—

Halko, Martinsson, Tropp (2011) "Finding Structure with Randomness". The key insight: if Y = A*Omega has range ~= range(A), a small SVD of Qᵀ*A gives the dominant singular triplets cheaply.

Definition in file svd.cpp.