Fix numerical issue on SVD with near-identity matrix (#1369)

* 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>
This commit is contained in:
Vollkornaffe 2024-03-28 15:26:11 +01:00 committed by GitHub
parent 749a9fee17
commit c475c4000c
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
2 changed files with 37 additions and 1 deletions

View File

@ -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.

View File

@ -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