irlba
A C++ library for IRLBA
Loading...
Searching...
No Matches
compute.hpp
Go to the documentation of this file.
1#ifndef IRLBA_COMPUTE_HPP
2#define IRLBA_COMPUTE_HPP
3
4#include "Eigen/Dense"
5#include "utils.hpp"
6#include "lanczos.hpp"
7
8#include <cmath>
9#include <cstdint>
10#include <random>
11#include <stdexcept>
12
18namespace irlba {
19
23namespace internal {
24
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);
28
29 if constexpr(internal::is_eigen<Matrix_>::value) {
30 svd.compute(matrix);
31 } else {
33 svd.compute(adjusted);
34 }
35
36 outD.resize(requested_number);
37 outD = svd.singularValues().head(requested_number);
38
39 outU.resize(matrix.rows(), requested_number);
40 outU = svd.matrixU().leftCols(requested_number);
41
42 outV.resize(matrix.cols(), requested_number);
43 outV = svd.matrixV().leftCols(requested_number);
44
45 return;
46}
47
48}
79template<class Matrix_, class EigenMatrix_, class EigenVector_>
80std::pair<bool, int> compute(const Matrix_& matrix, Eigen::Index number, EigenMatrix_& outU, EigenMatrix_& outV, EigenVector_& outD, const Options& options) {
81 Eigen::Index smaller = std::min(matrix.rows(), matrix.cols());
82 Eigen::Index requested_number = number;
84 if (options.cap_number) {
86 } else {
87 throw std::runtime_error("requested number of singular values cannot be greater than the smaller matrix dimension");
88 }
89 } else if (requested_number == smaller && !options.exact_for_large_number) {
90 throw std::runtime_error("requested number of singular values must be less than the smaller matrix dimension for IRLBA iterations");
91 }
92
93 // Falling back to an exact SVD for small matrices or if the requested number is too large
94 // (not enough of a workspace). Hey, I don't make the rules.
95 if ((options.exact_for_small_matrix && smaller < 6) || (options.exact_for_large_number && requested_number * 2 >= smaller)) {
96 internal::exact(matrix, requested_number, outU, outV, outD);
97 return std::make_pair(true, 0);
98 }
99
100 Eigen::Index work = std::min(requested_number + options.extra_work, smaller);
101
102 EigenMatrix_ V(matrix.cols(), work);
103 std::mt19937_64 eng(options.seed);
104 if (options.initial) {
105 auto init = reinterpret_cast<EigenVector_*>(options.initial);
106 if (init->size() != V.rows()) {
107 throw std::runtime_error("initialization vector does not have expected number of rows");
108 }
109 V.col(0) = *init;
110 } else {
111 internal::fill_with_random_normals(V, 0, eng);
112 }
113 V.col(0) /= V.col(0).norm();
114
115 bool converged = false;
116 int iter = 0;
117 Eigen::Index k = 0;
118 Eigen::JacobiSVD<EigenMatrix_> svd(work, work, Eigen::ComputeThinU | Eigen::ComputeThinV);
119
120 internal::LanczosWorkspace<EigenVector_, Matrix_> lptmp(matrix);
121
122 EigenMatrix_ W(matrix.rows(), work);
123 EigenMatrix_ Wtmp(matrix.rows(), work);
124 EigenMatrix_ Vtmp(matrix.cols(), work);
125
127 B.setZero(work, work);
129 EigenVector_ F(matrix.cols());
130
132 typename EigenMatrix_::Scalar svtol = options.singular_value_ratio_tolerance;
133 typename EigenMatrix_::Scalar tol = options.convergence_tolerance;
134 typename EigenMatrix_::Scalar svtol_actual = (svtol >= 0 ? svtol : tol);
135
136 for (; iter < options.max_iterations; ++iter) {
137 // Technically, this is only a 'true' Lanczos bidiagonalization
138 // when k = 0. All other times, we're just recycling the machinery,
139 // see the text below Equation 3.11 in Baglama and Reichel.
140 internal::run_lanczos_bidiagonalization(matrix, W, V, B, eng, lptmp, k, options);
141
142// if (iter < 2) {
143// std::cout << "B is currently:\n" << B << std::endl;
144// std::cout << "W is currently:\n" << W << std::endl;
145// std::cout << "V is currently:\n" << V << std::endl;
146// }
147
148 svd.compute(B);
149 const auto& BS = svd.singularValues();
150 const auto& BU = svd.matrixU();
151 const auto& BV = svd.matrixV();
152
153 // Checking for convergence.
154 if (B(work - 1, work - 1) == 0) { // a.k.a. the final value of 'S' from the Lanczos iterations.
155 converged = true;
156 break;
157 }
158
159 const auto& F = lptmp.F; // i.e., the residuals, see Algorithm 2.1 of Baglama and Reichel.
160 auto R_F = F.norm();
161
162 // Computes the convergence criterion defined in on the LHS of
163 // Equation 2.13 of Baglama and Riechel. (e_m is literally the unit
164 // vector along the m-th dimension, so it ends up just being the
165 // m-th row of the U matrix.) We expose it here as we will be
166 // re-using the same values to update B, see below.
167 res = R_F * BU.row(work - 1);
168
169 Eigen::Index n_converged = 0;
170 if (iter > 0) {
171 auto Smax = *std::max_element(BS.begin(), BS.end());
172 auto threshold = Smax * tol;
173
174 for (Eigen::Index j = 0; j < work; ++j) {
175 if (std::abs(res[j]) <= threshold) {
176 auto ratio = std::abs(BS[j] - prevS[j]) / BS[j];
177 if (ratio <= svtol_actual) {
178 ++n_converged;
179 }
180 }
181 }
182
184 converged = true;
185 break;
186 }
187 }
188 prevS = BS;
189
190 // Setting 'k'. This looks kinda weird, but this is deliberate,
191 // see the text below Algorithm 3.1 of Baglama and Reichel.
192 //
193 // - Our n_converged is their k'.
194 // - Our requested_number is their k.
195 // - Our work is their m.
196 // - Our k is their k + k''.
199 }
200 if (k > work - 3) {
201 k = work - 3;
202 }
203 if (k < 1) {
204 k = 1;
205 }
206
207 Vtmp.leftCols(k).noalias() = V * BV.leftCols(k);
208 V.leftCols(k) = Vtmp.leftCols(k);
209
210 // See Equation 3.2 of Baglama and Reichel, where our 'V' is their
211 // 'P', and our 'F / R_F' is their 'p_{m+1}' (Equation 2.2). 'F'
212 // was orthogonal to the old 'V' and so it should still be
213 // orthogonal to the new left-most columns of 'V'; the input
214 // expectations of the Lanczos bidiagonalization are still met, so
215 // V is ok to use in the next run_lanczos_bidiagonalization().
216 V.col(k) = F / R_F;
217
218 Wtmp.leftCols(k).noalias() = W * BU.leftCols(k);
219 W.leftCols(k) = Wtmp.leftCols(k);
220
221 B.setZero(work, work);
222 for (Eigen::Index l = 0; l < k; ++l) {
223 B(l, l) = BS[l];
224
225 // This assignment looks weird but is deliberate, see Equation
226 // 3.6 of Baglama and Reichel. Happily, this is the same value
227 // used to determine convergence in 2.13, so we just re-use it.
228 // (See the equation just above Equation 3.5; I think they
229 // misplaced a tilde on the final 'u', given no other 'u' has
230 // a B_m superscript as well as a tilde.)
231 B(l, k) = res[l];
232 }
233 }
234
235 // See Equation 2.11 of Baglama and Reichel for how to get from B's
236 // singular triplets to mat's singular triplets.
237 outD.resize(requested_number);
238 outD = svd.singularValues().head(requested_number);
239
240 outU.resize(matrix.rows(), requested_number);
241 outU.noalias() = W * svd.matrixU().leftCols(requested_number);
242
243 outV.resize(matrix.cols(), requested_number);
244 outV.noalias() = V * svd.matrixV().leftCols(requested_number);
245
246 return std::make_pair(converged, iter + 1);
247}
248
254template<class EigenMatrix_, class EigenVector_>
286
300template<class EigenMatrix_ = Eigen::MatrixXd, class EigenVector_ = Eigen::VectorXd, class Matrix_>
304 output.converged = stats.first;
305 output.iterations = stats.second;
306 return output;
307}
308
338template<class Matrix_, class EigenMatrix_, class EigenVector_>
339std::pair<bool, int> compute(const Matrix_& matrix, bool center, bool scale, Eigen::Index number, EigenMatrix_& outU, EigenMatrix_& outV, EigenVector_& outD, const Options& options) {
340 if (!scale && !center) {
342 }
343
344 auto nr = matrix.rows();
345 auto nc = matrix.cols();
346
347 Eigen::VectorXd center0;
348 if (center) {
349 if (nr < 1) {
350 throw std::runtime_error("cannot center with no observations");
351 }
352 center0.resize(nc);
353 }
354
355 Eigen::VectorXd scale0;
356 if (scale) {
357 if (nr < 2) {
358 throw std::runtime_error("cannot scale with fewer than two observations");
359 }
360 scale0.resize(nc);
361 }
362
363 for (Eigen::Index i = 0; i < nc; ++i) {
364 double mean = 0;
365 if (center) {
366 mean = matrix.col(i).sum() / nr;
367 center0[i] = mean;
368 }
369 if (scale) {
370 EigenVector_ current = matrix.col(i); // force it to be a Vector, even if it's a sparse matrix.
371 typename EigenMatrix_::Scalar var = 0;
372 for (auto x : current) {
373 var += (x - mean)*(x - mean);
374 }
375
376 if (var) {
377 scale0[i] = std::sqrt(var/(nr - 1));
378 } else {
379 scale0[i] = 1;
380 }
381 }
382 }
383
384 if (center) {
386 if (scale) {
389 } else {
391 }
392 } else {
393 auto scaled = make_Scaled<true>(matrix, scale0, true);
395 }
396}
397
413template<class EigenMatrix_ = Eigen::MatrixXd, class EigenVector_ = Eigen::VectorXd, class Matrix_>
417 output.converged = stats.first;
418 output.iterations = stats.second;
419 return output;
420}
421
422}
423
424#endif
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