From 8f800d905e383e0b17a7a04c576ee0874589a148 Mon Sep 17 00:00:00 2001 From: Vollkornaffe Date: Fri, 8 Mar 2024 17:17:08 +0100 Subject: [PATCH 1/4] 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. --- src/linalg/householder.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs index 08212c675..0b3f1bab8 100644 --- a/src/linalg/householder.rs +++ b/src/linalg/householder.rs @@ -34,6 +34,7 @@ pub fn reflection_axis_mut>( if !factor.is_zero() { column.unscale_mut(factor.sqrt()); + let _ = column.normalize_mut(); (-signed_norm, true) } else { // TODO: not sure why we don't have a - sign here. From c1602994eb0d24f447dd2371df87be82d8753928 Mon Sep 17 00:00:00 2001 From: Vollkornaffe Date: Fri, 8 Mar 2024 17:51:17 +0100 Subject: [PATCH 2/4] chore: Add test that fails before fix --- tests/linalg/eigen.rs | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index a5dcf8357..3fdd8e53b 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() { ); } +// Test for #1368 +#[test] +fn very_small_deviation_from_identity() { + let m = Matrix3::::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 From 4dad6b31e6e626da636a243492807a411a7be6e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Thu, 28 Mar 2024 15:07:05 +0100 Subject: [PATCH 3/4] chore: add comment providing details on the householder fix. --- src/linalg/householder.rs | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs index 0b3f1bab8..dbae93f96 100644 --- a/src/linalg/householder.rs +++ b/src/linalg/householder.rs @@ -34,7 +34,17 @@ pub fn reflection_axis_mut>( 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. From 5e70d049074ee785d3b06fdd8a45e302189a3483 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Thu, 28 Mar 2024 15:18:33 +0100 Subject: [PATCH 4/4] =?UTF-8?q?chore:=E2=80=AFrename=20regression=20test?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- tests/linalg/eigen.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index 3fdd8e53b..7a944c444 100644 --- a/tests/linalg/eigen.rs +++ b/tests/linalg/eigen.rs @@ -116,9 +116,9 @@ fn symmetric_eigen_singular_24x24() { ); } -// Test for #1368 +// Regression test for #1368 #[test] -fn very_small_deviation_from_identity() { +fn very_small_deviation_from_identity_issue_1368() { let m = Matrix3::::new( 1.0, 3.1575704e-23,