Switch to modified Gram-Schmidt for orthogonalization.

This is actually faster as we can overwrite 'vec' as we go along, rather than
having to store the projections in a separate allocation to avoid mutating it.

Author: LTLA

include/irlba/lanczos.hpp

@@ -15,24 +15,28 @@
 namespace irlba {
 
 template<class EigenMatrix_, class EigenVector_>
-void orthogonalize_vector(const EigenMatrix_& mat, EigenVector_& vec, size_t ncols, EigenVector_& tmp) {
-    tmp.head(ncols).noalias() = mat.leftCols(ncols).adjoint() * vec;
-    vec.noalias() -= mat.leftCols(ncols) * tmp.head(ncols);
+void orthogonalize_vector(const EigenMatrix_& mat, EigenVector_& vec, Eigen::Index ncols) {
+    // Original package uses classical Gram Schmidt but modified GS is actually faster when we
+    // sum the timing across all orthogonalization steps in the Lanczos process.
+    // Possibly because it allows us to modify 'vec' for each of mat's progressive column vectors, 
+    // rather than requiring a temporary space to store the projections from the original 'vec'.
+    for (Eigen::Index c = 0; c < ncols; ++c) {
+        // No need to divide by the norm of mat.col(c), as this is always 1.
+        vec -= mat.col(c) * vec.dot(mat.col(c)); 
+    }
 }
 
 template<class EigenVector_, class Matrix_>
 struct LanczosWorkspace {
     LanczosWorkspace(const Matrix_& mat) : 
         F(mat.cols()), 
         W_next(mat.rows()), 
-        orthog_tmp(mat.cols()), 
         work(mat.new_known_workspace()),
         awork(mat.new_known_adjoint_workspace())
     {}
 
     EigenVector_ F; 
     EigenVector_ W_next;
-    EigenVector_ orthog_tmp;
 
     I<decltype(std::declval<Matrix_>().new_known_workspace())> work;
     I<decltype(std::declval<Matrix_>().new_known_adjoint_workspace())> awork;
@@ -68,15 +72,14 @@ int run_lanczos_bidiagonalization(
     const Eigen::Index work = W.cols();
     auto& F = inter.F;
     auto& W_next = inter.W_next;
-    auto& otmp = inter.orthog_tmp;
 
     F = V.col(start);
     inter.work->multiply(F, W_next); // i.e., W_next = mat * F;
     int mult = 1;
 
     // If start = 0, there's nothing to orthogonalize against.
     if (start) {
-        orthogonalize_vector(W, W_next, start, otmp);
+        orthogonalize_vector(W, W_next, start);
     }
 
     Float S = W_next.norm();
@@ -95,14 +98,14 @@ int run_lanczos_bidiagonalization(
         ++mult;
 
         F -= S * V.col(j); // equivalent to daxpy.
-        orthogonalize_vector(V, F, j + 1, otmp);
+        orthogonalize_vector(V, F, j + 1);
 
         if (j + 1 < work) {
             Float R_F = F.norm();
 
             if (R_F < eps) {
                 fill_with_random_normals(F, eng);
-                orthogonalize_vector(V, F, j + 1, otmp);
+                orthogonalize_vector(V, F, j + 1);
                 R_F = F.norm();
                 F /= R_F;
                 R_F = 0;
@@ -121,12 +124,12 @@ int run_lanczos_bidiagonalization(
             // Full re-orthogonalization, using the left-most 'j + 1' columns of W.
             // Recall that W_next will be the 'j + 2'-th column, i.e., W.col(j + 1) in
             // 0-indexed terms, so we want to orthogonalize to all previous columns.
-            orthogonalize_vector(W, W_next, j + 1, otmp);
+            orthogonalize_vector(W, W_next, j + 1);
 
             S = W_next.norm();
             if (S < eps) {
                 fill_with_random_normals(W_next, eng);
-                orthogonalize_vector(W, W_next, j + 1, otmp);
+                orthogonalize_vector(W, W_next, j + 1);
                 S = W_next.norm();
                 W_next /= S;
                 S = 0;

tests/src/lanczos.cpp

@@ -17,37 +17,28 @@ TEST(Lanczos, Orthogonalize) {
          -0.1683330,  0.1634341, -0.75523159, -0.5473305105,
           0.8749166, -0.2434161, -0.09989224, -0.0004957394;
 
-    auto orthogonalize_vector = [](const Eigen::MatrixXd& mat, Eigen::VectorXd& vec, size_t ncols) -> void {
-        Eigen::VectorXd tmp(mat.cols());
-        irlba::orthogonalize_vector(mat, vec, ncols, tmp);
-    };
-
     // Generated by rnorm(5)
     Eigen::VectorXd v(5);
     v << -0.24054848, -0.04785069, -0.76491749, -0.65634291, 0.62815141;
 
     auto copy = v;
-    orthogonalize_vector(m, copy, 4);
+    irlba::orthogonalize_vector(m, copy, 4);
 
     // Checking that the L2 norm is not all-zero.
     double l2 = 0;
     for (auto x : copy) { l2 += x*x; }
-    EXPECT_TRUE(l2 > 0.1);
+    EXPECT_GT(l2, 0.1);
 
     // Checking that we do have orthogonality.
     for (size_t i = 0; i < 4; ++i) {
         auto col = m.col(i);
-        auto cit = col.begin();
-        double sum = 0;
-        for (auto x : copy) {
-            sum += x * (*cit);
-        }
-        EXPECT_TRUE(sum < 0.0000000001);
+        double sum = std::inner_product(col.begin(), col.end(), copy.begin(), 0.0);
+        EXPECT_LT(std::abs(sum), 0.000001);
     }
 
     // Checking that the specification of columns has an effect.
     auto copy2 = v;
-    orthogonalize_vector(m, copy, 2);
+    irlba::orthogonalize_vector(m, copy2, 2);
     EXPECT_NE(copy, copy2);
 }