Switch back to an exact SVD if too many singular values are requested.

include/irlba/irlba.hpp

@@ -49,6 +49,16 @@ class Irlba {
          */
         static constexpr int maxit = 1000;
 
+        /**
+         * See `set_exact_for_small_matrix()` for more details.
+         */
+        static constexpr bool exact_for_small_matrix = true;
+
+        /**
+         * See `set_exact_for_large_number()` for more details.
+         */
+        static constexpr bool exact_for_large_number = true;
+
         /**
          * See `set_seed()` for more details.
          */
@@ -62,6 +72,9 @@ class Irlba {
     int maxit = Defaults::maxit;
     uint64_t seed = Defaults::seed;
 
+    bool exact_for_small_matrix = Defaults::exact_for_small_matrix;
+    bool exact_for_large_number = Defaults::exact_for_large_number;
+
     ConvergenceTest convtest;
 
 public:
@@ -152,6 +165,28 @@ class Irlba {
         return *this;
     }
 
+    /**
+     * @param s Whether to perform an exact SVD if the matrix is too small (fewer than 6 elements in any dimension).
+     * This is more efficient and avoids inaccuracies from an insufficient workspace.
+     *
+     * @return A reference to the `Irlba` instance.
+     */
+    Irlba& set_exact_for_small_matrix(bool s = Defaults::exact_for_small_matrix) {
+        exact_for_small_matrix = s;
+        return *this;
+    }
+
+    /**
+     * @param s Whether to perform an exact SVD if the requested number of singular values is too large (greater than or equal to half the smaller matrix dimension).
+     * This is more efficient and avoids inaccuracies from an insufficient workspace.
+     *
+     * @return A reference to the `Irlba` instance.
+     */
+    Irlba& set_exact_for_large_number(bool s = Defaults::exact_for_large_number) {
+        exact_for_large_number = s;
+        return *this;
+    }
+
 public:
     /** 
      * Run IRLBA on an input matrix to perform an approximate SVD, with arbitrary centering and scaling operations.
@@ -307,12 +342,15 @@ class Irlba {
         Eigen::VectorXd* init)
     const {
         const int smaller = std::min(mat.rows(), mat.cols());
-        if (number >= smaller) {
-            throw std::runtime_error("requested number of singular values must be less than smaller matrix dimension");
+        if (number > smaller) {
+            throw std::runtime_error("requested number of singular values cannot be greater than the smaller matrix dimension");
+        } else if (number == smaller && !exact_for_large_number) {
+            throw std::runtime_error("requested number of singular values must be less than the smaller matrix dimension");
         }
 
-        // Falling back to an exact SVD for small matrices.
-        if (smaller < 6) {
+        // 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 ((exact_for_small_matrix && smaller < 6) || (exact_for_large_number && number * 2 >= smaller)) {
             exact(mat, outU, outV, outD);
             return std::make_pair(true, 0);
         }

tests/src/irlba.cpp

@@ -9,12 +9,12 @@
 #include "Eigen/Dense"
 #include <random>
 
-TEST(IrlbaTest, Exact) {
+TEST(IrlbaTest, CompareToExact) {
     // 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_matrix(20, 10);
+    auto A = create_random_matrix(50, 20);
 
     irlba::Irlba irb;
     auto res = irb.set_number(rank).run(A);
@@ -138,13 +138,31 @@ TEST(IrlbaTest, SmallExact) {
 
     irlba::Irlba irb;
     auto res = irb.set_number(2).run(small);
-      
+
     Eigen::BDCSVD svd(small, Eigen::ComputeThinU | Eigen::ComputeThinV);
     EXPECT_EQ(svd.singularValues().head(2), res.D);
     EXPECT_EQ(svd.matrixU().leftCols(2), res.U);
     EXPECT_EQ(svd.matrixV().leftCols(2), res.V);
 }
 
+TEST(IrlbaTest, LargeExact) {
+    Eigen::MatrixXd mat = create_random_matrix(20, 50);
+
+    irlba::Irlba irb;
+    auto res = irb.set_number(15).run(mat);
+
+    Eigen::BDCSVD svd(mat, Eigen::ComputeThinU | Eigen::ComputeThinV);
+    EXPECT_EQ(svd.singularValues().head(15), res.D);
+    EXPECT_EQ(svd.matrixU().leftCols(15), res.U);
+    EXPECT_EQ(svd.matrixV().leftCols(15), res.V);
+
+    // Works with the maximum number.
+    auto res2 = irb.set_number(20).run(mat);
+    EXPECT_EQ(svd.singularValues(), res2.D);
+    EXPECT_EQ(svd.matrixU(), res2.U);
+    EXPECT_EQ(svd.matrixV(), res2.V);
+}
+
 TEST(IrlbaTest, SmallExactCenterScale) {
     Eigen::MatrixXd small = create_random_matrix(10, 3);
 
@@ -200,7 +218,15 @@ TEST(IrlbaTester, Fails) {
         irb.set_number(100).run(A);
     } catch (const std::exception& e) {
         std::string message(e.what());
-        EXPECT_EQ(message.find("requested"), 0);
+        EXPECT_TRUE(message.find("must be greater than") != std::string::npos);
+    }
+
+    // Requested number of SVs > smaller dimension of the matrix.
+    try {
+        irb.set_exact_for_large_number(false).set_number(10).run(A);
+    } catch (const std::exception& e) {
+        std::string message(e.what());
+        EXPECT_TRUE(message.find("must be less than") != std::string::npos);
     }
 
     // Initialization vector is not of the right length.