compute.hpp

Go to the documentation of this file.

#ifndef IRLBA_COMPUTE_HPP
#define IRLBA_COMPUTE_HPP
 
#include "Eigen/Dense"
#include "utils.hpp"
#include "lanczos.hpp"
 
#include <cmath>
#include <cstdint>
#include <random>
#include <stdexcept>
 
namespace irlba {
 
namespace internal {
 
template<class Matrix_, class EigenMatrix_, class EigenVector_>
void exact(const Matrix_& matrix, int requested_number, EigenMatrix_& outU, EigenMatrix_& outV, EigenVector_& outD) {
    Eigen::JacobiSVD<EigenMatrix_> svd(matrix.rows(), matrix.cols(), Eigen::ComputeThinU | Eigen::ComputeThinV);
 
    if constexpr(internal::is_eigen<Matrix_>::value) {
        svd.compute(matrix);
    } else {
        auto adjusted = wrapped_realize<EigenMatrix_>(matrix);
        svd.compute(adjusted);
    }
 
    outD.resize(requested_number);
    outD = svd.singularValues().head(requested_number);
 
    outU.resize(matrix.rows(), requested_number);
    outU = svd.matrixU().leftCols(requested_number);
 
    outV.resize(matrix.cols(), requested_number);
    outV = svd.matrixV().leftCols(requested_number);
 
    return;
}
 
}
template<class Matrix_, class EigenMatrix_, class EigenVector_>
std::pair<bool, int> compute(const Matrix_& matrix, Eigen::Index number, EigenMatrix_& outU, EigenMatrix_& outV, EigenVector_& outD, const Options& options) {
    Eigen::Index smaller = std::min(matrix.rows(), matrix.cols());
    Eigen::Index requested_number = number;
    if (requested_number > smaller) {
        if (options.cap_number) {
            requested_number = smaller;
        } else {
            throw std::runtime_error("requested number of singular values cannot be greater than the smaller matrix dimension");
        }
    } else if (requested_number == smaller && !options.exact_for_large_number) {
        throw std::runtime_error("requested number of singular values must be less than the smaller matrix dimension for IRLBA iterations");
    }
 
    // Falling back to an exact SVD for small matrices or if the requested number is too large 
    // (not enough of a workspace). Hey, I don't make the rules.
    if ((options.exact_for_small_matrix && smaller < 6) || (options.exact_for_large_number && requested_number * 2 >= smaller)) {
        internal::exact(matrix, requested_number, outU, outV, outD);
        return std::make_pair(true, 0);
    }
 
    Eigen::Index work = std::min(requested_number + options.extra_work, smaller);
 
    EigenMatrix_ V(matrix.cols(), work);
    std::mt19937_64 eng(options.seed);
    if (options.initial) {
        auto init = reinterpret_cast<EigenVector_*>(options.initial);
        if (init->size() != V.rows()) {
            throw std::runtime_error("initialization vector does not have expected number of rows");
        }
        V.col(0) = *init;
    } else {
        internal::fill_with_random_normals(V, 0, eng);
    }
    V.col(0) /= V.col(0).norm();
 
    bool converged = false;
    int iter = 0;
    Eigen::Index k = 0;
    Eigen::JacobiSVD<EigenMatrix_> svd(work, work, Eigen::ComputeThinU | Eigen::ComputeThinV);
 
    internal::LanczosWorkspace<EigenVector_, Matrix_> lptmp(matrix);
 
    EigenMatrix_ W(matrix.rows(), work);
    EigenMatrix_ Wtmp(matrix.rows(), work);
    EigenMatrix_ Vtmp(matrix.cols(), work);
 
    EigenMatrix_ B(work, work);
    B.setZero(work, work);
    EigenVector_ res(work);
    EigenVector_ F(matrix.cols());
 
    EigenVector_ prevS(work);
    typename EigenMatrix_::Scalar svtol = options.singular_value_ratio_tolerance;
    typename EigenMatrix_::Scalar tol = options.convergence_tolerance;
    typename EigenMatrix_::Scalar svtol_actual = (svtol >= 0 ? svtol : tol);
 
    for (; iter < options.max_iterations; ++iter) {
        // Technically, this is only a 'true' Lanczos bidiagonalization
        // when k = 0. All other times, we're just recycling the machinery,
        // see the text below Equation 3.11 in Baglama and Reichel.
        internal::run_lanczos_bidiagonalization(matrix, W, V, B, eng, lptmp, k, options);
 
//            if (iter < 2) {
//                std::cout << "B is currently:\n" << B << std::endl;
//                std::cout << "W is currently:\n" << W << std::endl;
//                std::cout << "V is currently:\n" << V << std::endl;
//            }
 
        svd.compute(B);
        const auto& BS = svd.singularValues();
        const auto& BU = svd.matrixU();
        const auto& BV = svd.matrixV();
 
        // Checking for convergence.
        if (B(work - 1, work - 1) == 0) { // a.k.a. the final value of 'S' from the Lanczos iterations.
            converged = true;
            break;
        }
 
        const auto& F = lptmp.F; // i.e., the residuals, see Algorithm 2.1 of Baglama and Reichel.
        auto R_F = F.norm();
 
        // Computes the convergence criterion defined in on the LHS of
        // Equation 2.13 of Baglama and Riechel. (e_m is literally the unit
        // vector along the m-th dimension, so it ends up just being the
        // m-th row of the U matrix.) We expose it here as we will be
        // re-using the same values to update B, see below.
        res = R_F * BU.row(work - 1);
 
        Eigen::Index n_converged = 0;
        if (iter > 0) {
            auto Smax = *std::max_element(BS.begin(), BS.end());
            auto threshold = Smax * tol;
 
            for (Eigen::Index j = 0; j < work; ++j) {
                if (std::abs(res[j]) <= threshold) {
                    auto ratio = std::abs(BS[j] - prevS[j]) / BS[j];
                    if (ratio <= svtol_actual) {
                        ++n_converged;
                    }
                }
            }
 
            if (n_converged >= requested_number) {
                converged = true;
                break;
            }
        }
        prevS = BS;
 
        // Setting 'k'. This looks kinda weird, but this is deliberate,
        // see the text below Algorithm 3.1 of Baglama and Reichel.
        //
        // - Our n_converged is their k'.
        // - Our requested_number is their k.
        // - Our work is their m.
        // - Our k is their k + k''.
        if (n_converged + requested_number > k) {
            k = n_converged + requested_number;
        }
        if (k > work - 3) {
            k = work - 3;
        }
        if (k < 1) {
            k = 1;
        }
 
        Vtmp.leftCols(k).noalias() = V * BV.leftCols(k);
        V.leftCols(k) = Vtmp.leftCols(k);
 
        // See Equation 3.2 of Baglama and Reichel, where our 'V' is their
        // 'P', and our 'F / R_F' is their 'p_{m+1}' (Equation 2.2).  'F'
        // was orthogonal to the old 'V' and so it should still be
        // orthogonal to the new left-most columns of 'V'; the input
        // expectations of the Lanczos bidiagonalization are still met, so
        // V is ok to use in the next run_lanczos_bidiagonalization().
        V.col(k) = F / R_F; 
 
        Wtmp.leftCols(k).noalias() = W * BU.leftCols(k);
        W.leftCols(k) = Wtmp.leftCols(k);
 
        B.setZero(work, work);
        for (Eigen::Index l = 0; l < k; ++l) {
            B(l, l) = BS[l];
 
            // This assignment looks weird but is deliberate, see Equation
            // 3.6 of Baglama and Reichel. Happily, this is the same value
            // used to determine convergence in 2.13, so we just re-use it.
            // (See the equation just above Equation 3.5; I think they 
            // misplaced a tilde on the final 'u', given no other 'u' has
            // a B_m superscript as well as a tilde.)
            B(l, k) = res[l]; 
        }
    }
 
    // See Equation 2.11 of Baglama and Reichel for how to get from B's
    // singular triplets to mat's singular triplets.
    outD.resize(requested_number);
    outD = svd.singularValues().head(requested_number);
 
    outU.resize(matrix.rows(), requested_number);
    outU.noalias() = W * svd.matrixU().leftCols(requested_number);
 
    outV.resize(matrix.cols(), requested_number);
    outV.noalias() = V * svd.matrixV().leftCols(requested_number);
 
    return std::make_pair(converged, iter + 1);
}
 
template<class EigenMatrix_, class EigenVector_>
struct Results {
    EigenMatrix_ U;
 
    EigenMatrix_ V;
 
    EigenVector_ D;
 
    int iterations;
 
    bool converged;
};
 
template<class EigenMatrix_ = Eigen::MatrixXd, class EigenVector_ = Eigen::VectorXd, class Matrix_>
Results<EigenMatrix_, EigenVector_> compute(const Matrix_& matrix, Eigen::Index number, const Options& options) {
    Results<EigenMatrix_, EigenVector_> output;
    auto stats = compute(matrix, number, output.U, output.V, output.D, options);
    output.converged = stats.first;
    output.iterations = stats.second;
    return output;
}
 
template<class Matrix_, class EigenMatrix_, class EigenVector_>
std::pair<bool, int> compute(const Matrix_& matrix, bool center, bool scale, Eigen::Index number, EigenMatrix_& outU, EigenMatrix_& outV, EigenVector_& outD, const Options& options) {
    if (!scale && !center) {
        return compute(matrix, number, outU, outV, outD, options);
    }
 
    auto nr = matrix.rows();
    auto nc = matrix.cols();
 
    Eigen::VectorXd center0;
    if (center) {
        if (nr < 1) {
            throw std::runtime_error("cannot center with no observations");    
        }
        center0.resize(nc);
    }
 
    Eigen::VectorXd scale0;
    if (scale) {
        if (nr < 2) {
            throw std::runtime_error("cannot scale with fewer than two observations");    
        }
        scale0.resize(nc);
    }
 
    for (Eigen::Index i = 0; i < nc; ++i) {
        double mean = 0;
        if (center) {
            mean = matrix.col(i).sum() / nr;
            center0[i] = mean;
        }
        if (scale) {
            EigenVector_ current = matrix.col(i); // force it to be a Vector, even if it's a sparse matrix.
            typename EigenMatrix_::Scalar var = 0;
            for (auto x : current) {
                var += (x - mean)*(x - mean);
            }
 
            if (var) {
                scale0[i] = std::sqrt(var/(nr - 1));
            } else {
                scale0[i] = 1;
            }
        }
    }
 
    if (center) {
        Centered<Matrix_, EigenVector_> centered(matrix, center0);
        if (scale) {
            auto centered_scaled = make_Scaled<true>(centered, scale0, true);
            return compute(centered_scaled, number, outU, outV, outD, options);
        } else {
            return compute(centered, number, outU, outV, outD, options);
        }
    } else {
        auto scaled = make_Scaled<true>(matrix, scale0, true);
        return compute(scaled, number, outU, outV, outD, options);
    }
}
 
template<class EigenMatrix_ = Eigen::MatrixXd, class EigenVector_ = Eigen::VectorXd, class Matrix_>
Results<EigenMatrix_, EigenVector_> compute(const Matrix_& matrix, bool center, bool scale, Eigen::Index number, const Options& options) {
    Results<EigenMatrix_, EigenVector_> output;
    auto stats = compute(matrix, center, scale, number, output.U, output.V, output.D, options);
    output.converged = stats.first;
    output.iterations = stats.second;
    return output;
}
 
}
 
#endif