Protect against zero-extent matrices, which cause Eigen's SVD to abort.

Author: LTLA

CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.14...3.25)
 
 project(irlba
-    VERSION 3.0.0
+    VERSION 3.0.1
     DESCRIPTION "C++ port of the IRLBA algorithm"
     LANGUAGES CXX)
 

include/irlba/compute.hpp

@@ -159,6 +159,15 @@ std::pair<bool, int> compute(
         throw std::runtime_error("requested number of singular values must be less than the smaller matrix dimension for IRLBA iterations");
     }
 
+    // Optimization if no SVs are requested. Also protect the exact SVD code
+    // from a zero-extent matrix, as this causes Eigen to shit itself.
+    if (requested_number == 0) {
+        outD.resize(requested_number);
+        outU.resize(matrix.rows(), requested_number);
+        outV.resize(matrix.cols(), requested_number);
+        return std::make_pair(true, 0);
+    }
+
     // 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 (
@@ -169,10 +178,9 @@ std::pair<bool, int> compute(
         return std::make_pair(true, 0);
     }
 
-    Eigen::Index work = choose_requested_plus_extra_work_or_smaller(requested_number, options.extra_work, smaller);
-    if (work == 0) {
-        throw std::runtime_error("number of requested dimensions must be positive");
-    }
+    // We know work must be positive at this point, as we would have returned early if requested_number = 0.
+    // Similar, if smaller = 0, we would have either thrown earlier or set requested_number = 0.
+    const Eigen::Index work = choose_requested_plus_extra_work_or_smaller(requested_number, options.extra_work, smaller);
 
     // Don't worry about sanitizing dimensions for Eigen constructors,
     // as the former are stored as Eigen::Index and the latter accepts Eigen::Index inputs.

tests/src/compute.cpp

@@ -22,7 +22,7 @@ TEST(Compute, CompareToExact) {
     irlba::Options opt;
     opt.exact_for_large_number = false;
     irlba::SimpleMatrix<Eigen::VectorXd, Eigen::MatrixXd, decltype(&A)> wrapped(&A);
-    auto res = irlba::compute(wrapped, 5, opt);
+    auto res = irlba::compute(wrapped, rank, opt);
 
     Eigen::JacobiSVD<decltype(A), Eigen::ComputeThinU | Eigen::ComputeThinV> svd(A);
     expect_equal_vectors(res.D, svd.singularValues().head(rank), 1e-8);
@@ -49,6 +49,35 @@ TEST(Compute, ChooseRequestedPlusExtraWorkOrSmaller) {
     EXPECT_EQ(irlba::choose_requested_plus_extra_work_or_smaller(10, 20, 35), 30);
 }
 
+TEST(Compute, ZeroExtent) {
+    irlba::Options opt;
+    opt.cap_number = true;
+
+    {
+        auto A = create_random_matrix(0, 10);
+        irlba::SimpleMatrix<Eigen::VectorXd, Eigen::MatrixXd, decltype(&A)> wrapped(&A);
+        auto res = irlba::compute(wrapped, 5, opt);
+
+        EXPECT_EQ(res.D.size(), 0);
+        EXPECT_EQ(res.U.rows(), 0);
+        EXPECT_EQ(res.U.cols(), 0);
+        EXPECT_EQ(res.V.rows(), 10);
+        EXPECT_EQ(res.V.cols(), 0);
+    }
+
+    {
+        auto A = create_random_matrix(10, 0);
+        irlba::SimpleMatrix<Eigen::VectorXd, Eigen::MatrixXd, decltype(&A)> wrapped(&A);
+        auto res = irlba::compute(wrapped, 5, opt);
+
+        EXPECT_EQ(res.D.size(), 0);
+        EXPECT_EQ(res.U.rows(), 10);
+        EXPECT_EQ(res.U.cols(), 0);
+        EXPECT_EQ(res.V.rows(), 0);
+        EXPECT_EQ(res.V.cols(), 0);
+    }
+}
+
 TEST(Compute, UpdateK) {
     auto ref_update_k = [](int k, const int requested_number, const int n_converged, const int work) -> int {
         if (k < requested_number + n_converged) {
@@ -214,17 +243,6 @@ TEST(Compute, Fails) {
     }
     EXPECT_TRUE(message.find("must be less than") != std::string::npos);
 
-    // Requested number of SVs + extra work is zero.
-    message.clear();
-    try {
-        irlba::Options opt;
-        opt.extra_work = 0;
-        irlba::compute(wrapped, 0, opt);
-    } catch (const std::exception& e) {
-        message = e.what();
-    }
-    EXPECT_TRUE(message.find("must be positive") != std::string::npos);
-
     // Initialization vector is not of the right length.
     message.clear();
     try {