Add bidiagonal SVD API to BDCSVD and remove dead debug code

libeigen/eigen!2238 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
2026-04-10 11:34:33 +08:00 · 2026-03-28 20:38:31 -07:00
parent ba9871e46b
commit 624ab58e8d
3 changed files with 354 additions and 286 deletions
--- a/test/bdcsvd.cpp
+++ b/test/bdcsvd.cpp
@@ -85,6 +85,212 @@ void bdcsvd_check_convergence(const MatrixType& input) {
  VERIFY_IS_APPROX(input, D);
 }

+// Verify SVD of bidiagonal matrix given as diagonal + superdiagonal vectors.
+template <typename RealScalar>
+void verify_bidiagonal_svd(const Matrix<RealScalar, Dynamic, 1>& diag,
+                           const Matrix<RealScalar, Dynamic, 1>& superdiag) {
+  typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixXr;
+  typedef Matrix<RealScalar, Dynamic, 1> VectorXr;
+  const Index n = diag.size();
+
+  BDCSVD<MatrixXr, ComputeFullU | ComputeFullV> bdcsvd(diag, superdiag);
+  VERIFY(bdcsvd.info() == Success);
+
+  const VectorXr& sv = bdcsvd.singularValues();
+
+  // Singular values must be non-negative.
+  for (Index i = 0; i < sv.size(); ++i) {
+    VERIFY(sv(i) >= RealScalar(0));
+  }
+
+  // Singular values must be sorted descending.
+  for (Index i = 1; i < sv.size(); ++i) {
+    VERIFY(sv(i - 1) >= sv(i));
+  }
+
+  // Orthogonality of U and V.
+  VERIFY_IS_APPROX(bdcsvd.matrixU().transpose() * bdcsvd.matrixU(), MatrixXr::Identity(n, n));
+  VERIFY_IS_APPROX(bdcsvd.matrixV().transpose() * bdcsvd.matrixV(), MatrixXr::Identity(n, n));
+
+  // Reconstruction: U * S * V^T should equal the original bidiagonal.
+  MatrixXr B = MatrixXr::Zero(n, n);
+  B.diagonal() = diag;
+  if (n > 1) B.diagonal(1) = superdiag;
+  MatrixXr recon = bdcsvd.matrixU() * sv.asDiagonal() * bdcsvd.matrixV().transpose();
+  VERIFY_IS_APPROX(recon, B);
+
+  // Cross-validate singular values against JacobiSVD.
+  JacobiSVD<MatrixXr> jacobi(B);
+  VERIFY_IS_APPROX(sv, jacobi.singularValues());
+}
+
+// Verify that bidiagonal API and matrix API produce matching singular values.
+template <typename RealScalar>
+void verify_bidiagonal_vs_matrix_svd(const Matrix<RealScalar, Dynamic, 1>& diag,
+                                     const Matrix<RealScalar, Dynamic, 1>& superdiag) {
+  typedef Matrix<RealScalar, Dynamic, Dynamic> MatrixXr;
+  const Index n = diag.size();
+
+  // Build dense bidiagonal matrix.
+  MatrixXr B = MatrixXr::Zero(n, n);
+  B.diagonal() = diag;
+  if (n > 1) B.diagonal(1) = superdiag;
+
+  BDCSVD<MatrixXr> bidiag_svd(diag, superdiag);
+  BDCSVD<MatrixXr> matrix_svd(B);
+
+  VERIFY(bidiag_svd.info() == Success);
+  VERIFY(matrix_svd.info() == Success);
+  VERIFY_IS_APPROX(bidiag_svd.singularValues(), matrix_svd.singularValues());
+}
+
+template <typename RealScalar>
+void bdcsvd_bidiagonal_hard_cases() {
+  using std::abs;
+  using std::cos;
+  using std::pow;
+  using std::sin;
+  typedef Matrix<RealScalar, Dynamic, 1> VectorXr;
+
+  Eigen::internal::set_is_malloc_allowed(true);
+
+  const RealScalar eps = NumTraits<RealScalar>::epsilon();
+
+  // Test sizes: cover n=1, very small, below/above algoSwap (16), and larger.
+  const int sizes[] = {1, 2, 3, 5, 10, 16, 20, 50, 100};
+  const int numSizes = sizeof(sizes) / sizeof(sizes[0]);
+
+  for (int si = 0; si < numSizes; ++si) {
+    const Index n = sizes[si];
+    VectorXr diag(n), superdiag(n > 1 ? n - 1 : 0);
+
+    // 1. Identity: d=[1,...,1], e=[0,...,0]
+    diag.setOnes();
+    superdiag.setZero();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+    verify_bidiagonal_vs_matrix_svd<RealScalar>(diag, superdiag);
+
+    // 2. Zero: d=[0,...,0], e=[0,...,0]
+    diag.setZero();
+    superdiag.setZero();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 3. Scalar (only meaningful for n=1, but runs for all)
+    if (n == 1) {
+      diag(0) = RealScalar(3.14);
+      verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+    }
+
+    // 4. Golub-Kahan: d=[1,...,1], e=[1,...,1]
+    diag.setOnes();
+    if (n > 1) superdiag.setOnes();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 5. Kahan matrix: d_i = s^(i-1), e_i = -c*s^(i-1)
+    // Clamp exponents so condition number stays bounded by 1/eps.
+    {
+      const RealScalar theta = RealScalar(0.3);
+      const RealScalar s = sin(theta);
+      const RealScalar c = cos(theta);
+      using std::log;
+      const RealScalar maxPower = -log(eps) / (-log(s));
+      for (Index i = 0; i < n; ++i) diag(i) = pow(s, numext::mini(RealScalar(i), maxPower));
+      for (Index i = 0; i < n - 1; ++i) superdiag(i) = -c * pow(s, numext::mini(RealScalar(i), maxPower));
+      verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+    }
+
+    // 6. Geometric decay diagonal: d_i = 0.5^i, e=[0,...,0]
+    // Clamp so condition number stays bounded by 1/eps.
+    {
+      using std::log;
+      const RealScalar base = RealScalar(0.5);
+      const RealScalar maxPower = -log(eps) / (-log(base));
+      for (Index i = 0; i < n; ++i) diag(i) = pow(base, numext::mini(RealScalar(i), maxPower));
+      superdiag.setZero();
+      verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+    }
+
+    // 7. Geometric decay superdiagonal: d=[1,...,1], e_i = 0.5^i
+    diag.setOnes();
+    for (Index i = 0; i < n - 1; ++i) superdiag(i) = pow(RealScalar(0.5), RealScalar(i));
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 8. Clustered at 1: d_i = 1 + i*eps, e=[0,...,0]
+    for (Index i = 0; i < n; ++i) diag(i) = RealScalar(1) + RealScalar(i) * eps;
+    superdiag.setZero();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 9. Two clusters: half ≈ 1, half ≈ eps
+    for (Index i = 0; i < n; ++i) diag(i) = (i < n / 2) ? RealScalar(1) : eps;
+    superdiag.setZero();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 10. Single tiny singular value: d=[1,...,1,eps], e=[eps^2,...]
+    diag.setOnes();
+    diag(n - 1) = eps;
+    for (Index i = 0; i < n - 1; ++i) superdiag(i) = eps * eps;
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 11. Graded: d_i = 10^(-i), e_i = 10^(-i)
+    for (Index i = 0; i < n; ++i) diag(i) = pow(RealScalar(10), -RealScalar(i));
+    for (Index i = 0; i < n - 1; ++i) superdiag(i) = pow(RealScalar(10), -RealScalar(i));
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 12. Nearly diagonal: random diag, eps * random superdiag
+    diag = VectorXr::Random(n).cwiseAbs() + VectorXr::Constant(n, RealScalar(0.1));
+    for (Index i = 0; i < n - 1; ++i) superdiag(i) = eps * (RealScalar(0.5) + abs(internal::random<RealScalar>()));
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 13. All equal: d=[c,...,c], e=[c,...,c]
+    diag.setConstant(RealScalar(2.5));
+    if (n > 1) superdiag.setConstant(RealScalar(2.5));
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 14. Wilkinson: d_i = |n/2 - i|, e=[1,...,1]
+    for (Index i = 0; i < n; ++i) diag(i) = abs(RealScalar(n / 2) - RealScalar(i));
+    if (n > 1) superdiag.setOnes();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 15. Overflow/underflow: alternating big/tiny diagonal, tiny/big superdiagonal
+    {
+      const RealScalar big = (std::numeric_limits<RealScalar>::max)() / RealScalar(1000);
+      const RealScalar tiny = (std::numeric_limits<RealScalar>::min)() * RealScalar(1000);
+      for (Index i = 0; i < n; ++i) diag(i) = (i % 2 == 0) ? big : tiny;
+      for (Index i = 0; i < n - 1; ++i) superdiag(i) = (i % 2 == 0) ? tiny : big;
+      verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+    }
+
+    // 16. Prescribed condition number: d_i = kappa^(-i/(n-1)), e_i = eps * random
+    if (n > 1) {
+      const RealScalar kappa = RealScalar(1) / eps;
+      for (Index i = 0; i < n; ++i) diag(i) = pow(kappa, -RealScalar(i) / RealScalar(n - 1));
+      for (Index i = 0; i < n - 1; ++i) superdiag(i) = eps * abs(internal::random<RealScalar>());
+      verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+    }
+
+    // 17. Rank-deficient: d=[1,..,0,..,0,..,1], e=[0,...,0]
+    for (Index i = 0; i < n; ++i) diag(i) = (i < n / 3 || i >= 2 * n / 3) ? RealScalar(1) : RealScalar(0);
+    superdiag.setZero();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 18. Arrowhead stress: d_i = linspace(1, n), e_i = 1/(i+1)
+    for (Index i = 0; i < n; ++i) diag(i) = RealScalar(1) + RealScalar(i);
+    for (Index i = 0; i < n - 1; ++i) superdiag(i) = RealScalar(1) / RealScalar(i + 1);
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 19. Repeated singular values: d=[1,2,3,1,2,3,...], e=[0,...,0]
+    for (Index i = 0; i < n; ++i) diag(i) = RealScalar((i % 3) + 1);
+    superdiag.setZero();
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+
+    // 20. Glued identity: d=[1,...,1], e=0 except e[n/2-1]=eps
+    diag.setOnes();
+    superdiag.setZero();
+    if (n > 2) superdiag(n / 2 - 1) = eps;
+    verify_bidiagonal_svd<RealScalar>(diag, superdiag);
+  }
+}
+
 EIGEN_DECLARE_TEST(bdcsvd) {
  CALL_SUBTEST_1((bdcsvd_verify_assert<Matrix3f>()));
  CALL_SUBTEST_2((bdcsvd_verify_assert<Matrix4d>()));
@@ -174,4 +380,8 @@ EIGEN_DECLARE_TEST(bdcsvd) {

  // Convergence for large constant matrix (https://gitlab.com/libeigen/eigen/-/issues/2491)
  CALL_SUBTEST_49(bdcsvd_check_convergence<MatrixXf>(MatrixXf::Constant(500, 500, 1)));
+
+  // Bidiagonal SVD hard test cases
+  CALL_SUBTEST_50((bdcsvd_bidiagonal_hard_cases<float>()));
+  CALL_SUBTEST_51((bdcsvd_bidiagonal_hard_cases<double>()));
 }