From a2aae083ce2ebd5c36a759d11cf95bc2264a5f47 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Fri, 3 Jan 2025 22:45:28 +0100 Subject: [PATCH 1/8] add reproducible_matrix function for const sized matrices --- benches/lib.rs | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/benches/lib.rs b/benches/lib.rs index 48633317d..429ba0e3b 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -6,8 +6,9 @@ extern crate rand_package as rand; #[macro_use] extern crate criterion; -use na::DMatrix; +use na::{DMatrix, Matrix, Scalar}; use rand::Rng; +use rand_distr::{Distribution, Standard}; use rand_isaac::IsaacRng; pub mod core; @@ -20,6 +21,19 @@ fn reproducible_dmatrix(nrows: usize, ncols: usize) -> DMatrix { DMatrix::::from_fn(nrows, ncols, |_, _| rng.gen()) } +fn reproducible_matrix( +) -> Matrix, na::Const, na::ArrayStorage> +where + [[T; R]; C]: Default, + Standard: Distribution, +{ + use rand::SeedableRng; + let mut rng = IsaacRng::seed_from_u64(0); + let mut m: Matrix, na::Const, na::ArrayStorage> = Matrix::default(); + m.iter_mut().for_each(|x| *x = rng.gen()); + return m; +} + criterion_main!( core::matrix, core::vector, @@ -34,4 +48,5 @@ criterion_main!( linalg::solve, linalg::svd, linalg::symmetric_eigen, + linalg::matrix_inversion ); From 0c0d292366ec970387f5e40ec83995faf8607622 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Fri, 3 Jan 2025 22:45:43 +0100 Subject: [PATCH 2/8] add matrix inversion benchmark for 4x4 case --- benches/linalg/inverse.rs | 19 +++++++++++++++++++ benches/linalg/mod.rs | 3 ++- 2 files changed, 21 insertions(+), 1 deletion(-) create mode 100644 benches/linalg/inverse.rs diff --git a/benches/linalg/inverse.rs b/benches/linalg/inverse.rs new file mode 100644 index 000000000..2dc31f419 --- /dev/null +++ b/benches/linalg/inverse.rs @@ -0,0 +1,19 @@ +use na::LU; + +use crate::reproducible_matrix; + +fn inverse_4(bh: &mut criterion::Criterion) { + let m = reproducible_matrix::(); + bh.bench_function("4x4_matrix_inverse_det", move |bh| { + bh.iter(|| std::hint::black_box(m.try_inverse().unwrap())) + }); +} + +fn inverse_lu(bh: &mut criterion::Criterion) { + let m = reproducible_matrix::(); + bh.bench_function("4x4_matrix_inverse_lu", move |bh| { + bh.iter(|| std::hint::black_box(LU::new(m).try_inverse().unwrap())) + }); +} + +criterion_group!(matrix_inversion, inverse_4, inverse_lu); diff --git a/benches/linalg/mod.rs b/benches/linalg/mod.rs index cbd4df4a7..528bb3ebd 100644 --- a/benches/linalg/mod.rs +++ b/benches/linalg/mod.rs @@ -2,17 +2,18 @@ pub use self::bidiagonal::bidiagonal; pub use self::cholesky::cholesky; pub use self::full_piv_lu::full_piv_lu; pub use self::hessenberg::hessenberg; +pub use self::inverse::matrix_inversion; pub use self::lu::lu; pub use self::qr::qr; pub use self::schur::schur; pub use self::solve::solve; pub use self::svd::svd; pub use self::symmetric_eigen::symmetric_eigen; - mod bidiagonal; mod cholesky; mod full_piv_lu; mod hessenberg; +mod inverse; mod lu; mod qr; mod schur; From a64963021ceba8278264d9d42e1e5f2b6aa68e30 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Fri, 3 Jan 2025 22:54:40 +0100 Subject: [PATCH 3/8] add numerical check to SquareMatrix::do_inverse4 method --- src/linalg/inverse.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/linalg/inverse.rs b/src/linalg/inverse.rs index b11eae13a..4d96514fd 100644 --- a/src/linalg/inverse.rs +++ b/src/linalg/inverse.rs @@ -278,6 +278,9 @@ where for j in 0..4 { for i in 0..4 { out[(i, j)] *= inv_det.clone(); + if !out[(i, j)].is_finite() { + return false; + } } } true From f696cd529238ba311c57e1a8f2643f01490128b6 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Fri, 3 Jan 2025 23:10:56 +0100 Subject: [PATCH 4/8] add test for uncatched non numeric values for near singular matrices --- tests/linalg/inverse.rs | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tests/linalg/inverse.rs b/tests/linalg/inverse.rs index ab641a791..83659c640 100644 --- a/tests/linalg/inverse.rs +++ b/tests/linalg/inverse.rs @@ -140,3 +140,25 @@ fn matrix5_try_inverse_scaled_identity() { assert_relative_eq!(a_inv, expected_inverse); } + +#[test] +#[rustfmt::skip] +fn matrix4_try_inverse_singular(){ + let m1 = Matrix4::new( + 1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 0.0 ); + assert!(m1.try_inverse().is_none()); +} + +#[test] +#[rustfmt::skip] +fn matrix4_try_inverse_near_singular(){ + let m1 = Matrix4::new( + 1.0, 1e1, 1e1, 1e1, + 0.0, 1.0, 1e1, 1e1, + 0.0, 0.0, 1.0, 1e1, + 0.0, 0.0, 0.0, f64::MIN_POSITIVE ); + assert!(m1.try_inverse().is_none()); +} From 0661f2b0250cd0f7953b344731b9fbaac464f5e4 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Fri, 3 Jan 2025 23:17:54 +0100 Subject: [PATCH 5/8] fix format (whitespace) --- benches/linalg/mod.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/benches/linalg/mod.rs b/benches/linalg/mod.rs index 528bb3ebd..4477654d2 100644 --- a/benches/linalg/mod.rs +++ b/benches/linalg/mod.rs @@ -9,6 +9,7 @@ pub use self::schur::schur; pub use self::solve::solve; pub use self::svd::svd; pub use self::symmetric_eigen::symmetric_eigen; + mod bidiagonal; mod cholesky; mod full_piv_lu; From 2b155017c390128310c5ca22a72bf970de96e114 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Sat, 4 Jan 2025 00:21:40 +0100 Subject: [PATCH 6/8] Update documentation for the desired behaviour, ensure, is_finite checks are done for all matrix sizes --- src/linalg/inverse.rs | 12 +++++++---- tests/linalg/inverse.rs | 47 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 54 insertions(+), 5 deletions(-) diff --git a/src/linalg/inverse.rs b/src/linalg/inverse.rs index 4d96514fd..451ba8f16 100644 --- a/src/linalg/inverse.rs +++ b/src/linalg/inverse.rs @@ -10,6 +10,11 @@ use crate::linalg::lu; impl> SquareMatrix { /// Attempts to invert this square matrix. /// + /// # Returns + /// + /// - `Some`: The inverted matrix, if the inversion succeeds and all entries are finite. + /// - `None`: If the inversion fails or produces `Inf` or `NaN` values. + /// /// # Panics /// /// Panics if `self` isn’t a square matrix. @@ -20,7 +25,7 @@ impl> SquareMatrix { DefaultAllocator: Allocator, { let mut me = self.into_owned(); - if me.try_inverse_mut() { + if me.try_inverse_mut() && me.iter().all(|x| x.is_finite()) { Some(me) } else { None @@ -31,6 +36,8 @@ impl> SquareMatrix { impl> SquareMatrix { /// Attempts to invert this square matrix in-place. Returns `false` and leaves `self` untouched if /// inversion fails. + /// + /// Calling this function on near singular matrices may lead to`Inf` or `NaN` values! /// /// # Panics /// @@ -278,9 +285,6 @@ where for j in 0..4 { for i in 0..4 { out[(i, j)] *= inv_det.clone(); - if !out[(i, j)].is_finite() { - return false; - } } } true diff --git a/tests/linalg/inverse.rs b/tests/linalg/inverse.rs index 83659c640..aa4054fb2 100644 --- a/tests/linalg/inverse.rs +++ b/tests/linalg/inverse.rs @@ -141,6 +141,49 @@ fn matrix5_try_inverse_scaled_identity() { assert_relative_eq!(a_inv, expected_inverse); } +#[test] +#[rustfmt::skip] +fn matrix2_try_inverse_singular(){ + let m1 = Matrix2::new( + 1.0, 0.0, + 0.0,0.0 + ); + assert!(m1.try_inverse().is_none()); +} + +#[test] +#[rustfmt::skip] +fn matrix2_try_inverse_near_singular(){ + let m1 = Matrix2::new( + 1.0, 1e1, + 0.0, f64::MIN_POSITIVE); + if let Some(inv) = m1.try_inverse(){ + assert!(inv.iter().all(|x|x.is_finite())); + } +} + +#[test] +#[rustfmt::skip] +fn matrix3_try_inverse_singular(){ + let m1 = Matrix3::new( + 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0); + assert!(m1.try_inverse().is_none()); +} + +#[test] +#[rustfmt::skip] +fn matrix3_try_inverse_near_singular(){ + let m1 = Matrix3::new( + 1.0, 1e1, 1e1, + 0.0, 1.0, 1e1, + 0.0, 0.0, f64::MIN_POSITIVE); + if let Some(inv) = m1.try_inverse(){ + assert!(inv.iter().all(|x|x.is_finite())); + } +} + #[test] #[rustfmt::skip] fn matrix4_try_inverse_singular(){ @@ -160,5 +203,7 @@ fn matrix4_try_inverse_near_singular(){ 0.0, 1.0, 1e1, 1e1, 0.0, 0.0, 1.0, 1e1, 0.0, 0.0, 0.0, f64::MIN_POSITIVE ); - assert!(m1.try_inverse().is_none()); + if let Some(inv) = m1.try_inverse(){ + assert!(inv.iter().all(|x|x.is_finite())); + } } From fa45ee478ebc357318dc736beab56b3b3ca2f0c4 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Sat, 4 Jan 2025 00:43:31 +0100 Subject: [PATCH 7/8] remove whitespace --- src/linalg/inverse.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg/inverse.rs b/src/linalg/inverse.rs index 451ba8f16..cc4793f11 100644 --- a/src/linalg/inverse.rs +++ b/src/linalg/inverse.rs @@ -36,7 +36,7 @@ impl> SquareMatrix { impl> SquareMatrix { /// Attempts to invert this square matrix in-place. Returns `false` and leaves `self` untouched if /// inversion fails. - /// + /// /// Calling this function on near singular matrices may lead to`Inf` or `NaN` values! /// /// # Panics From 405d7e8010cda8be3b7ab03399e17502db9a0ae8 Mon Sep 17 00:00:00 2001 From: Mattis Kieffer <62772198+mat-kie@users.noreply.github.com> Date: Sat, 4 Jan 2025 00:44:39 +0100 Subject: [PATCH 8/8] remove unnececary assumption on array layout --- benches/lib.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benches/lib.rs b/benches/lib.rs index 429ba0e3b..d63b7b50d 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -24,7 +24,7 @@ fn reproducible_dmatrix(nrows: usize, ncols: usize) -> DMatrix { fn reproducible_matrix( ) -> Matrix, na::Const, na::ArrayStorage> where - [[T; R]; C]: Default, + na::ArrayStorage: Default, Standard: Distribution, { use rand::SeedableRng;