Merge pull request #1314 from rasmusgo/fix-svd-near-zero

Fix bug in SVD related to values near zero
This commit is contained in:
Sébastien Crozet 2023-11-12 23:18:59 +01:00 committed by GitHub
commit 2e99320d01
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
4 changed files with 62 additions and 10 deletions

View File

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

View File

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

View File

@ -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);
}

View File

@ -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);
}