From b6e094c82ff1a03a68ed8d1b742dd51f065982c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Br=C3=B6nneg=C3=A5rd?= <1162652+rasmusgo@users.noreply.github.com> Date: Thu, 9 Nov 2023 00:48:13 +0100 Subject: [PATCH 1/4] Fix spelling in givens.rs --- src/linalg/givens.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs index 30e2bed38..dbd254391 100644 --- a/src/linalg/givens.rs +++ b/src/linalg/givens.rs @@ -17,7 +17,7 @@ pub struct GivensRotation { // Matrix = UnitComplex * Matrix impl GivensRotation { - /// 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 GivensRotation { } } - /// 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() From 03c24fb369f5f3873a9a93d6e75661b39199e132 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Br=C3=B6nneg=C3=A5rd?= <1162652+rasmusgo@users.noreply.github.com> Date: Thu, 9 Nov 2023 00:53:42 +0100 Subject: [PATCH 2/4] Add unit tests for issue 1313 --- tests/linalg/bidiagonal.rs | 39 +++++++++++++++++++++++++++++++++----- tests/linalg/svd.rs | 14 ++++++++++++++ 2 files changed, 48 insertions(+), 5 deletions(-) diff --git a/tests/linalg/bidiagonal.rs b/tests/linalg/bidiagonal.rs index aaee393fb..800967335 100644 --- a/tests/linalg/bidiagonal.rs +++ b/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); +} diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 900901ad3..d8b23d024 100644 --- a/tests/linalg/svd.rs +++ b/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); +} From 7ea9ecee081fc80f2c0c3a23f0985abddfdd5d05 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rasmus=20Br=C3=B6nneg=C3=A5rd?= <1162652+rasmusgo@users.noreply.github.com> Date: Thu, 9 Nov 2023 01:20:44 +0100 Subject: [PATCH 3/4] Test for axes with zero magnitude --- src/linalg/bidiagonal.rs | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs index 3036e9f15..fd32e744e 100644 --- a/src/linalg/bidiagonal.rs +++ b/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,10 @@ 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.magnitude().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 +267,10 @@ 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.magnitude().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..); From 469390f4b922e8c0d455ead00dfd99e8ea610671 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sun, 12 Nov 2023 23:12:52 +0100 Subject: [PATCH 4/4] Check norm_squared instead of mangnitude. --- src/linalg/bidiagonal.rs | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs index fd32e744e..4faa9c5a0 100644 --- a/src/linalg/bidiagonal.rs +++ b/src/linalg/bidiagonal.rs @@ -228,8 +228,9 @@ where for i in (0..dim - shift).rev() { let axis = self.uv.view_range(i + shift.., i); + // Sometimes, the axis might have a zero magnitude. - if axis.magnitude().is_zero() { + if axis.norm_squared().is_zero() { continue; } let refl = Reflection::new(Unit::new_unchecked(axis), T::zero()); @@ -267,8 +268,9 @@ 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); + // Sometimes, the axis might have a zero magnitude. - if axis_packed.magnitude().is_zero() { + if axis_packed.norm_squared().is_zero() { continue; } let refl = Reflection::new(Unit::new_unchecked(axis_packed), T::zero());