1#ifndef IRLBA_COMPUTE_HPP
2#define IRLBA_COMPUTE_HPP
25template<
class Matrix_,
class EigenMatrix_,
class EigenVector_>
26void exact(
const Matrix_& matrix,
int requested_number, EigenMatrix_& outU, EigenMatrix_& outV, EigenVector_& outD) {
27 Eigen::BDCSVD<EigenMatrix_> svd(matrix.rows(), matrix.cols(), Eigen::ComputeThinU | Eigen::ComputeThinV);
29 if constexpr(internal::is_eigen<Matrix_>::value) {
79template<
class Matrix_,
class EigenMatrix_,
class EigenVector_>
87 throw std::runtime_error(
"requested number of singular values cannot be greater than the smaller matrix dimension");
90 throw std::runtime_error(
"requested number of singular values must be less than the smaller matrix dimension for IRLBA iterations");
97 return std::make_pair(
true, 0);
106 if (
init->size() != V.rows()) {
107 throw std::runtime_error(
"initialization vector does not have expected number of rows");
111 internal::fill_with_random_normals(V, 0,
eng);
113 V.col(0) /= V.col(0).norm();
115 bool converged =
false;
118 Eigen::JacobiSVD<EigenMatrix_>
svd(
work,
work, Eigen::ComputeThinU | Eigen::ComputeThinV);
120 internal::LanczosWorkspace<EigenVector_, Matrix_>
lptmp(
matrix);
132 typename EigenMatrix_::Scalar
svtol =
options.singular_value_ratio_tolerance;
133 typename EigenMatrix_::Scalar
tol =
options.convergence_tolerance;
149 const auto&
BS =
svd.singularValues();
150 const auto&
BU =
svd.matrixU();
151 const auto&
BV =
svd.matrixV();
171 auto Smax = *std::max_element(
BS.begin(),
BS.end());
174 for (Eigen::Index
j = 0;
j <
work; ++
j) {
207 Vtmp.leftCols(
k).noalias() = V *
BV.leftCols(
k);
208 V.leftCols(
k) =
Vtmp.leftCols(
k);
218 Wtmp.leftCols(
k).noalias() =
W *
BU.leftCols(
k);
219 W.leftCols(
k) =
Wtmp.leftCols(
k);
222 for (Eigen::Index
l = 0;
l <
k; ++
l) {
246 return std::make_pair(converged,
iter + 1);
254template<
class EigenMatrix_,
class EigenVector_>
300template<
class EigenMatrix_ = Eigen::MatrixXd,
class EigenVector_ = Eigen::VectorXd,
class Matrix_>
338template<
class Matrix_,
class EigenMatrix_,
class EigenVector_>
350 throw std::runtime_error(
"cannot center with no observations");
358 throw std::runtime_error(
"cannot scale with fewer than two observations");
363 for (Eigen::Index
i = 0;
i <
nc; ++
i) {
371 typename EigenMatrix_::Scalar
var = 0;
413template<
class EigenMatrix_ = Eigen::MatrixXd,
class EigenVector_ = Eigen::VectorXd,
class Matrix_>
Wrapper for a centered matrix.
Definition wrappers.hpp:205
Implements IRLBA for approximate SVD.
Definition compute.hpp:18
std::pair< bool, int > compute(const Matrix_ &matrix, Eigen::Index number, EigenMatrix_ &outU, EigenMatrix_ &outV, EigenVector_ &outD, const Options &options)
Definition compute.hpp:80
EigenMatrix_ wrapped_realize(const Matrix_ &matrix)
Definition wrappers.hpp:185
Options for the IRLBA algorithm.
Definition Options.hpp:16
Result of the IRLBA-based decomposition.
Definition compute.hpp:255
EigenVector_ D
Definition compute.hpp:274
int iterations
Definition compute.hpp:279
EigenMatrix_ U
Definition compute.hpp:261
EigenMatrix_ V
Definition compute.hpp:268
bool converged
Definition compute.hpp:284