svd.cpp

Go to the documentation of this file.

/// @file svd/svd.cpp
/// @brief SVD dispatcher + randomized truncated SVD.
 
#include "linalg/svd/svd.hpp"
#include "backends/lapack/impl.hpp"
#include "backends/seq/impl.hpp"
#include "linalg/factorization/qr.hpp"
 
namespace num {
 
SVDResult svd(const Matrix& A_in, Backend backend, real tol, idx max_sweeps) {
    switch (backend) {
        case Backend::lapack:
            return backends::lapack::svd(A_in);
        default:
            return backends::seq::svd(A_in, tol, max_sweeps);
    }
}
 
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;
 
    Rng local_rng;
    if (!rng)
        rng = &local_rng;
 
    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);
 
    Matrix Y(m, l, 0.0);
    matmul(A, Omega, Y, backend);
 
    QRResult qr_res = qr(Y);
    const Matrix& Q = qr_res.Q;
 
    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);
        }
 
    SVDResult small = svd(B, backend);
 
    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);
 
    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