From 4837520e24553b77ede2910b52a2f7fb3415ec01 Mon Sep 17 00:00:00 2001 From: Tobias Wood Date: Tue, 24 May 2022 22:41:23 +0100 Subject: [PATCH] Implements support for Complex --- include/irlba/irlba.hpp | 118 +++++++++++++++++++------------------ include/irlba/lanczos.hpp | 24 ++++---- include/irlba/utils.hpp | 28 +++++---- include/irlba/wrappers.hpp | 52 ++++++++-------- tests/src/NormalSampler.h | 11 ++++ tests/src/compare.h | 13 ++-- tests/src/irlba.cpp | 16 +++++ tests/src/utils.cpp | 12 +++- tests/src/wrappers.cpp | 8 +-- 9 files changed, 165 insertions(+), 117 deletions(-) diff --git a/include/irlba/irlba.hpp b/include/irlba/irlba.hpp index f683a58..b3b28ef 100644 --- a/include/irlba/irlba.hpp +++ b/include/irlba/irlba.hpp @@ -25,8 +25,12 @@ namespace irlba { * This is heavily derived from the C code in the [**irlba** package](https://github.com/bwlewis/irlba), * with refactoring into C++ to use Eigen instead of LAPACK for much of the matrix algebra. */ +template class Irlba { public: + using Scalar = typename MatrixType::Scalar; + using VectorType = Eigen::Vector; + using RealVectorType = Eigen::Vector::Real, MatrixType::RowsAtCompileTime>; /** * @brief Default parameter values. */ @@ -52,14 +56,14 @@ class Irlba { static constexpr uint64_t seed = std::mt19937_64::default_seed; }; private: - LanczosBidiagonalization lp; + LanczosBidiagonalization lp; int number = Defaults::number; int extra_work = Defaults::extra_work; int maxit = Defaults::maxit; uint64_t seed = Defaults::seed; - ConvergenceTest convtest; + ConvergenceTest convtest; public: /** @@ -121,7 +125,7 @@ class Irlba { * * @return A reference to the `Irlba` instance. */ - Irlba& set_invariant_tolerance(double e = LanczosBidiagonalization::Defaults::epsilon) { + Irlba& set_invariant_tolerance(double e = LanczosBidiagonalization::Defaults::epsilon) { lp.set_epsilon(e); return *this; } @@ -133,7 +137,7 @@ class Irlba { * * @return A reference to the `Irlba` instance. */ - Irlba& set_convergence_tolerance(double t = ConvergenceTest::Defaults::tol) { + Irlba& set_convergence_tolerance(double t = ConvergenceTest::Defaults::tol) { convtest.set_tol(t); return *this; } @@ -145,7 +149,7 @@ class Irlba { * * @return A reference to the `Irlba` instance. */ - Irlba& set_singular_value_ratio_tolerance(double t = ConvergenceTest::Defaults::svtol) { + Irlba& set_singular_value_ratio_tolerance(double t = ConvergenceTest::Defaults::svtol) { convtest.set_svtol(t); return *this; } @@ -187,14 +191,14 @@ class Irlba { const M& mat, bool center, bool scale, - Eigen::MatrixXd& outU, - Eigen::MatrixXd& outV, - Eigen::VectorXd& outD, + MatrixType& outU, + MatrixType& outV, + RealVectorType& outD, Engine* eng = null_rng(), - Eigen::VectorXd* init = NULL) + VectorType* init = NULL) { if (scale || center) { - Eigen::VectorXd center0, scale0; + VectorType center0, scale0; if (center) { if (mat.rows() < 1) { @@ -211,14 +215,14 @@ class Irlba { } for (Eigen::Index i = 0; i < mat.cols(); ++i) { - double mean = 0; + Scalar mean = 0; if (center) { mean = mat.col(i).sum() / mat.rows(); center0[i] = mean; } if (scale) { - Eigen::VectorXd current = mat.col(i); // force it to be a VectorXd, even if it's a sparse matrix. - double var = 0; + VectorType current = mat.col(i); // force it to be a VectorXd, even if it's a sparse matrix. + Scalar var = 0; for (auto x : current) { var += (x - mean)*(x - mean); } @@ -232,7 +236,7 @@ class Irlba { } if (center) { - Centered centered(&mat, ¢er0); + Centered centered(&mat, ¢er0); if (scale) { Scaled centered_scaled(¢ered, &scale0); return run(centered_scaled, outU, outV, outD, eng, init); @@ -240,7 +244,7 @@ class Irlba { return run(centered, outU, outV, outD, eng, init); } } else { - Scaled scaled(&mat, &scale0); + Scaled scaled(&mat, &scale0); return run(scaled, outU, outV, outD, eng, init); } } else { @@ -277,31 +281,31 @@ class Irlba { * - A `rows()` method that returns the number of rows. * - A `cols()` method that returns the number of columns. * - One of the following for matrix-vector multiplication: - * - `multiply(rhs, out)`, which should compute the product of the matrix with `rhs`, a `Eigen::VectorXd`-equivalent of length equal to the number of columns; - * and stores the result in `out`, an `Eigen::VectorXd` of length equal to the number of rows. - * - A `*` method where the right-hand side is an `Eigen::VectorXd` (or equivalent expression) of length equal to the number of columsn, - * and returns an `Eigen::VectorXd`-equivalent of length equal to the number of rows. + * - `multiply(rhs, out)`, which should compute the product of the matrix with `rhs`, a `VectorType`-equivalent of length equal to the number of columns; + * and stores the result in `out`, an `VectorType` of length equal to the number of rows. + * - A `*` method where the right-hand side is an `VectorType` (or equivalent expression) of length equal to the number of columsn, + * and returns an `VectorType`-equivalent of length equal to the number of rows. * - One of the following for matrix transpose-vector multiplication: - * - `adjoint_multiply(rhs, out)`, which should compute the product of the matrix transpose with `rhs`, a `Eigen::VectorXd`-equivalent of length equal to the number of rows; - * and stores the result in `out`, an `Eigen::VectorXd` of length equal to the number of columns. + * - `adjoint_multiply(rhs, out)`, which should compute the product of the matrix transpose with `rhs`, a `VectorType`-equivalent of length equal to the number of rows; + * and stores the result in `out`, an `VectorType` of length equal to the number of columns. * - An `adjoint()` method that returns an instance of any class that has a `*` method for matrix-vector multiplication. - * The method should accept an `Eigen::VectorXd`-equivalent of length equal to the number of rows, - * and return an `Eigen::VectorXd`-equvialent of length equal to the number of columns. - * - A `realize()` method that returns an `Eigen::MatrixXd` object representing the modified matrix. - * This can be omitted if an `Eigen::MatrixXd` can be copy-constructed from the class. + * The method should accept an `VectorType`-equivalent of length equal to the number of rows, + * and return an `VectorType`-equvialent of length equal to the number of columns. + * - A `realize()` method that returns an `MatrixType` object representing the modified matrix. + * This can be omitted if an `MatrixType` can be copy-constructed from the class. * * See the `Centered` and `Scaled` classes for more details. * * If the smallest dimension of `mat` is below 6, this method falls back to performing an exact SVD. */ - template + template std::pair run( - const Matrix& mat, - Eigen::MatrixXd& outU, - Eigen::MatrixXd& outV, - Eigen::VectorXd& outD, + const M& mat, + MatrixType& outU, + MatrixType& outV, + RealVectorType& outD, Engine* eng = null_rng(), - Eigen::VectorXd* init = NULL) + VectorType* init = NULL) { if (eng == NULL) { std::mt19937_64 rng(seed); @@ -316,10 +320,10 @@ class Irlba { std::pair run_internal( const M& mat, Engine& eng, - Eigen::MatrixXd& outU, - Eigen::MatrixXd& outV, - Eigen::VectorXd& outD, - Eigen::VectorXd* init) + MatrixType& outU, + MatrixType& outV, + RealVectorType& outD, + VectorType* init) { const int smaller = std::min(mat.rows(), mat.cols()); if (number >= smaller) { @@ -334,7 +338,7 @@ class Irlba { const int work = std::min(number + extra_work, smaller); - Eigen::MatrixXd V(mat.cols(), work); + MatrixType V(mat.cols(), work); if (init) { if (init->size() != V.rows()) { throw std::runtime_error("initialization vector does not have expected number of rows"); @@ -347,19 +351,19 @@ class Irlba { bool converged = false; int iter = 0, k =0; - Eigen::JacobiSVD svd(work, work, Eigen::ComputeThinU | Eigen::ComputeThinV); + Eigen::JacobiSVD svd(work, work, Eigen::ComputeThinU | Eigen::ComputeThinV); auto lptmp = lp.initialize(mat); - Eigen::MatrixXd W(mat.rows(), work); - Eigen::MatrixXd Wtmp(mat.rows(), work); - Eigen::MatrixXd Vtmp(mat.cols(), work); + MatrixType W(mat.rows(), work); + MatrixType Wtmp(mat.rows(), work); + MatrixType Vtmp(mat.cols(), work); - Eigen::MatrixXd B(work, work); + MatrixType B(work, work); B.setZero(work, work); - Eigen::VectorXd res(work); - Eigen::VectorXd F(mat.cols()); + VectorType res(work); + VectorType F(mat.cols()); - Eigen::VectorXd prevS(work); + VectorType prevS(work); for (; iter < maxit; ++iter) { // Technically, this is only a 'true' Lanczos bidiagonalization @@ -379,7 +383,7 @@ class Irlba { 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. + if (B(work-1, work-1) == Scalar(0)) { // a.k.a. the final value of 'S' from the Lanczos iterations. converged = true; break; } @@ -454,18 +458,18 @@ class Irlba { } private: - template - void exact(const M& mat, Eigen::MatrixXd& outU, Eigen::MatrixXd& outV, Eigen::VectorXd& outD) { - Eigen::BDCSVD svd(mat.rows(), mat.cols(), Eigen::ComputeThinU | Eigen::ComputeThinV); + template + void exact(const M& mat, MatrixType& outU, MatrixType& outV, RealVectorType& outD) { + Eigen::BDCSVD svd(mat.rows(), mat.cols(), Eigen::ComputeThinU | Eigen::ComputeThinV); - if constexpr(std::is_same::value) { + if constexpr(std::is_same::value) { svd.compute(mat); } else { - if constexpr(has_realize_method::value) { - Eigen::MatrixXd adjusted = mat.realize(); + if constexpr(has_realize_method::value) { + MatrixType adjusted = mat.realize(); svd.compute(adjusted); } else { - Eigen::MatrixXd adjusted(mat); + MatrixType adjusted(mat); svd.compute(adjusted); } } @@ -492,20 +496,20 @@ class Irlba { * The number of rows in `U` is equal to the number of rows in the input matrix, * and the number of columns is equal to the number of requested vectors. */ - Eigen::MatrixXd U; + MatrixType U; /** * The right singular vectors, stored as columns of `U`. * The number of rows in `U` is equal to the number of columns in the input matrix, * and the number of columns is equal to the number of requested vectors. */ - Eigen::MatrixXd V; + MatrixType V; /** * The requested number of singular values, ordered by decreasing value. * These correspond to the vectors in `U` and `V`. */ - Eigen::VectorXd D; + RealVectorType D; /** * Whether the algorithm converged. @@ -536,7 +540,7 @@ class Irlba { * @return A `Results` object containing the singular vectors and values, as well as some statistics on convergence. */ template - Results run(const M& mat, bool center, bool scale, Engine* eng = null_rng(), Eigen::VectorXd* init = NULL) { + Results run(const M& mat, bool center, bool scale, Engine* eng = null_rng(), VectorType* init = NULL) { Results output; auto stats = run(mat, center, scale, output.U, output.V, output.D, eng, init); output.converged = stats.first; @@ -560,7 +564,7 @@ class Irlba { * @return A `Results` object containing the singular vectors and values, as well as some statistics on convergence. */ template - Results run(const M& mat, Engine* eng = null_rng(), Eigen::VectorXd* init = NULL) { + Results run(const M& mat, Engine* eng = null_rng(), VectorType* init = NULL) { Results output; auto stats = run(mat, output.U, output.V, output.D, eng, init); output.converged = stats.first; diff --git a/include/irlba/lanczos.hpp b/include/irlba/lanczos.hpp index 3eab342..125d700 100644 --- a/include/irlba/lanczos.hpp +++ b/include/irlba/lanczos.hpp @@ -18,8 +18,12 @@ namespace irlba { /** * @brief Perform Lanczos bidiagonalization on an input matrix. */ +template class LanczosBidiagonalization { public: + using Scalar = typename MatrixType::Scalar; + using VectorType = Eigen::Vector; + using RealVectorType = Eigen::Vector::Real, MatrixType::RowsAtCompileTime>; struct Defaults { /** * See `set_epsilon()` for details. @@ -60,16 +64,16 @@ class LanczosBidiagonalization { * * @return Vector of residuals of length equal to the number of columns of `mat` in `run()`. */ - const Eigen::VectorXd& residuals() const { + const VectorType& residuals() const { return F; } /** * @cond */ - Eigen::VectorXd F; - Eigen::VectorXd W_next; - Eigen::VectorXd orthog_tmp; + VectorType F; + VectorType W_next; + VectorType orthog_tmp; /** * @endcond */ @@ -110,9 +114,9 @@ class LanczosBidiagonalization { template void run( const M& mat, - Eigen::MatrixXd& W, - Eigen::MatrixXd& V, - Eigen::MatrixXd& B, + MatrixType& W, + MatrixType& V, + MatrixType& B, Engine& eng, Intermediates& inter, int start = 0) @@ -125,7 +129,7 @@ class LanczosBidiagonalization { auto& otmp = inter.orthog_tmp; F = V.col(start); - if constexpr(has_multiply_method::value) { + if constexpr(has_multiply_method::value) { W_next.noalias() = mat * F; } else { mat.multiply(F, W_next); @@ -145,7 +149,7 @@ class LanczosBidiagonalization { // The Lanczos iterations themselves. for (int j = start; j < work; ++j) { - if constexpr(has_adjoint_multiply_method::value) { + if constexpr(has_adjoint_multiply_method::value) { F.noalias() = mat.adjoint() * W.col(j); } else { mat.adjoint_multiply(W.col(j), F); @@ -172,7 +176,7 @@ class LanczosBidiagonalization { B(j, j) = S; B(j, j + 1) = R_F; - if constexpr(has_multiply_method::value) { + if constexpr(has_multiply_method::value) { W_next.noalias() = mat * F; } else { mat.multiply(F, W_next); diff --git a/include/irlba/utils.hpp b/include/irlba/utils.hpp index cc5295b..c8a31b9 100644 --- a/include/irlba/utils.hpp +++ b/include/irlba/utils.hpp @@ -25,7 +25,8 @@ namespace irlba { * @return `vec` is modified to contain `vec - mat0 * t(mat0) * vec`, where `mat0` is defined as the first `ncols` columns of `mat`. * This ensures that it is orthogonal to each column of `mat0`. */ -inline void orthogonalize_vector(const Eigen::MatrixXd& mat, Eigen::VectorXd& vec, size_t ncols, Eigen::VectorXd& tmp) { +template +inline void orthogonalize_vector(const MatrixType& mat, VectorType& vec, size_t ncols, VectorType& tmp) { tmp.head(ncols).noalias() = mat.leftCols(ncols).adjoint() * vec; vec.noalias() -= mat.leftCols(ncols) * tmp.head(ncols); return; @@ -95,8 +96,11 @@ void fill_with_random_normals(Matrix& mat, int column, Engine& eng) { /** * @brief Test for IRLBA convergence. */ +template struct ConvergenceTest { public: + using Scalar = typename VectorType::Scalar; + using RealVectorType = Eigen::Vector::Real, VectorType::RowsAtCompileTime>; struct Defaults { /** * See `set_tol()` for more details. @@ -148,7 +152,7 @@ struct ConvergenceTest { * * @return The number of singular values/vectors that have achieved convergence. */ - int run(const Eigen::VectorXd& sv, const Eigen::VectorXd& residuals, const Eigen::VectorXd& last) { + int run(const RealVectorType& sv, const VectorType& residuals, const VectorType& last) { int counter = 0; double Smax = *std::max_element(sv.begin(), sv.end()); double svtol_actual = (svtol >= 0 ? svtol : tol); @@ -176,34 +180,34 @@ constexpr std::mt19937_64* null_rng() { return NULL; } /** * @cond */ -template +template struct has_multiply_method { static constexpr bool value = false; }; -template -struct has_multiply_method() * std::declval()), 0)> { +template +struct has_multiply_method() * std::declval()), 0)> { static constexpr bool value = true; }; -template +template struct has_adjoint_multiply_method { static constexpr bool value = false; }; -template -struct has_adjoint_multiply_method().adjoint() * std::declval()), 0)> { +template +struct has_adjoint_multiply_method().adjoint() * std::declval()), 0)> { static constexpr bool value = true; }; -template +template struct has_realize_method { static constexpr bool value = false; }; -template -struct has_realize_method().realize(), 0)> { - static constexpr bool value = std::is_same().realize()), Eigen::MatrixXd>::value; +template +struct has_realize_method().realize(), 0)> { + static constexpr bool value = std::is_same().realize()), M2>::value; }; /** * @endcond diff --git a/include/irlba/wrappers.hpp b/include/irlba/wrappers.hpp index f0eb717..8562889 100644 --- a/include/irlba/wrappers.hpp +++ b/include/irlba/wrappers.hpp @@ -34,14 +34,14 @@ namespace irlba { * This modification involves centering all columns, i.e., subtracting the mean of each column from the values of that column. * Naively doing such an operation would involve loss of sparsity, which we avoid by deferring the subtraction into the subspace defined by `rhs`. */ -template +template struct Centered { /** * @param m Underlying matrix to be column-centered. * @param c Vector of length equal to the number of columns of `m`, * containing the value to subtract from each column. */ - Centered(const Matrix* m, const Eigen::VectorXd* c) : mat(m), center(c) {} + Centered(const InputType* m, const VectorType* c) : mat(m), center(c) {} /** * @return Number of rows in the matrix. @@ -63,8 +63,8 @@ struct Centered { * @return `out` is filled with the product of this matrix and `rhs`. */ template - void multiply(const Right& rhs, Eigen::VectorXd& out) const { - if constexpr(has_multiply_method::value) { + void multiply(const Right& rhs, VectorType& out) const { + if constexpr(has_multiply_method::value) { out.noalias() = *mat * rhs; } else { mat->multiply(rhs, out); @@ -87,8 +87,8 @@ struct Centered { * @return `out` is filled with the product of the transpose of this matrix and `rhs`. */ template - void adjoint_multiply(const Right& rhs, Eigen::VectorXd& out) const { - if constexpr(has_adjoint_multiply_method::value) { + void adjoint_multiply(const Right& rhs, VectorType& out) const { + if constexpr(has_adjoint_multiply_method::value) { out.noalias() = mat->adjoint() * rhs; } else { mat->adjoint_multiply(rhs, out); @@ -103,8 +103,8 @@ struct Centered { * @return A realized version of the centered matrix, * where the centering has been explicitly applied. */ - Eigen::MatrixXd realize() const { - auto subtractor = [&](Eigen::MatrixXd& m) -> void { + MatrixType realize() const { + auto subtractor = [&](MatrixType& m) -> void { for (Eigen::Index c = 0; c < m.cols(); ++c) { for (Eigen::Index r = 0; r < m.rows(); ++r) { m(r, c) -= (*center)[c]; @@ -112,20 +112,20 @@ struct Centered { } }; - if constexpr(has_realize_method::value) { - Eigen::MatrixXd output = mat->realize(); + if constexpr(has_realize_method::value) { + MatrixType output = mat->realize(); subtractor(output); return output; } else { - Eigen::MatrixXd output(*mat); + MatrixType output(*mat); subtractor(output); return output; } } private: - const Matrix* mat; - const Eigen::VectorXd* center; + const InputType* mat; + const VectorType* center; }; /** @@ -136,14 +136,14 @@ struct Centered { * This modification involves scaling all columns, i.e., dividing the values of each column by the standard deviation of that column to achieve unit variance. * Naively doing such an operation would involve a copy of the matrix, which we avoid by deferring the scaling into the subspace defined by `rhs`. */ -template +template struct Scaled { /** * @param m Underlying matrix to be column-scaled. * @param s Vector of length equal to the number of columns of `m`, * containing the value to scale each column. */ - Scaled(const Matrix* m, const Eigen::VectorXd* s) : mat(m), scale(s) {} + Scaled(const InputType* m, const VectorType* s) : mat(m), scale(s) {} /** * @return Number of rows in the matrix. @@ -165,8 +165,8 @@ struct Scaled { * @return `out` is filled with the product of this matrix and `rhs`. */ template - void multiply(const Right& rhs, Eigen::VectorXd& out) const { - if constexpr(has_multiply_method::value) { + void multiply(const Right& rhs, VectorType& out) const { + if constexpr(has_multiply_method::value) { out.noalias() = *mat * rhs.cwiseQuotient(*scale); } else { mat->multiply(rhs.cwiseQuotient(*scale), out); @@ -184,8 +184,8 @@ struct Scaled { * @return `out` is filled with the product of the transpose of this matrix and `rhs`. */ template - void adjoint_multiply(const Right& rhs, Eigen::VectorXd& out) const { - if constexpr(has_adjoint_multiply_method::value) { + void adjoint_multiply(const Right& rhs, VectorType& out) const { + if constexpr(has_adjoint_multiply_method::value) { out.noalias() = mat->adjoint() * rhs; } else { mat->adjoint_multiply(rhs, out); @@ -198,8 +198,8 @@ struct Scaled { * @return A realized version of the scaled matrix, * where the scaling has been explicitly applied. */ - Eigen::MatrixXd realize() const { - auto scaler = [&](Eigen::MatrixXd& m) -> void { + MatrixType realize() const { + auto scaler = [&](MatrixType& m) -> void { for (Eigen::Index c = 0; c < m.cols(); ++c) { for (Eigen::Index r = 0; r < m.rows(); ++r) { m(r, c) /= (*scale)[c]; @@ -207,20 +207,20 @@ struct Scaled { } }; - if constexpr(has_realize_method::value) { - Eigen::MatrixXd output = mat->realize(); + if constexpr(has_realize_method::value) { + MatrixType output = mat->realize(); scaler(output); return output; } else { - Eigen::MatrixXd output(*mat); + MatrixType output(*mat); scaler(output); return output; } } private: - const Matrix* mat; - const Eigen::VectorXd* scale; + const InputType* mat; + const VectorType* scale; }; } diff --git a/tests/src/NormalSampler.h b/tests/src/NormalSampler.h index 8b1785a..8a5a83e 100644 --- a/tests/src/NormalSampler.h +++ b/tests/src/NormalSampler.h @@ -22,6 +22,17 @@ inline Eigen::MatrixXd create_random_matrix(size_t nr, size_t nc, int seed = 42) return A; } +inline Eigen::MatrixXcd create_random_complex_matrix(size_t nr, size_t nc, int seed = 42) { + Eigen::MatrixXcd A(nr, nc); + NormalSampler norm(seed); + for (size_t i = 0; i < nc; ++i) { + for (size_t j = 0; j < nr; ++j) { + A(j, i) = std::complex(norm(), norm()); + } + } + return A; +} + inline Eigen::VectorXd create_random_vector(size_t n, int seed = 50) { NormalSampler norm(seed); Eigen::VectorXd output(n); diff --git a/tests/src/compare.h b/tests/src/compare.h index 6436c38..fa79f77 100644 --- a/tests/src/compare.h +++ b/tests/src/compare.h @@ -4,12 +4,13 @@ #include #include "Eigen/Dense" -inline bool same_same(double left, double right, double tol) { +template +inline bool same_same(T left, T right, double tol) { return std::abs(left - right) <= (std::abs(left) + std::abs(right)) * tol; } -template -void expect_equal_column_vectors(const Eigen::MatrixXd& left, const Eigen::MatrixXd& right, double tol=1e-8) { +template +void expect_equal_column_vectors(const LeftType& left, const RightType& right, double tol=1e-8) { ASSERT_EQ(left.cols(), right.cols()); ASSERT_EQ(left.rows(), right.rows()); @@ -30,7 +31,8 @@ void expect_equal_column_vectors(const Eigen::MatrixXd& left, const Eigen::Matri return; } -inline void expect_equal_matrix(const Eigen::MatrixXd& left, const Eigen::MatrixXd& right, double tol=1e-8) { +template +inline void expect_equal_matrix(const LeftType& left, const RightType& right, double tol=1e-8) { ASSERT_EQ(left.cols(), right.cols()); ASSERT_EQ(left.rows(), right.rows()); for (size_t i = 0; i < left.cols(); ++i) { @@ -40,7 +42,8 @@ inline void expect_equal_matrix(const Eigen::MatrixXd& left, const Eigen::Matrix } } -inline void expect_equal_vectors(const Eigen::VectorXd& left, const Eigen::VectorXd& right, double tol=1e-8) { +template +inline void expect_equal_vectors(const LeftType& left, const RightType& right, double tol=1e-8) { ASSERT_EQ(left.size(), right.size()); for (size_t i = 0; i < left.size(); ++i) { EXPECT_TRUE(same_same(left[i], right[i], tol)); diff --git a/tests/src/irlba.cpp b/tests/src/irlba.cpp index af25943..97f4207 100644 --- a/tests/src/irlba.cpp +++ b/tests/src/irlba.cpp @@ -25,6 +25,22 @@ TEST(IrlbaTest, Exact) { expect_equal_column_vectors(res.V, svd.matrixV().leftCols(rank), 1e-8); } +TEST(IrlbaComplexTest, Exact) { + // For the test, the key is that rank + workspace > min(nr, nc), in which + // case we can be pretty confident of getting a near-exact match of the + // true SVD. Otherwise it's more approximate and the test is weaker. + int rank = 5; + auto A = create_random_complex_matrix(20, 10); + + irlba::Irlba irb; + auto res = irb.set_number(rank).run(A); + + Eigen::BDCSVD svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV); + expect_equal_vectors(res.D, svd.singularValues().head(rank), 1e-8); + expect_equal_column_vectors(res.U, svd.matrixU().leftCols(rank), 1e-8); + expect_equal_column_vectors(res.V, svd.matrixV().leftCols(rank), 1e-8); +} + class IrlbaTester : public ::testing::TestWithParam > { protected: template diff --git a/tests/src/utils.cpp b/tests/src/utils.cpp index 39095b0..44af712 100644 --- a/tests/src/utils.cpp +++ b/tests/src/utils.cpp @@ -1,5 +1,6 @@ #include #include "irlba/utils.hpp" +#include "irlba/wrappers.hpp" #include "Eigen/Dense" #include "Eigen/Sparse" #include "NormalSampler.h" @@ -98,9 +99,14 @@ struct dummy_class { }; TEST(UtilsTest, ConvertibleCheck) { - EXPECT_FALSE(irlba::has_realize_method::value); - EXPECT_FALSE(irlba::has_realize_method >::value); - EXPECT_TRUE(irlba::has_realize_method::value); + bool const realize_xd_xd = irlba::has_realize_method::value; + EXPECT_FALSE(realize_xd_xd); + bool const realize_sp_sp = irlba::has_realize_method, Eigen::SparseMatrix>::value; + EXPECT_FALSE(realize_sp_sp); + bool const realize_dummy = irlba::has_realize_method::value; + EXPECT_TRUE(realize_dummy); + bool const realize_scaled = irlba::has_realize_method, Eigen::MatrixXd>::value; + EXPECT_TRUE(realize_scaled); } TEST(UtilsTest, NormalSampler) { diff --git a/tests/src/wrappers.cpp b/tests/src/wrappers.cpp index 9744aca..b420fe3 100644 --- a/tests/src/wrappers.cpp +++ b/tests/src/wrappers.cpp @@ -17,14 +17,14 @@ TEST(WrapperTest, Centering) { auto C = create_random_vector(10, 1234); Eigen::VectorXd output(20); centered.multiply(C, output); - expect_equal_matrix(ref * C, output); + expect_equal_matrix(ref * C, output); } { auto C = create_random_vector(20, 1234); Eigen::VectorXd output(10); centered.adjoint_multiply(C, output); - expect_equal_matrix(ref.adjoint() * C, output); + expect_equal_matrix(ref.adjoint() * C, output); } } @@ -41,13 +41,13 @@ TEST(WrapperTest, Scaling) { auto C = create_random_vector(10, 1234); Eigen::VectorXd output(20); centered.multiply(C, output); - expect_equal_matrix(ref * C, output); + expect_equal_matrix(ref * C, output); } { auto C = create_random_vector(20, 1234); Eigen::VectorXd output(10); centered.adjoint_multiply(C, output); - expect_equal_matrix(ref.adjoint() * C, output); + expect_equal_matrix(ref.adjoint() * C, output); } }