fix: SVD loss of accuracy with small non-zero singular values

Author: ABorgna

src/linalg/svd.rs

@@ -159,6 +159,20 @@ where
         let mut diagonal = bi_matrix.diagonal();
         let mut off_diagonal = bi_matrix.off_diagonal();
 
+        // Compute the Frobenius norm of the bidiagonal matrix for relative convergence checks,
+        // so we can use a relative tolerance rather than absolute.
+        // This avoids stability issues when diagonal elements are very small but not exactly zero.
+        let mut norm_sqr = T::RealField::zero();
+        for i in 0..diagonal.len() {
+            norm_sqr += diagonal[i].clone() * diagonal[i].clone();
+        }
+        for i in 0..off_diagonal.len() {
+            norm_sqr += off_diagonal[i].clone() * off_diagonal[i].clone();
+        }
+        let b_norm = norm_sqr.sqrt();
+        // Use relative epsilon: eps * ||B|| for zero checks
+        let eps_rel = eps.clone() * b_norm;
+
         let mut niter = 0;
         let (mut start, mut end) = Self::delimit_subproblem(
             &mut diagonal,
@@ -168,6 +182,7 @@ where
             bi_matrix.is_upper_diagonal(),
             dim - 1,
             eps.clone(),
+            eps_rel.clone(),
         );
 
         while end != start {
@@ -317,6 +332,7 @@ where
                 bi_matrix.is_upper_diagonal(),
                 end,
                 eps.clone(),
+                eps_rel.clone(),
             );
             start = sub.0;
             end = sub.1;
@@ -372,6 +388,7 @@ where
         is_upper_diagonal: bool,
         end: usize,
         eps: T::RealField,
+        eps_rel: T::RealField,
     ) -> (usize, usize) {
         let mut n = end;
 
@@ -383,7 +400,7 @@ where
                     <= eps.clone() * (diagonal[n].clone().norm1() + diagonal[m].clone().norm1())
             {
                 off_diagonal[m] = T::RealField::zero();
-            } else if diagonal[m].clone().norm1() <= eps {
+            } else if diagonal[m].clone().norm1() <= eps_rel {
                 diagonal[m] = T::RealField::zero();
                 Self::cancel_horizontal_off_diagonal_elt(
                     diagonal,
@@ -405,7 +422,7 @@ where
                         m - 1,
                     );
                 }
-            } else if diagonal[n].clone().norm1() <= eps {
+            } else if diagonal[n].clone().norm1() <= eps_rel {
                 diagonal[n] = T::RealField::zero();
                 Self::cancel_vertical_off_diagonal_elt(
                     diagonal,
@@ -435,9 +452,7 @@ where
             {
                 off_diagonal[m] = T::RealField::zero();
                 break;
-            }
-            // TODO: write a test that enters this case.
-            else if diagonal[m].clone().norm1() <= eps {
+            } else if diagonal[m].clone().norm1() <= eps_rel {
                 diagonal[m] = T::RealField::zero();
                 Self::cancel_horizontal_off_diagonal_elt(
                     diagonal,

tests/linalg/svd.rs

@@ -1,4 +1,5 @@
 use crate::utils::is_sorted_descending;
+use approx::{AbsDiffEq, RelativeEq};
 use na::{DMatrix, Matrix6};
 
 #[cfg(feature = "proptest-support")]
@@ -513,3 +514,47 @@ fn svd_regression_issue_1313() {
     let m2 = svd.recompose().unwrap();
     assert_relative_eq!(&m, &m2, epsilon = 1e-5);
 }
+
+#[test]
+// Accuracy bug reported in issue #1172 of nalgebra (https://github.com/dimforge/nalgebra/issues/1172)
+fn svd_regression_issue_1172() {
+    use nalgebra::{Complex, Matrix4};
+    type M4C = Matrix4<Complex<f64>>;
+
+    let m = M4C::new(
+        Complex::new(0.4846888711394364, -0.0000000000000002529226450498658),
+        Complex::new(0.4997655143494952, -0.00000000000000001731471891552503),
+        Complex::new(0.00000000000000001527512211317369, 0.49976551434949495),
+        Complex::new(0.00000000000000009643636194372324, -0.48468887113943676),
+        Complex::new(0.4997655143494954, -0.00000000000000023999992468308935),
+        Complex::new(0.5153111288605629, -0.0000000000000002516108359162637),
+        Complex::new(-0.000000000000000005044961400919933, 0.5153111288605631),
+        Complex::new(0.000000000000000042486770760464153, -0.49976551434949495),
+        Complex::new(-0.000000000000000032383811520876863, -0.49976551434949484),
+        Complex::new(0.000000000000000014710206963221994, -0.5153111288605633),
+        Complex::new(0.5153111288605635, 0.00000000000000016459157448729957),
+        Complex::new(-0.4997655143494946, -0.00000000000000028045339796436483),
+        Complex::new(0.00000000000000007455816688442185, 0.4846888711394367),
+        Complex::new(0.0000000000000000641248522548626, 0.49976551434949484),
+        Complex::new(-0.49976551434949473, -0.00000000000000018487105599880945),
+        Complex::new(0.48468887113943704, 0.00000000000000015953066497724206),
+    );
+    let svd = m.svd(true, true);
+    let sings = svd.singular_values;
+    let u = svd.u.unwrap();
+    let v_t = svd.v_t.unwrap();
+    let sigma = M4C::from_diagonal(&sings.cast::<Complex<f64>>());
+    let m1 = u * sigma * v_t;
+
+    // Should be accurate to machine precision
+    assert_relative_eq!(m, m1, epsilon = 1e-12);
+
+    // Singular values should be sorted and non-negative
+    assert!(sings.iter().all(|&x| x >= 0.0));
+
+    // The largest singular value should be close to 2.0
+    assert!(sings[0].abs_diff_eq(&2.0, 1e-12));
+
+    // The smallest singular value should be small
+    assert!(sings[3] < 1e-12);
+}