Skip to content

Commit

Permalink
Merge pull request #1314 from rasmusgo/fix-svd-near-zero
Browse files Browse the repository at this point in the history
Fix bug in SVD related to values near zero
  • Loading branch information
sebcrozet authored Nov 12, 2023
2 parents f8cd2d4 + 469390f commit 2e99320
Show file tree
Hide file tree
Showing 4 changed files with 62 additions and 10 deletions.
13 changes: 11 additions & 2 deletions src/linalg/bidiagonal.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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..);
Expand Down Expand Up @@ -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..);
Expand Down
6 changes: 3 additions & 3 deletions src/linalg/givens.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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(),
Expand Down Expand Up @@ -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()
Expand Down
39 changes: 34 additions & 5 deletions tests/linalg/bidiagonal.rs
Original file line number Diff line number Diff line change
@@ -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)]
Expand Down Expand Up @@ -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() {
Expand All @@ -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);
}
14 changes: 14 additions & 0 deletions tests/linalg/svd.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

0 comments on commit 2e99320

Please sign in to comment.