From c475c4000c12f585f65731b32ef4bf1ecdcf90ad Mon Sep 17 00:00:00 2001
From: Vollkornaffe <larshelge.scheel@gmail.com>
Date: Thu, 28 Mar 2024 15:26:11 +0100
Subject: [PATCH] Fix numerical issue on SVD with near-identity matrix (#1369)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* fix: Normalize the column once more

The column may not be normalized if the `factor` is on a scale of 1e-40.
Possibly, f32 just runs out of precision.

There is likely a better solution to the problem.

* chore: Add test that fails before fix

* chore: add comment providing details on the householder fix.

* chore: rename regression test

---------

Co-authored-by: Sébastien Crozet <sebcrozet@dimforge.com>
---
 src/linalg/householder.rs | 11 +++++++++++
 tests/linalg/eigen.rs     | 27 ++++++++++++++++++++++++++-
 2 files changed, 37 insertions(+), 1 deletion(-)

diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs
index 08212c67..dbae93f9 100644
--- a/src/linalg/householder.rs
+++ b/src/linalg/householder.rs
@@ -34,6 +34,17 @@ pub fn reflection_axis_mut<T: ComplexField, D: Dim, S: StorageMut<T, D>>(
 
     if !factor.is_zero() {
         column.unscale_mut(factor.sqrt());
+
+        // Normalize again, making sure the vector is unit-sized.
+        // If `factor` had a very small value, the first normalization
+        // (dividing by `factor.sqrt()`) might end up with a slightly
+        // non-unit vector (especially when using 32-bits float).
+        // Decompositions strongly rely on that unit-vector property,
+        // so we run a second normalization (that is much more numerically
+        // stable since the norm is close to 1) to ensure it has a unit
+        // size.
+        let _ = column.normalize_mut();
+
         (-signed_norm, true)
     } else {
         // TODO: not sure why we don't have a - sign here.
diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs
index a5dcf835..7a944c44 100644
--- a/tests/linalg/eigen.rs
+++ b/tests/linalg/eigen.rs
@@ -1,4 +1,4 @@
-use na::DMatrix;
+use na::{DMatrix, Matrix3};
 
 #[cfg(feature = "proptest-support")]
 mod proptest_tests {
@@ -116,6 +116,31 @@ fn symmetric_eigen_singular_24x24() {
     );
 }
 
+// Regression test for #1368
+#[test]
+fn very_small_deviation_from_identity_issue_1368() {
+    let m = Matrix3::<f32>::new(
+        1.0,
+        3.1575704e-23,
+        8.1146196e-23,
+        3.1575704e-23,
+        1.0,
+        1.7471054e-22,
+        8.1146196e-23,
+        1.7471054e-22,
+        1.0,
+    );
+
+    for v in m
+        .try_symmetric_eigen(f32::EPSILON, 0)
+        .unwrap()
+        .eigenvalues
+        .into_iter()
+    {
+        assert_relative_eq!(*v, 1.);
+    }
+}
+
 //  #[cfg(feature = "arbitrary")]
 //  quickcheck! {
 // TODO: full eigendecomposition is not implemented yet because of its complexity when some