BDCSVD: improve accuracy and workspace use

libeigen/eigen!2702

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/src/SVD/BDCSVDImpl.h

@@ -113,7 +113,13 @@ void bdcsvd_impl<RealScalar_>::allocate(Index diagSize, bool compU, bool compV)
 
   if (m_compV) m_naiveV = MatrixXr::Zero(diagSize, diagSize);
 
-  m_workspace.resize((diagSize + 1) * (diagSize + 1) * 3);
+  // Vector updates need the three matrix-sized packing buffers used by
+  // structured_update(). Values-only decompositions only need five vectors:
+  // diag, shifts, mus, zhat, and diagShifted.
+  if (m_compU || m_compV)
+    m_workspace.resize((diagSize + 1) * (diagSize + 1) * 3);
+  else
+    m_workspace.resize(5 * diagSize);
   m_workspaceI.resize(3 * diagSize);
 }
 
@@ -193,8 +199,6 @@ void bdcsvd_impl<RealScalar_>::computeBaseCase(SVDType& svd, Index n, Index firs
 template <typename RealScalar_>
 void bdcsvd_impl<RealScalar_>::divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift) {
   // requires rows = cols + 1;
-  using std::abs;
-  using std::sqrt;
   const Index n = lastCol - firstCol + 1;
   const Index k = n / 2;
   const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
@@ -202,7 +206,6 @@ void bdcsvd_impl<RealScalar_>::divide(Index firstCol, Index lastCol, Index first
   RealScalar betaK;
   RealScalar r0;
   RealScalar lambda, phi, c0, s0;
-  VectorType l, f;
   // We use the other algorithm which is more efficient for small
   // matrices.
   if (n < m_algoswap) {
@@ -231,14 +234,10 @@ void bdcsvd_impl<RealScalar_>::divide(Index firstCol, Index lastCol, Index first
     lambda = m_naiveU(1, firstCol + k);
     phi = m_naiveU(0, lastCol + 1);
   }
-  r0 = sqrt((abs(alphaK * lambda) * abs(alphaK * lambda)) + abs(betaK * phi) * abs(betaK * phi));
-  if (m_compU) {
-    l = m_naiveU.row(firstCol + k).segment(firstCol, k);
-    f = m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1);
-  } else {
-    l = m_naiveU.row(1).segment(firstCol, k);
-    f = m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1);
-  }
+  // LAPACK's xLASD2 likewise uses xLAPY2 for this merge coupling. The
+  // scaled hypotenuse avoids destructive underflow when both products are
+  // below sqrt(min()).
+  r0 = numext::hypot(alphaK * lambda, betaK * phi);
   if (m_compV) m_naiveV(firstRowW + k, firstColW) = Literal(1);
   if (r0 < considerZero) {
     c0 = Literal(1);
@@ -248,8 +247,22 @@ void bdcsvd_impl<RealScalar_>::divide(Index firstCol, Index lastCol, Index first
     s0 = betaK * phi / r0;
   }
 
+  m_computed(firstCol + shift, firstCol + shift) = r0;
   if (m_compU) {
-    MatrixXr q1(m_naiveU.col(firstCol + k).segment(firstCol, k + 1));
+    m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) =
+        alphaK * m_naiveU.row(firstCol + k).segment(firstCol, k).transpose();
+    m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) =
+        betaK * m_naiveU.row(firstCol + k + 1).segment(firstCol + k + 1, n - k - 1).transpose();
+  } else {
+    m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) =
+        alphaK * m_naiveU.row(1).segment(firstCol, k).transpose();
+    m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) =
+        betaK * m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).transpose();
+  }
+
+  if (m_compU) {
+    Map<VectorType, Aligned> q1(m_workspace.data(), k + 1);
+    q1 = m_naiveU.col(firstCol + k).segment(firstCol, k + 1);
     // we shift Q1 to the right
     for (Index i = firstCol + k - 1; i >= firstCol; i--)
       m_naiveU.col(i + 1).segment(firstCol, k + 1) = m_naiveU.col(i).segment(firstCol, k + 1);
@@ -278,10 +291,6 @@ void bdcsvd_impl<RealScalar_>::divide(Index firstCol, Index lastCol, Index first
     m_naiveU.row(0).segment(firstCol + k + 1, n - k - 1).setZero();
   }
 
-  m_computed(firstCol + shift, firstCol + shift) = r0;
-  m_computed.col(firstCol + shift).segment(firstCol + shift + 1, k) = alphaK * l.transpose();
-  m_computed.col(firstCol + shift).segment(firstCol + shift + k + 1, n - k - 1) = betaK * f.transpose();
-
   // Second part: try to deflate singular values in combined matrix
   deflation(firstCol, lastCol, k, firstRowW, firstColW, shift);
 
@@ -300,8 +309,10 @@ void bdcsvd_impl<RealScalar_>::divide(Index firstCol, Index lastCol, Index first
 
   if (m_compV) structured_update(m_naiveV.block(firstRowW, firstColW, n, n), VofSVD, (n + 1) / 2);
 
-  m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero();
-  m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals;
+  // Recursive children leave this block diagonal; this merge only adds its
+  // first column. Clear that column instead of rewriting the full n-by-n block.
+  m_computed.col(firstCol + shift).segment(firstCol + shift, n).setZero();
+  m_computed.diagonal().segment(firstCol + shift, n) = singVals;
 }  // end divide
 
 // Compute SVD of m_computed.block(firstCol, firstCol, n + 1, n); this block only has non-zeros in
@@ -385,8 +396,10 @@ typename bdcsvd_impl<RealScalar_>::RealScalar bdcsvd_impl<RealScalar_>::secularE
 template <typename RealScalar_>
 void bdcsvd_impl<RealScalar_>::computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm,
                                                VectorType& singVals, ArrayRef shifts, ArrayRef mus) {
+  // See Ren-Cang Li, "Solving Secular Equations Stably and Efficiently",
+  // LAPACK Working Note 89 (1994), and LAPACK's xLASD4/xLASD5 for the
+  // stability rationale behind pole-relative shifts and safeguarded steps.
   using std::abs;
-  using std::sqrt;
   using std::swap;
 
   Index n = col0.size();
@@ -499,9 +512,9 @@ void bdcsvd_impl<RealScalar_>::computeSingVals(const ArrayRef& col0, const Array
       if (numext::equal_strict(shift, left)) {
         // to avoid overflow, we must have mu > max(real_min, |z(k)|/sqrt(real_max)),
         // the factor 2 is to be more conservative
-        leftShifted =
-            numext::maxi<RealScalar>((std::numeric_limits<RealScalar>::min)(),
-                                     Literal(2) * abs(col0(k)) / sqrt((std::numeric_limits<RealScalar>::max)()));
+        leftShifted = numext::maxi<RealScalar>(
+            (std::numeric_limits<RealScalar>::min)(),
+            Literal(2) * abs(col0(k)) / numext::sqrt((std::numeric_limits<RealScalar>::max)()));
 
         // check that we did it right:
         eigen_internal_assert(
@@ -512,8 +525,9 @@ void bdcsvd_impl<RealScalar_>::computeSingVals(const ArrayRef& col0, const Array
       } else {
         leftShifted = -(right - left) * RealScalar(0.51);
         if (k + 1 < n)
-          rightShifted = -numext::maxi<RealScalar>((std::numeric_limits<RealScalar>::min)(),
-                                                   abs(col0(k + 1)) / sqrt((std::numeric_limits<RealScalar>::max)()));
+          rightShifted =
+              -numext::maxi<RealScalar>((std::numeric_limits<RealScalar>::min)(),
+                                        abs(col0(k + 1)) / numext::sqrt((std::numeric_limits<RealScalar>::max)()));
         else
           rightShifted = -(std::numeric_limits<RealScalar>::min)();
       }
@@ -557,7 +571,7 @@ template <typename RealScalar_>
 void bdcsvd_impl<RealScalar_>::perturbCol0(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm,
                                            const VectorType& singVals, const ArrayRef& shifts, const ArrayRef& mus,
                                            ArrayRef zhat) {
-  using std::sqrt;
+  using std::abs;
   Index n = col0.size();
   Index m = perm.size();
   if (m == 0) {
@@ -595,7 +609,9 @@ void bdcsvd_impl<RealScalar_>::perturbCol0(const ArrayRef& col0, const ArrayRef&
           prod *= ((singVals(j) + dk) / ((diag(i) + dk))) * ((mus(j) + diff) / ((diag(i) - dk)));
         }
       }
-      RealScalar tmp = sqrt(prod);
+      // This product is non-negative in exact arithmetic. As in LAPACK's
+      // xLASD8, take abs before sqrt to tolerate a negative rounding residue.
+      RealScalar tmp = numext::sqrt(abs(prod));
       zhat(k) = col0(k) > Literal(0) ? RealScalar(tmp) : RealScalar(-tmp);
     }
   }
@@ -615,6 +631,7 @@ void bdcsvd_impl<RealScalar_>::computeSingVecs(const ArrayRef& zhat, const Array
       if (m_compV) V.col(k) = VectorType::Unit(n, k);
     } else {
       U.col(k).setZero();
+      if (m_compV) V.col(k).setZero();
       for (Index l = 0; l < m; ++l) {
         Index i = perm(l);
         RealScalar diff = diag(i) - shifts(k);
@@ -622,22 +639,18 @@ void bdcsvd_impl<RealScalar_>::computeSingVecs(const ArrayRef& zhat, const Array
         diff -= mus(k);
         EIGEN_OPTIMIZATION_BARRIER(diff)
         U(i, k) = zhat(i) / diff / ((diag(i) + singVals[k]));
+        if (m_compV && l > 0) V(i, k) = diag(i) * zhat(i) / diff / ((diag(i) + singVals[k]));
       }
       U(n, k) = Literal(0);
-      U.col(k).normalize();
+      // LAPACK's xLASD3 normalizes these vectors with xNRM2. Use the scaled
+      // normalization unconditionally: under -ffast-math, compilers may
+      // assume that the overflowing result of norm() is finite and discard
+      // an isfinite-based fallback.
+      U.col(k).stableNormalize();
 
       if (m_compV) {
-        V.col(k).setZero();
-        for (Index l = 1; l < m; ++l) {
-          Index i = perm(l);
-          RealScalar diff = diag(i) - shifts(k);
-          EIGEN_OPTIMIZATION_BARRIER(diff)
-          diff -= mus(k);
-          EIGEN_OPTIMIZATION_BARRIER(diff)
-          V(i, k) = diag(i) * zhat(i) / diff / ((diag(i) + singVals[k]));
-        }
         V(0, k) = Literal(-1);
-        V.col(k).normalize();
+        V.col(k).stableNormalize();
       }
     }
   }
@@ -649,8 +662,6 @@ void bdcsvd_impl<RealScalar_>::computeSingVecs(const ArrayRef& zhat, const Array
 // We use a rotation to zero out zi applied to the left of M, and set di = 0.
 template <typename RealScalar_>
 void bdcsvd_impl<RealScalar_>::deflation43(Index firstCol, Index shift, Index i, Index size) {
-  using std::abs;
-  using std::sqrt;
   Index start = firstCol + shift;
   RealScalar c = m_computed(start, start);
   RealScalar s = m_computed(start + i, start);
@@ -676,9 +687,6 @@ void bdcsvd_impl<RealScalar_>::deflation43(Index firstCol, Index shift, Index i,
 template <typename RealScalar_>
 void bdcsvd_impl<RealScalar_>::deflation44(Index firstColu, Index firstColm, Index firstRowW, Index firstColW, Index i,
                                            Index j, Index size) {
-  using std::abs;
-  using std::sqrt;
-
   RealScalar s = m_computed(firstColm + i, firstColm);
   RealScalar c = m_computed(firstColm + j, firstColm);
   RealScalar r = numext::hypot(c, s);
@@ -705,7 +713,6 @@ template <typename RealScalar_>
 void bdcsvd_impl<RealScalar_>::deflation(Index firstCol, Index lastCol, Index k, Index firstRowW, Index firstColW,
                                          Index shift) {
   using std::abs;
-  using std::sqrt;
   const Index length = lastCol + 1 - firstCol;
 
   Block<MatrixXr, Dynamic, 1> col0(m_computed, firstCol + shift, firstCol + shift, length, 1);

benchmarks/SVD/bench_svd.cpp

@@ -49,6 +49,19 @@ static void BM_BDCSVD(benchmark::State& state) {
   state.SetItemsProcessed(state.iterations());
 }
 
+template <typename Scalar, int Options>
+static void BM_BDCSVDBidiagonal(benchmark::State& state) {
+  const Index size = state.range(0);
+  Matrix<Scalar, Dynamic, 1> diagonal = Matrix<Scalar, Dynamic, 1>::Random(size);
+  Matrix<Scalar, Dynamic, 1> superdiagonal = Matrix<Scalar, Dynamic, 1>::Random(size - 1);
+  BDCSVD<Mat<Scalar>, Options> svd(size, size);
+  for (auto _ : state) {
+    svd.compute(diagonal, superdiagonal);
+    benchmark::DoNotOptimize(svd.singularValues().data());
+  }
+  state.SetItemsProcessed(state.iterations());
+}
+
 // ---------- Size configurations ----------
 
 // ---------- Register benchmarks ----------
@@ -67,6 +80,9 @@ static void BM_BDCSVD(benchmark::State& state) {
     ->Args({100, 4})->Args({1000, 4})->Args({1000, 10})->Args({1000, 100}) \
     ->Args({10000, 10})->Args({10000, 100})
 
+// Direct bidiagonal input isolates the divide-and-conquer phase.
+#define BDC_BIDIAG_SIZES ->Arg(16)->Arg(32)->Arg(64)->Arg(128)->Arg(256)->Arg(512)->Arg(1024)
+
 // JacobiSVD — float
 BENCHMARK(BM_JacobiSVD<float, ComputeThinU | ComputeThinV>) JACOBI_SIZES ->Name("JacobiSVD_float_ThinUV");
 BENCHMARK(BM_JacobiSVD<float, 0>) JACOBI_SIZES ->Name("JacobiSVD_float_ValuesOnly");
@@ -83,8 +99,13 @@ BENCHMARK(BM_BDCSVD<float, 0>) BDC_SIZES ->Name("BDCSVD_float_ValuesOnly");
 BENCHMARK(BM_BDCSVD<double, ComputeThinU | ComputeThinV>) BDC_SIZES ->Name("BDCSVD_double_ThinUV");
 BENCHMARK(BM_BDCSVD<double, 0>) BDC_SIZES ->Name("BDCSVD_double_ValuesOnly");
 
+BENCHMARK(BM_BDCSVDBidiagonal<double, ComputeThinU | ComputeThinV>)
+    BDC_BIDIAG_SIZES ->Name("BDCSVD_Bidiagonal_double_ThinUV");
+BENCHMARK(BM_BDCSVDBidiagonal<double, 0>) BDC_BIDIAG_SIZES ->Name("BDCSVD_Bidiagonal_double_ValuesOnly");
+
 #undef JACOBI_SIZES
 #undef BDC_SIZES
+#undef BDC_BIDIAG_SIZES
 // clang-format on
 
 // JacobiSVD — QR preconditioner comparison (double, 64x64, ThinUV)

cmake/EigenSmokeTestList.cmake

@@ -16,6 +16,7 @@ set(ei_smoke_test_list
   autodiff_scalar_1
   bandmatrix
   bdcsvd_9
+  bdcsvd_fastmath
   bessel_functions_1
   bfloat16_float
   blasutil_1

test/CMakeLists.txt

@@ -385,6 +385,7 @@ if(COMPILER_SUPPORT_WNO_NAN_INF)
 endif()
 
 ei_add_test(fastmath "${EIGEN_FASTMATH_FLAGS}")
+ei_add_test(bdcsvd_fastmath "${EIGEN_FASTMATH_FLAGS}")
 
 # # ei_add_test(denseLM)
 

test/bdcsvd.cpp

@@ -188,23 +188,54 @@ void bdcsvd_mixed_option_enum_regression() {
 #endif
 
 #if defined(EIGEN_TEST_PART_53) || defined(EIGEN_TEST_PART_ALL)
-void bdcsvd_public_missing_predecessor() {
-  Matrix<double, 6, 6> matrix = Matrix<double, 6, 6>::Zero();
+void bdcsvd_extreme_scale_regressions() {
+  typedef Matrix<double, 6, 6> Matrix6d;
+  const double kTolerance = 16 * Matrix6d::RowsAtCompileTime * NumTraits<double>::epsilon();
+
+  const auto verify_decomposition = [kTolerance](const Matrix6d& matrix) {
+    BDCSVD<Matrix6d, ComputeFullU | ComputeFullV> svd;
+    svd.setSwitchSize(3);
+    svd.compute(matrix);
+
+    VERIFY(svd.info() == Success);
+
+    const Matrix6d reconstruction = svd.matrixU() * svd.singularValues().asDiagonal() * svd.matrixV().transpose();
+    VERIFY((reconstruction - matrix).stableNorm() <= kTolerance * matrix.stableNorm());
+
+    const Matrix6d identity = Matrix6d::Identity();
+    VERIFY((svd.matrixU().transpose() * svd.matrixU() - identity).stableNorm() <= kTolerance);
+    VERIFY((svd.matrixV().transpose() * svd.matrixV() - identity).stableNorm() <= kTolerance);
+
+    // Also exercise the values-only path, which uses a compact m_naiveU and
+    // linear workspace during divide-and-conquer merges.
+    BDCSVD<Matrix6d> valuesOnlySvd;
+    valuesOnlySvd.setSwitchSize(3);
+    valuesOnlySvd.compute(matrix);
+    VERIFY(valuesOnlySvd.info() == Success);
+    VERIFY_IS_APPROX(valuesOnlySvd.singularValues(), svd.singularValues());
+  };
+
+  Matrix6d matrix = Matrix6d::Zero();
   const double kSubnormal1040 = std::numeric_limits<double>::denorm_min() * 17179869184.0;  // 2^-1040
   const double kSubnormal1060 = std::numeric_limits<double>::denorm_min() * 16384.0;        // 2^-1060
   const double kSmallestNormal = (std::numeric_limits<double>::min)();                      // 2^-1022
   const double kNormal1000 = kSmallestNormal * 4194304.0;                                   // 2^-1000
 
-  // `perm` filters subnormals below `considerZero`, but `perturbCol0` still
-  // treats exact subnormal entries as non-zero.
+  // The merge combines normal and subnormal couplings. Squaring the two
+  // coupling terms directly used to underflow, which later left perturbCol0
+  // without a predecessor and made the decomposition report NumericalIssue.
   matrix.diagonal() << kSubnormal1040, -kSubnormal1060, kSmallestNormal, 0.5, 1.0, kSmallestNormal;
   matrix.diagonal(1) << -kNormal1000, kNormal1000, kSubnormal1040, kSubnormal1060, -8.0;
+  verify_decomposition(matrix);
 
-  BDCSVD<Matrix<double, 6, 6>> svd;
-  svd.setSwitchSize(3);
-  svd.compute(matrix);
-
-  VERIFY(svd.info() == NumericalIssue);
+  // A singular-vector coefficient grows to about 2^570 here. Its squared
+  // norm overflows even though the vector has a finite, well-scaled
+  // normalization.
+  using std::ldexp;
+  matrix.setZero();
+  matrix.diagonal() << 0.0, 0.0, ldexp(1.0, -487), -1.0, 0.0, 0.0;
+  matrix.diagonal(1) << 0.0, ldexp(1.0, -453), -ldexp(1.0, -627), 0.0, 0.0;
+  verify_decomposition(matrix);
 }
 #endif
 
@@ -309,6 +340,6 @@ EIGEN_DECLARE_TEST(bdcsvd) {
   // Bidiagonal SVD hard test cases
   CALL_SUBTEST_51((bdcsvd_bidiagonal_hard_cases<float>()));
   CALL_SUBTEST_52((bdcsvd_bidiagonal_hard_cases<double>()));
-  CALL_SUBTEST_53((bdcsvd_public_missing_predecessor()));
+  CALL_SUBTEST_53((bdcsvd_extreme_scale_regressions()));
   CALL_SUBTEST_54((bdcsvd_fast_math_regression_1588()));
 }

test/bdcsvd_fastmath.cpp

@@ -0,0 +1,29 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include "main.h"
+#include <Eigen/SVD>
+
+EIGEN_DECLARE_TEST(bdcsvd_fastmath) {
+  typedef Matrix<double, 6, 6> Matrix6d;
+  const double kTolerance = 16 * Matrix6d::RowsAtCompileTime * NumTraits<double>::epsilon();
+
+  // A finite singular-vector coefficient grows to about 2^570. Its ordinary
+  // squared norm overflows, and GCC's -ffast-math assumes isfinite() is true.
+  Matrix6d matrix = Matrix6d::Zero();
+  using std::ldexp;
+  matrix.diagonal() << 0.0, 0.0, ldexp(1.0, -487), -1.0, 0.0, 0.0;
+  matrix.diagonal(1) << 0.0, ldexp(1.0, -453), -ldexp(1.0, -627), 0.0, 0.0;
+
+  BDCSVD<Matrix6d, ComputeFullU | ComputeFullV> svd;
+  svd.setSwitchSize(3);
+  svd.compute(matrix);
+
+  VERIFY(svd.info() == Success);
+  const Matrix6d reconstruction = svd.matrixU() * svd.singularValues().asDiagonal() * svd.matrixV().transpose();
+  VERIFY((reconstruction - matrix).stableNorm() <= kTolerance * matrix.stableNorm());
+
+  const Matrix6d identity = Matrix6d::Identity();
+  VERIFY((svd.matrixU().transpose() * svd.matrixU() - identity).stableNorm() <= kTolerance);
+  VERIFY((svd.matrixV().transpose() * svd.matrixV() - identity).stableNorm() <= kTolerance);
+}