Fix bug in SVD related to values near zero

src/linalg/bidiagonal.rs

@@ -8,6 +8,7 @@ use simba::scalar::ComplexField;
 
 use crate::geometry::Reflection;
 use crate::linalg::householder;
+use crate::num::Zero;
 use std::mem::MaybeUninit;
 
 /// The bidiagonalization of a general matrix.
@@ -227,7 +228,11 @@ where
 
         for i in (0..dim - shift).rev() {
             let axis = self.uv.view_range(i + shift.., i);
-            // TODO: sometimes, the axis might have a zero magnitude.
+
+            // Sometimes, the axis might have a zero magnitude.
+            if axis.norm_squared().is_zero() {
+                continue;
+            }
             let refl = Reflection::new(Unit::new_unchecked(axis), T::zero());
 
             let mut res_rows = res.view_range_mut(i + shift.., i..);
@@ -263,7 +268,11 @@ where
             let axis = self.uv.view_range(i, i + shift..);
             let mut axis_packed = axis_packed.rows_range_mut(i + shift..);
             axis_packed.tr_copy_from(&axis);
-            // TODO: sometimes, the axis might have a zero magnitude.
+
+            // Sometimes, the axis might have a zero magnitude.
+            if axis_packed.norm_squared().is_zero() {
+                continue;
+            }
             let refl = Reflection::new(Unit::new_unchecked(axis_packed), T::zero());
 
             let mut res_rows = res.view_range_mut(i.., i + shift..);

src/linalg/givens.rs

@@ -17,7 +17,7 @@ pub struct GivensRotation<T: ComplexField> {
 
 // Matrix = UnitComplex * Matrix
 impl<T: ComplexField> GivensRotation<T> {
-    /// The Givents rotation that does nothing.
+    /// The Givens rotation that does nothing.
     pub fn identity() -> Self {
         Self {
             c: T::RealField::one(),
@@ -88,13 +88,13 @@ impl<T: ComplexField> GivensRotation<T> {
         }
     }
 
-    /// The cos part of this roration.
+    /// The cos part of this rotation.
     #[must_use]
     pub fn c(&self) -> T::RealField {
         self.c.clone()
     }
 
-    /// The sin part of this roration.
+    /// The sin part of this rotation.
     #[must_use]
     pub fn s(&self) -> T {
         self.s.clone()

tests/linalg/bidiagonal.rs

@@ -1,6 +1,6 @@
-#![cfg(feature = "proptest-support")]
-
-macro_rules! gen_tests(
+#[cfg(feature = "proptest-support")]
+mod proptest_tests {
+    macro_rules! gen_tests(
     ($module: ident, $scalar: expr) => {
         mod $module {
             #[allow(unused_imports)]
@@ -54,8 +54,9 @@ macro_rules! gen_tests(
     }
 );
 
-gen_tests!(complex, complex_f64());
-gen_tests!(f64, PROPTEST_F64);
+    gen_tests!(complex, complex_f64());
+    gen_tests!(f64, PROPTEST_F64);
+}
 
 #[test]
 fn bidiagonal_identity() {
@@ -74,3 +75,31 @@ fn bidiagonal_identity() {
     let (u, d, v_t) = bidiagonal.unpack();
     assert_eq!(m, &u * d * &v_t);
 }
+
+#[test]
+fn bidiagonal_regression_issue_1313() {
+    let s = 6.123234e-16_f32;
+    let mut m = nalgebra::dmatrix![
+        10.0,   0.0, 0.0,  0.0, -10.0, 0.0, 0.0, 0.0;
+        s,     10.0, 0.0, 10.0,     s, 0.0, 0.0, 0.0;
+        20.0, -20.0, 0.0, 20.0,  20.0, 0.0, 0.0, 0.0;
+    ];
+    m.unscale_mut(m.camax());
+    let bidiagonal = m.clone().bidiagonalize();
+    let (u, d, v_t) = bidiagonal.unpack();
+    let m2 = &u * d * &v_t;
+    assert_relative_eq!(m, m2, epsilon = 1e-6);
+}
+
+#[test]
+fn bidiagonal_regression_issue_1313_minimal() {
+    let s = 6.123234e-17_f32;
+    let m = nalgebra::dmatrix![
+        1.0,   0.0, -1.0;
+        s,     1.0,     s;
+    ];
+    let bidiagonal = m.clone().bidiagonalize();
+    let (u, d, v_t) = bidiagonal.unpack();
+    let m2 = &u * &d * &v_t;
+    assert_relative_eq!(m, m2, epsilon = 1e-6);
+}

tests/linalg/svd.rs

@@ -499,3 +499,17 @@ fn svd_regression_issue_1072() {
         epsilon = 1e-9
     );
 }
+
+#[test]
+// Exercises bug reported in issue #1313 of nalgebra (https://github.com/dimforge/nalgebra/issues/1313)
+fn svd_regression_issue_1313() {
+    let s = 6.123234e-16_f32;
+    let m = nalgebra::dmatrix![
+        10.0,   0.0, 0.0,  0.0, -10.0, 0.0, 0.0, 0.0;
+           s,  10.0, 0.0, 10.0,     s, 0.0, 0.0, 0.0;
+        20.0, -20.0, 0.0, 20.0,  20.0, 0.0, 0.0, 0.0;
+    ];
+    let svd = m.clone().svd(true, true);
+    let m2 = svd.recompose().unwrap();
+    assert_relative_eq!(&m, &m2, epsilon = 1e-5);
+}