From 325618ba22660b5ee22f308b4f06ac5f0196a9ae Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 9 Mar 2022 02:13:12 -0500 Subject: [PATCH 1/3] Fix SVD instability bug --- src/linalg/svd.rs | 22 +++++++--------------- tests/linalg/svd.rs | 24 ++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 15 deletions(-) diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 3f945a65..42a6abb3 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -2,7 +2,6 @@ use serde::{Deserialize, Serialize}; use std::any::TypeId; -use approx::AbsDiffEq; use num::{One, Zero}; use crate::allocator::Allocator; @@ -94,14 +93,7 @@ where /// The singular values are not guaranteed to be sorted in any particular order. /// If a descending order is required, consider using `new` instead. pub fn new_unordered(matrix: OMatrix, compute_u: bool, compute_v: bool) -> Self { - Self::try_new_unordered( - matrix, - compute_u, - compute_v, - T::RealField::default_epsilon(), - 0, - ) - .unwrap() + Self::try_new_unordered(matrix, compute_u, compute_v, crate::convert(1e-15), 0).unwrap() } /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. @@ -888,13 +880,13 @@ fn compute_2x2_uptrig_svd( v_t = Some(csv.clone()); } - if compute_u { - let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1.clone(); - let su = (m22 * csv.s()) / v1.clone(); - let (csu, sgn_u) = GivensRotation::new(cu, su); + let cu = (m11.scale(csv.c()) + m12 * csv.s()) / v1.clone(); + let su = (m22 * csv.s()) / v1.clone(); + let (csu, sgn_u) = GivensRotation::new(cu, su); + v1 *= sgn_u.clone(); + v2 *= sgn_u; - v1 *= sgn_u.clone(); - v2 *= sgn_u; + if compute_u { u = Some(csu); } } diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index deb3d38d..7eabe3d0 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -460,3 +460,27 @@ fn svd_sorted() { epsilon = 1.0e-5 ); } + +#[test] +// Exercises bug reported in issue #983 of nalgebra +fn svd_consistent() { + let m = nalgebra::dmatrix![ + 10.74785316637712f64, -5.994983325167452, -6.064492921857296; + -4.149751381521569, 20.654504205822462, -4.470436210703133; + -22.772715014220207, -1.4554372570788008, 18.108113992170573 + ] + .transpose(); + let svd1 = m.clone().svd(true, true); + let svd2 = m.clone().svd(false, true); + let svd3 = m.clone().svd(true, false); + let svd4 = m.svd(false, false); + + assert_relative_eq!(svd1.singular_values, svd2.singular_values, epsilon = 1e-5); + assert_relative_eq!(svd1.singular_values, svd3.singular_values, epsilon = 1e-5); + assert_relative_eq!(svd1.singular_values, svd4.singular_values, epsilon = 1e-5); + assert_relative_eq!( + svd1.singular_values, + nalgebra::dvector![3.16188022e+01, 2.23811978e+01, 0.], + epsilon = 1e-5 + ); +} From 1acd48f6f1b684317e7cd893d973ba67ab4f7290 Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 9 Mar 2022 21:04:43 -0500 Subject: [PATCH 2/3] Address review comments --- src/linalg/svd.rs | 10 +++++++++- tests/linalg/svd.rs | 12 ++++++------ 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 42a6abb3..06bae4a3 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -2,6 +2,7 @@ use serde::{Deserialize, Serialize}; use std::any::TypeId; +use approx::AbsDiffEq; use num::{One, Zero}; use crate::allocator::Allocator; @@ -93,7 +94,14 @@ where /// The singular values are not guaranteed to be sorted in any particular order. /// If a descending order is required, consider using `new` instead. pub fn new_unordered(matrix: OMatrix, compute_u: bool, compute_v: bool) -> Self { - Self::try_new_unordered(matrix, compute_u, compute_v, crate::convert(1e-15), 0).unwrap() + Self::try_new_unordered( + matrix, + compute_u, + compute_v, + T::RealField::default_epsilon() * crate::convert(5.0), + 0, + ) + .unwrap() } /// Attempts to compute the Singular Value Decomposition of `matrix` using implicit shift. diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 7eabe3d0..8e83df81 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -462,8 +462,8 @@ fn svd_sorted() { } #[test] -// Exercises bug reported in issue #983 of nalgebra -fn svd_consistent() { +// Exercises bug reported in issue #983 of nalgebra (https://github.com/dimforge/nalgebra/issues/983) +fn svd_regression_issue_983() { let m = nalgebra::dmatrix![ 10.74785316637712f64, -5.994983325167452, -6.064492921857296; -4.149751381521569, 20.654504205822462, -4.470436210703133; @@ -475,12 +475,12 @@ fn svd_consistent() { let svd3 = m.clone().svd(true, false); let svd4 = m.svd(false, false); - assert_relative_eq!(svd1.singular_values, svd2.singular_values, epsilon = 1e-5); - assert_relative_eq!(svd1.singular_values, svd3.singular_values, epsilon = 1e-5); - assert_relative_eq!(svd1.singular_values, svd4.singular_values, epsilon = 1e-5); + assert_relative_eq!(svd1.singular_values, svd2.singular_values, epsilon = 1e-9); + assert_relative_eq!(svd1.singular_values, svd3.singular_values, epsilon = 1e-9); + assert_relative_eq!(svd1.singular_values, svd4.singular_values, epsilon = 1e-9); assert_relative_eq!( svd1.singular_values, nalgebra::dvector![3.16188022e+01, 2.23811978e+01, 0.], - epsilon = 1e-5 + epsilon = 1e-6 ); } From a27d121a7a1783c822ff47c46aee0475bda15852 Mon Sep 17 00:00:00 2001 From: YuhanLiin Date: Wed, 9 Mar 2022 21:10:45 -0500 Subject: [PATCH 3/3] Add regression test for #1072 --- tests/linalg/svd.rs | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 8e83df81..900901ad 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -484,3 +484,18 @@ fn svd_regression_issue_983() { epsilon = 1e-6 ); } + +#[test] +// Exercises bug reported in issue #1072 of nalgebra (https://github.com/dimforge/nalgebra/issues/1072) +fn svd_regression_issue_1072() { + let x = nalgebra::dmatrix![-6.206610118536945f64, -3.67612186839874; -1.2755730783423473, 6.047238193479124]; + let mut x_svd = x.svd(true, true); + x_svd.singular_values = nalgebra::dvector![1.0, 0.0]; + let y = x_svd.recompose().unwrap(); + let y_svd = y.svd(true, true); + assert_relative_eq!( + y_svd.singular_values, + nalgebra::dvector![1.0, 0.0], + epsilon = 1e-9 + ); +}