svd.cpp

Go to the documentation of this file.

/// @file svd/svd.cpp
/// @brief 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.
 
#include "linalg/svd/svd.hpp"
#include "linalg/factorization/qr.hpp"
#include <algorithm>
#include <cmath>
#include <stdexcept>
 
namespace num {
 
// Full SVD  -- one-sided Jacobi
 
SVDResult svd(const Matrix& A_in, Backend /*backend*/, real tol, idx max_sweeps) {
    constexpr real tiny = 1e-300;
    idx m = A_in.rows(), n = A_in.cols();
    idx r = std::min(m, n);   // rank of economy SVD
 
    // Work on a copy; rotate columns in-place
    Matrix A = A_in;
 
    // V accumulates right singular vectors (n x n), initialised to identity
    Matrix V(n, n, 0.0);
    for (idx i = 0; i < n; ++i) V(i,i) = 1.0;
 
    idx sweeps = 0;
    bool converged = false;
 
    for (idx sweep = 0; sweep < max_sweeps; ++sweep) {
        // Max relative off-orthogonality across column pairs
        real max_cos = 0;
        for (idx p = 0; p < r-1; ++p) {
            for (idx q = p+1; q < r; ++q) {
                // Inner products of columns p and q
                real alpha = 0, beta = 0, gamma = 0;
                for (idx i = 0; i < m; ++i) {
                    alpha += A(i,p) * A(i,p);
                    beta  += A(i,q) * A(i,q);
                    gamma += A(i,p) * A(i,q);
                }
                if (alpha < tiny || beta < tiny) continue;
 
                real cos_pq = std::abs(gamma) / std::sqrt(alpha * beta);
                max_cos = std::max(max_cos, cos_pq);
 
                if (cos_pq < tol) continue;  // already orthogonal enough
 
                // Jacobi rotation angle (identical structure to eig_sym)
                real zeta = (beta - alpha) / (2.0 * gamma);
                real t    = std::copysign(1.0, zeta)
                            / (std::abs(zeta) + std::sqrt(1.0 + zeta*zeta));
                real c    = 1.0 / std::sqrt(1.0 + t*t);
                real s    = c * t;
 
                // Update columns of A
                for (idx i = 0; i < m; ++i) {
                    real aip = A(i,p), aiq = A(i,q);
                    A(i,p) = c*aip - s*aiq;
                    A(i,q) = s*aip + c*aiq;
                }
 
                // Accumulate right singular vectors V
                for (idx i = 0; i < n; ++i) {
                    real vip = V(i,p), viq = V(i,q);
                    V(i,p) = c*vip - s*viq;
                    V(i,q) = s*vip + c*viq;
                }
            }
        }
 
        ++sweeps;
        if (max_cos < tol) { converged = true; break; }
    }
 
    // Singular values = column norms; left singular vectors = normalised cols
    Vector S(r);
    Matrix U(m, r, 0.0);
    for (idx j = 0; j < r; ++j) {
        real nrm = 0;
        for (idx i = 0; i < m; ++i) nrm += A(i,j)*A(i,j);
        S[j] = std::sqrt(nrm);
        if (S[j] > tiny)
            for (idx i = 0; i < m; ++i) U(i,j) = A(i,j) / S[j];
    }
 
    // Sort singular values descending; permute U, V columns to match
    for (idx i = 0; i < r-1; ++i) {
        idx max_j = i;
        for (idx j = i+1; j < r; ++j)
            if (S[j] > S[max_j]) max_j = j;
 
        if (max_j != i) {
            std::swap(S[i], S[max_j]);
            for (idx k = 0; k < m; ++k) std::swap(U(k,i), U(k,max_j));
            for (idx k = 0; k < n; ++k) std::swap(V(k,i), V(k,max_j));
        }
    }
 
    // Return economy Vt (r x n): transpose columns of V[:, 0..r-1]
    Matrix Vt(r, n, 0.0);
    for (idx i = 0; i < r; ++i)
        for (idx j = 0; j < n; ++j)
            Vt(i,j) = V(j,i);
 
    return {U, S, Vt, sweeps, converged};
}
 
// Randomized truncated SVD  -- Halko, Martinsson, Tropp (2011)
//
// Steps:
//   1. Draw Omega in R^{nx(k+p)}  (Gaussian random sketch)
//   2. Y = A * Omega
//   3. Q, _ = QR(Y)            (orthonormal basis for approx. range of A)
//   4. B = Qᵀ * A              (project to small k-dimensional subspace)
//   5. Ũ, Sigma, Vᵀ = svd(B)      (cheap: B is (k+p)xn, small)
//   6. U = Q * Ũ               (lift back to full dimension)
 
SVDResult svd_truncated(const Matrix& A, idx k,
                        Backend backend, idx oversampling, Rng* rng) {
    const idx m = A.rows(), n = A.cols();
    if (k == 0 || k > std::min(m, n))
        throw std::invalid_argument("svd_truncated: k out of range");
 
    const idx l = k + oversampling;   // sketch width
 
    // Optional local RNG if none supplied
    Rng local_rng;
    if (!rng) rng = &local_rng;
 
    // 1. Gaussian sketch matrix Omega in R^{nxl}
    Matrix Omega(n, l);
    for (idx j = 0; j < l; ++j)
        for (idx i = 0; i < n; ++i)
            Omega(i,j) = rng_normal(rng, 0.0, 1.0);
 
    // 2. Y = A * Omega  (m x l)
    Matrix Y(m, l, 0.0);
    matmul(A, Omega, Y, backend);
 
    // 3. Q = QR(Y).Q  (m x l, orthonormal columns)
    QRResult qr_res = qr(Y);
    const Matrix& Q = qr_res.Q;   // m x m; we only need first l columns
 
    // 4. B = Qᵀ * A  (l x n): B[i,j] = sum_k Q[k,i] * A[k,j]
    Matrix B(l, n, 0.0);
    for (idx i = 0; i < l; ++i)
        for (idx kk = 0; kk < m; ++kk) {
            const real q_ki = Q(kk, i);
            for (idx j = 0; j < n; ++j)
                B(i,j) += q_ki * A(kk,j);
        }
 
    // 5. SVD of B (small: l x n)
    SVDResult small = svd(B, backend);
 
    // 6. U = Q * Ũ[:, 0..k-1]  (m x k)
    Matrix U(m, k, 0.0);
    for (idx j = 0; j < k; ++j)
        for (idx i = 0; i < m; ++i)
            for (idx ii = 0; ii < l; ++ii)
                U(i,j) += Q(i,ii) * small.U(ii,j);
 
    // Slice top-k from small SVD
    Vector S(k);
    for (idx i = 0; i < k; ++i) S[i] = small.S[i];
 
    Matrix Vt(k, n, 0.0);
    for (idx i = 0; i < k; ++i)
        for (idx j = 0; j < n; ++j)
            Vt(i,j) = small.Vt(i,j);
 
    return {U, S, Vt, 0, true};
}
 
} // namespace num