From 8dda6714b5b8c515bc25004115125d622abe6f32 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Thu, 24 Sep 2020 23:21:13 -0600 Subject: [PATCH 01/14] Untested UDU implementation Pushing to trigger build Signed-off-by: Christopher Rabotin --- src/linalg/mod.rs | 2 ++ src/linalg/udu.rs | 68 +++++++++++++++++++++++++++++++++++++++++++++ tests/linalg/mod.rs | 1 + tests/linalg/udu.rs | 68 +++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 139 insertions(+) create mode 100644 src/linalg/udu.rs create mode 100644 tests/linalg/udu.rs diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index df9ae7de..075d91c2 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -25,6 +25,7 @@ mod solve; mod svd; mod symmetric_eigen; mod symmetric_tridiagonal; +mod udu; //// TODO: Not complete enough for publishing. //// This handles only cases where each eigenvalue has multiplicity one. @@ -45,3 +46,4 @@ pub use self::schur::*; pub use self::svd::*; pub use self::symmetric_eigen::*; pub use self::symmetric_tridiagonal::*; +pub use self::udu::*; diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs new file mode 100644 index 00000000..50c2e20d --- /dev/null +++ b/src/linalg/udu.rs @@ -0,0 +1,68 @@ +#[cfg(feature = "serde-serialize")] +use serde::{Deserialize, Serialize}; + +use crate::allocator::Allocator; +use crate::base::{DefaultAllocator, MatrixN}; +use crate::dimension::DimName; +use simba::scalar::ComplexField; + +/// UDU factorization +#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +#[derive(Clone, Debug)] +pub struct UDU +where + DefaultAllocator: Allocator, +{ + /// The upper triangular matrix resulting from the factorization + pub u: MatrixN, + /// The diagonal matrix resulting from the factorization + pub d: MatrixN, +} + +impl Copy for UDU +where + DefaultAllocator: Allocator, + MatrixN: Copy, +{ +} + +impl UDU +where + DefaultAllocator: Allocator, +{ + /// Computes the UDU^T factorization as per NASA's "Navigation Filter Best Practices", NTRS document ID 20180003657 + /// section 7.2 page 64. + /// NOTE: The provided matrix MUST be symmetric. + pub fn new(matrix: MatrixN) -> Self { + let mut d = MatrixN::::zeros(); + let mut u = MatrixN::::zeros(); + + let n = matrix.ncols(); + + d[(n, n)] = matrix[(n, n)]; + u[(n, n)] = N::one(); + + for j in (0..n - 1).rev() { + u[(j, n)] = matrix[(j, n)] / d[(n, n)]; + } + + for j in (1..n).rev() { + d[(j, j)] = matrix[(j, j)]; + for k in (j + 1..n).rev() { + d[(j, j)] = d[(j, j)] + d[(k, k)] * u[(j, k)].powi(2); + } + + u[(j, j)] = N::one(); + + for i in (0..j - 1).rev() { + u[(i, j)] = matrix[(i, j)]; + for k in j + 1..n { + u[(i, j)] = u[(i, j)] + d[(k, k)] * u[(i, k)] * u[(j, k)]; + } + u[(i, j)] /= d[(j, j)]; + } + } + + Self { u, d } + } +} diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index 66a0c06a..9c252bfd 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -14,3 +14,4 @@ mod schur; mod solve; mod svd; mod tridiagonal; +mod udu; diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs new file mode 100644 index 00000000..443178d5 --- /dev/null +++ b/tests/linalg/udu.rs @@ -0,0 +1,68 @@ +use na::Matrix3; +use na::UDU; + +#[test] +#[rustfmt::skip] +fn udu_simple() { + let m = Matrix3::new( + 2.0, -1.0, 0.0, + -1.0, 2.0, -1.0, + 0.0, -1.0, 2.0); + + let udu = UDU::new(m); + + println!("u = {}", udu.u); + println!("d = {}", udu.d); + + // Rebuild + let p = &udu.u * &udu.d * &udu.u.transpose(); + + println!("{}", p); + + assert!(relative_eq!(m, p, epsilon = 1.0e-7)); +} + +#[cfg(feature = "arbitrary")] +mod quickcheck_tests { + #[allow(unused_imports)] + use crate::core::helper::{RandComplex, RandScalar}; + + macro_rules! gen_tests( + ($module: ident, $scalar: ty) => { + mod $module { + use std::cmp; + use na::{DMatrix, Matrix4}; + #[allow(unused_imports)] + use crate::core::helper::{RandScalar, RandComplex}; + + quickcheck! { + fn udu(m: DMatrix<$scalar>) -> bool { + let mut m = m; + if m.len() == 0 { + m = DMatrix::<$scalar>::new_random(1, 1); + } + + let m = m.map(|e| e.0); + + let udu = UDU::new(m.clone()); + let p = &udu.u * &udu.d * &udu.u.transpose(); + + relative_eq!(m, p, epsilon = 1.0e-7) + } + + fn udu_static(m: Matrix4<$scalar>) -> bool { + let m = m.map(|e| e.0); + + let udu = UDU::new(m.clone()); + let p = &udu.u * &udu.d * &udu.u.transpose(); + + relative_eq!(m, p, epsilon = 1.0e-7) + } + } + } + } + ); + + gen_tests!(complex, RandComplex); + gen_tests!(f64, RandScalar); +} From d534c3bf9d0d5ea4515bc06a9ea908336cdf117a Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Thu, 24 Sep 2020 23:46:24 -0600 Subject: [PATCH 02/14] Trying to break the test to make sure it works Signed-off-by: Christopher Rabotin --- tests/linalg/udu.rs | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs index 443178d5..d9d73e81 100644 --- a/tests/linalg/udu.rs +++ b/tests/linalg/udu.rs @@ -10,16 +10,10 @@ fn udu_simple() { 0.0, -1.0, 2.0); let udu = UDU::new(m); - - println!("u = {}", udu.u); - println!("d = {}", udu.d); - // Rebuild let p = &udu.u * &udu.d * &udu.u.transpose(); - println!("{}", p); - - assert!(relative_eq!(m, p, epsilon = 1.0e-7)); + assert!(relative_eq!(m, 2.0*p, epsilon = 1.0e-7)); } #[cfg(feature = "arbitrary")] From a8d40423ea248d318838cfa941085b0a7883e4f0 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Fri, 25 Sep 2020 19:21:14 -0600 Subject: [PATCH 03/14] Fixed UDU algorithm Signed-off-by: Christopher Rabotin --- src/linalg/udu.rs | 35 +++++++++++++++++++---------------- tests/linalg/udu.rs | 8 ++++---- 2 files changed, 23 insertions(+), 20 deletions(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index 50c2e20d..30143e44 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -30,37 +30,40 @@ impl UDU where DefaultAllocator: Allocator, { - /// Computes the UDU^T factorization as per NASA's "Navigation Filter Best Practices", NTRS document ID 20180003657 - /// section 7.2 page 64. - /// NOTE: The provided matrix MUST be symmetric. - pub fn new(matrix: MatrixN) -> Self { + /// Computes the UDU^T factorization + /// NOTE: The provided matrix MUST be symmetric, and no verification is done in this regard. + /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360 + pub fn new(p: MatrixN) -> Self { let mut d = MatrixN::::zeros(); let mut u = MatrixN::::zeros(); - let n = matrix.ncols(); + let n = p.ncols() - 1; - d[(n, n)] = matrix[(n, n)]; + d[(n, n)] = p[(n, n)]; u[(n, n)] = N::one(); - for j in (0..n - 1).rev() { - u[(j, n)] = matrix[(j, n)] / d[(n, n)]; + for j in (0..n).rev() { + u[(j, n)] = p[(j, n)] / d[(n, n)]; } - for j in (1..n).rev() { - d[(j, j)] = matrix[(j, j)]; - for k in (j + 1..n).rev() { + for j in (0..n).rev() { + for k in j + 1..=n { d[(j, j)] = d[(j, j)] + d[(k, k)] * u[(j, k)].powi(2); } - u[(j, j)] = N::one(); + d[(j, j)] = p[(j, j)] - d[(j, j)]; - for i in (0..j - 1).rev() { - u[(i, j)] = matrix[(i, j)]; - for k in j + 1..n { - u[(i, j)] = u[(i, j)] + d[(k, k)] * u[(i, k)] * u[(j, k)]; + for i in (0..=j).rev() { + for k in j + 1..=n { + u[(i, j)] = u[(i, j)] + d[(k, k)] * u[(j, k)] * u[(i, k)]; } + + u[(i, j)] = p[(i, j)] - u[(i, j)]; + u[(i, j)] /= d[(j, j)]; } + + u[(j, j)] = N::one(); } Self { u, d } diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs index d9d73e81..0d457db5 100644 --- a/tests/linalg/udu.rs +++ b/tests/linalg/udu.rs @@ -11,9 +11,9 @@ fn udu_simple() { let udu = UDU::new(m); // Rebuild - let p = &udu.u * &udu.d * &udu.u.transpose(); + let p = udu.u * udu.d * udu.u.transpose(); - assert!(relative_eq!(m, 2.0*p, epsilon = 1.0e-7)); + assert!(relative_eq!(m, p, epsilon = 3.0e-16)); } #[cfg(feature = "arbitrary")] @@ -48,9 +48,9 @@ mod quickcheck_tests { let m = m.map(|e| e.0); let udu = UDU::new(m.clone()); - let p = &udu.u * &udu.d * &udu.u.transpose(); + let p = udu.u * udu.d * udu.u.transpose(); - relative_eq!(m, p, epsilon = 1.0e-7) + relative_eq!(m, p, epsilon = 3.0e-16) } } } From 5a7ed61e9b82134251ae37013acab293c3633674 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Fri, 25 Sep 2020 21:29:46 -0600 Subject: [PATCH 04/14] UDU impl: using 0-index nomenclature Signed-off-by: Christopher Rabotin --- src/linalg/udu.rs | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index 30143e44..cd2db206 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -37,24 +37,24 @@ where let mut d = MatrixN::::zeros(); let mut u = MatrixN::::zeros(); - let n = p.ncols() - 1; + let n = p.ncols(); - d[(n, n)] = p[(n, n)]; - u[(n, n)] = N::one(); + d[(n - 1, n - 1)] = p[(n - 1, n - 1)]; + u[(n - 1, n - 1)] = N::one(); - for j in (0..n).rev() { - u[(j, n)] = p[(j, n)] / d[(n, n)]; + for j in (0..n - 1).rev() { + u[(j, n - 1)] = p[(j, n - 1)] / d[(n - 1, n - 1)]; } - for j in (0..n).rev() { - for k in j + 1..=n { + for j in (0..n - 1).rev() { + for k in j + 1..n { d[(j, j)] = d[(j, j)] + d[(k, k)] * u[(j, k)].powi(2); } d[(j, j)] = p[(j, j)] - d[(j, j)]; for i in (0..=j).rev() { - for k in j + 1..=n { + for k in j + 1..n { u[(i, j)] = u[(i, j)] + d[(k, k)] * u[(j, k)] * u[(i, k)]; } From e9933e5c91c4fcd9693d2afd80bb0f7c43e4abf8 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sat, 26 Sep 2020 18:34:35 -0600 Subject: [PATCH 05/14] UDU: Expand to Dim from DimName Signed-off-by: Christopher Rabotin --- src/linalg/udu.rs | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index cd2db206..eada8e99 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -3,13 +3,13 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, MatrixN}; -use crate::dimension::DimName; +use crate::dimension::Dim; use simba::scalar::ComplexField; /// UDU factorization #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[derive(Clone, Debug)] -pub struct UDU +pub struct UDU where DefaultAllocator: Allocator, { @@ -19,14 +19,14 @@ where pub d: MatrixN, } -impl Copy for UDU +impl Copy for UDU where DefaultAllocator: Allocator, MatrixN: Copy, { } -impl UDU +impl UDU where DefaultAllocator: Allocator, { @@ -34,10 +34,11 @@ where /// NOTE: The provided matrix MUST be symmetric, and no verification is done in this regard. /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360 pub fn new(p: MatrixN) -> Self { - let mut d = MatrixN::::zeros(); - let mut u = MatrixN::::zeros(); - let n = p.ncols(); + let n_as_dim = D::from_usize(n); + + let mut d = MatrixN::::zeros_generic(n_as_dim, n_as_dim); + let mut u = MatrixN::::zeros_generic(n_as_dim, n_as_dim); d[(n - 1, n - 1)] = p[(n - 1, n - 1)]; u[(n - 1, n - 1)] = N::one(); From 7a49b9eecaac3e2075b8a91b08ee70a9d4e667ce Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 27 Sep 2020 15:28:50 -0600 Subject: [PATCH 06/14] UDU: d now stored in VectorN instead of MatrixN Signed-off-by: Christopher Rabotin --- src/linalg/udu.rs | 30 ++++++++++++++++++------------ tests/linalg/udu.rs | 6 +++--- 2 files changed, 21 insertions(+), 15 deletions(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index eada8e99..953359e8 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -2,7 +2,7 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; -use crate::base::{DefaultAllocator, MatrixN}; +use crate::base::{DefaultAllocator, MatrixN, VectorN, U1}; use crate::dimension::Dim; use simba::scalar::ComplexField; @@ -11,24 +11,25 @@ use simba::scalar::ComplexField; #[derive(Clone, Debug)] pub struct UDU where - DefaultAllocator: Allocator, + DefaultAllocator: Allocator + Allocator, { /// The upper triangular matrix resulting from the factorization pub u: MatrixN, /// The diagonal matrix resulting from the factorization - pub d: MatrixN, + pub d: VectorN, } impl Copy for UDU where - DefaultAllocator: Allocator, + DefaultAllocator: Allocator + Allocator, + VectorN: Copy, MatrixN: Copy, { } impl UDU where - DefaultAllocator: Allocator, + DefaultAllocator: Allocator + Allocator, { /// Computes the UDU^T factorization /// NOTE: The provided matrix MUST be symmetric, and no verification is done in this regard. @@ -37,31 +38,31 @@ where let n = p.ncols(); let n_as_dim = D::from_usize(n); - let mut d = MatrixN::::zeros_generic(n_as_dim, n_as_dim); + let mut d = VectorN::::zeros_generic(n_as_dim, U1); let mut u = MatrixN::::zeros_generic(n_as_dim, n_as_dim); - d[(n - 1, n - 1)] = p[(n - 1, n - 1)]; + d[n - 1] = p[(n - 1, n - 1)]; u[(n - 1, n - 1)] = N::one(); for j in (0..n - 1).rev() { - u[(j, n - 1)] = p[(j, n - 1)] / d[(n - 1, n - 1)]; + u[(j, n - 1)] = p[(j, n - 1)] / d[n - 1]; } for j in (0..n - 1).rev() { for k in j + 1..n { - d[(j, j)] = d[(j, j)] + d[(k, k)] * u[(j, k)].powi(2); + d[j] = d[j] + d[k] * u[(j, k)].powi(2); } - d[(j, j)] = p[(j, j)] - d[(j, j)]; + d[j] = p[(j, j)] - d[j]; for i in (0..=j).rev() { for k in j + 1..n { - u[(i, j)] = u[(i, j)] + d[(k, k)] * u[(j, k)] * u[(i, k)]; + u[(i, j)] = u[(i, j)] + d[k] * u[(j, k)] * u[(i, k)]; } u[(i, j)] = p[(i, j)] - u[(i, j)]; - u[(i, j)] /= d[(j, j)]; + u[(i, j)] /= d[j]; } u[(j, j)] = N::one(); @@ -69,4 +70,9 @@ where Self { u, d } } + + /// Returns the diagonal elements as a matrix + pub fn d_matrix(&self) -> MatrixN { + MatrixN::from_diagonal(&self.d) + } } diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs index 0d457db5..3304d73b 100644 --- a/tests/linalg/udu.rs +++ b/tests/linalg/udu.rs @@ -11,7 +11,7 @@ fn udu_simple() { let udu = UDU::new(m); // Rebuild - let p = udu.u * udu.d * udu.u.transpose(); + let p = udu.u * udu.d_matrix() * udu.u.transpose(); assert!(relative_eq!(m, p, epsilon = 3.0e-16)); } @@ -39,7 +39,7 @@ mod quickcheck_tests { let m = m.map(|e| e.0); let udu = UDU::new(m.clone()); - let p = &udu.u * &udu.d * &udu.u.transpose(); + let p = &udu.u * &udu.d_matrix() * &udu.u.transpose(); relative_eq!(m, p, epsilon = 1.0e-7) } @@ -48,7 +48,7 @@ mod quickcheck_tests { let m = m.map(|e| e.0); let udu = UDU::new(m.clone()); - let p = udu.u * udu.d * udu.u.transpose(); + let p = udu.u * udu.d_matrix() * udu.u.transpose(); relative_eq!(m, p, epsilon = 3.0e-16) } From f6c1aeb07f011a7c1cda295eaae09b54cf17657b Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 27 Sep 2020 16:18:47 -0600 Subject: [PATCH 07/14] UDU: add panic test for non symmetric matrix Signed-off-by: Christopher Rabotin --- tests/linalg/udu.rs | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs index 3304d73b..c7f416d7 100644 --- a/tests/linalg/udu.rs +++ b/tests/linalg/udu.rs @@ -16,6 +16,22 @@ fn udu_simple() { assert!(relative_eq!(m, p, epsilon = 3.0e-16)); } +#[test] +#[should_panic] +#[rustfmt::skip] +fn udu_non_sym_panic() { + let m = Matrix3::new( + 2.0, -1.0, 0.0, + 1.0, -2.0, 3.0, + -2.0, 1.0, 0.0); + + let udu = UDU::new(m); + // Rebuild + let p = udu.u * udu.d_matrix() * udu.u.transpose(); + + assert!(relative_eq!(m, p, epsilon = 3.0e-16)); +} + #[cfg(feature = "arbitrary")] mod quickcheck_tests { #[allow(unused_imports)] From 06861a9755e2f629bec0b25371affd7d0f8a697d Mon Sep 17 00:00:00 2001 From: Chris Date: Mon, 26 Oct 2020 14:39:34 -0600 Subject: [PATCH 08/14] Update src/linalg/udu.rs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Sébastien Crozet --- src/linalg/udu.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index 953359e8..182b8a6d 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -32,7 +32,7 @@ where DefaultAllocator: Allocator + Allocator, { /// Computes the UDU^T factorization - /// NOTE: The provided matrix MUST be symmetric, and no verification is done in this regard. + /// The input matrix `p` is assumed to be symmetric and this decomposition will only read the upper-triangular part of `p`. /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360 pub fn new(p: MatrixN) -> Self { let n = p.ncols(); From 4ff4911ac3559bc7e825a8bd2b05d1e0e475e69d Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Mon, 26 Oct 2020 22:06:37 -0600 Subject: [PATCH 09/14] Implement requested changes Signed-off-by: Christopher Rabotin --- src/linalg/udu.rs | 36 ++++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index 182b8a6d..bce2ea41 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -4,10 +4,21 @@ use serde::{Deserialize, Serialize}; use crate::allocator::Allocator; use crate::base::{DefaultAllocator, MatrixN, VectorN, U1}; use crate::dimension::Dim; +use crate::storage::Storage; use simba::scalar::ComplexField; /// UDU factorization #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] +#[cfg_attr( + feature = "serde-serialize", + serde(bound(serialize = "VectorN: Serialize, MatrixN: Serialize")) +)] +#[cfg_attr( + feature = "serde-serialize", + serde(bound( + deserialize = "VectorN: Deserialize<'de>, MatrixN: Deserialize<'de>" + )) +)] #[derive(Clone, Debug)] pub struct UDU where @@ -36,33 +47,30 @@ where /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360 pub fn new(p: MatrixN) -> Self { let n = p.ncols(); - let n_as_dim = D::from_usize(n); + let n_dim = p.data.shape().1; - let mut d = VectorN::::zeros_generic(n_as_dim, U1); - let mut u = MatrixN::::zeros_generic(n_as_dim, n_as_dim); + let mut d = VectorN::zeros_generic(n_dim, U1); + let mut u = MatrixN::zeros_generic(n_dim, n_dim); d[n - 1] = p[(n - 1, n - 1)]; - u[(n - 1, n - 1)] = N::one(); - - for j in (0..n - 1).rev() { - u[(j, n - 1)] = p[(j, n - 1)] / d[n - 1]; - } + u.column_mut(n - 1) + .axpy(N::one() / d[n - 1], &p.column(n - 1), N::zero()); for j in (0..n - 1).rev() { + let mut d_j = d[j]; for k in j + 1..n { - d[j] = d[j] + d[k] * u[(j, k)].powi(2); + d_j += d[k] * u[(j, k)].powi(2); } - d[j] = p[(j, j)] - d[j]; + d[j] = p[(j, j)] - d_j; for i in (0..=j).rev() { + let mut u_ij = u[(i, j)]; for k in j + 1..n { - u[(i, j)] = u[(i, j)] + d[k] * u[(j, k)] * u[(i, k)]; + u_ij += d[k] * u[(j, k)] * u[(i, k)]; } - u[(i, j)] = p[(i, j)] - u[(i, j)]; - - u[(i, j)] /= d[j]; + u[(i, j)] = (p[(i, j)] - u_ij) / d[j]; } u[(j, j)] = N::one(); From 89ca2fe5fbe7a8987096cb5d8ca0fe6da92189e8 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Wed, 28 Oct 2020 16:04:46 -0600 Subject: [PATCH 10/14] UDU only supported for Real matrices, not Complex Signed-off-by: Christopher Rabotin --- src/linalg/udu.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index bce2ea41..46b1497c 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -5,7 +5,7 @@ use crate::allocator::Allocator; use crate::base::{DefaultAllocator, MatrixN, VectorN, U1}; use crate::dimension::Dim; use crate::storage::Storage; -use simba::scalar::ComplexField; +use simba::scalar::RealField; /// UDU factorization #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] @@ -20,7 +20,7 @@ use simba::scalar::ComplexField; )) )] #[derive(Clone, Debug)] -pub struct UDU +pub struct UDU where DefaultAllocator: Allocator + Allocator, { @@ -30,7 +30,7 @@ where pub d: VectorN, } -impl Copy for UDU +impl Copy for UDU where DefaultAllocator: Allocator + Allocator, VectorN: Copy, @@ -38,7 +38,7 @@ where { } -impl UDU +impl UDU where DefaultAllocator: Allocator + Allocator, { From ab0d335b61f9869d8a958c7ac84ea80ecba8d42e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Crozet=20S=C3=A9bastien?= Date: Thu, 25 Feb 2021 13:16:04 +0100 Subject: [PATCH 11/14] Fix tests for the UDU decomposition. --- src/linalg/udu.rs | 9 ++++++--- tests/linalg/udu.rs | 21 +++++++-------------- 2 files changed, 13 insertions(+), 17 deletions(-) diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs index 46b1497c..346753b2 100644 --- a/src/linalg/udu.rs +++ b/src/linalg/udu.rs @@ -7,7 +7,7 @@ use crate::dimension::Dim; use crate::storage::Storage; use simba::scalar::RealField; -/// UDU factorization +/// UDU factorization. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize", @@ -42,8 +42,11 @@ impl UDU where DefaultAllocator: Allocator + Allocator, { - /// Computes the UDU^T factorization - /// The input matrix `p` is assumed to be symmetric and this decomposition will only read the upper-triangular part of `p`. + /// Computes the UDU^T factorization. + /// + /// The input matrix `p` is assumed to be symmetric and this decomposition will only read + /// the upper-triangular part of `p`. + /// /// Ref.: "Optimal control and estimation-Dover Publications", Robert F. Stengel, (1994) page 360 pub fn new(p: MatrixN) -> Self { let n = p.ncols(); diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs index c7f416d7..95cdc9c2 100644 --- a/tests/linalg/udu.rs +++ b/tests/linalg/udu.rs @@ -1,5 +1,4 @@ -use na::Matrix3; -use na::UDU; +use na::{Matrix3, UDU}; #[test] #[rustfmt::skip] @@ -40,19 +39,14 @@ mod quickcheck_tests { macro_rules! gen_tests( ($module: ident, $scalar: ty) => { mod $module { - use std::cmp; - use na::{DMatrix, Matrix4}; + use na::{UDU, DMatrix, Matrix4}; #[allow(unused_imports)] use crate::core::helper::{RandScalar, RandComplex}; quickcheck! { - fn udu(m: DMatrix<$scalar>) -> bool { - let mut m = m; - if m.len() == 0 { - m = DMatrix::<$scalar>::new_random(1, 1); - } - - let m = m.map(|e| e.0); + fn udu(n: usize) -> bool { + let n = std::cmp::max(1, std::cmp::min(n, 10)); + let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part(); let udu = UDU::new(m.clone()); let p = &udu.u * &udu.d_matrix() * &udu.u.transpose(); @@ -61,18 +55,17 @@ mod quickcheck_tests { } fn udu_static(m: Matrix4<$scalar>) -> bool { - let m = m.map(|e| e.0); + let m = m.map(|e| e.0).hermitian_part(); let udu = UDU::new(m.clone()); let p = udu.u * udu.d_matrix() * udu.u.transpose(); - relative_eq!(m, p, epsilon = 3.0e-16) + relative_eq!(m, p, epsilon = 1.0e-7) } } } } ); - gen_tests!(complex, RandComplex); gen_tests!(f64, RandScalar); } From aeb9f7ea3994499540e3a7eeb1b4eed46768a275 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Crozet=20S=C3=A9bastien?= Date: Thu, 25 Feb 2021 13:20:20 +0100 Subject: [PATCH 12/14] Add a matrix.udu() method to compute the UDU decomposition. --- src/linalg/decomposition.rs | 17 +++++++++++++++-- tests/linalg/udu.rs | 13 +++++++------ 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 7890748c..99f2bcfe 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -1,8 +1,8 @@ use crate::storage::Storage; use crate::{ Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff, - DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, Schur, SymmetricEigen, - SymmetricTridiagonal, LU, QR, SVD, U1, + DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, RealField, Schur, SymmetricEigen, + SymmetricTridiagonal, LU, QR, SVD, U1, UDU }; /// # Rectangular matrix decomposition @@ -134,6 +134,7 @@ impl> Matrix { /// | -------------------------|---------------------------|--------------| /// | Hessenberg | `Q * H * Qᵀ` | `Q` is a unitary matrix and `H` an upper-Hessenberg matrix. | /// | Cholesky | `L * Lᵀ` | `L` is a lower-triangular matrix. | +/// | UDU | `U * D * Uᵀ` | `U` is a upper-triangular matrix, and `D` a diagonal matrix. | /// | Schur decomposition | `Q * T * Qᵀ` | `Q` is an unitary matrix and `T` a quasi-upper-triangular matrix. | /// | Symmetric eigendecomposition | `Q ~ Λ ~ Qᵀ` | `Q` is an unitary matrix, and `Λ` is a real diagonal matrix. | /// | Symmetric tridiagonalization | `Q ~ T ~ Qᵀ` | `Q` is an unitary matrix, and `T` is a tridiagonal matrix. | @@ -149,6 +150,18 @@ impl> Matrix { Cholesky::new(self.into_owned()) } + /// Attempts to compute the UDU decomposition of this matrix. + /// + /// The input matrix `self` is assumed to be symmetric and this decomposition will only read + /// the upper-triangular part of `self`. + pub fn udu(self) -> UDU + where + N: RealField, + DefaultAllocator: Allocator + Allocator, + { + UDU::new(self.into_owned()) + } + /// Computes the Hessenberg decomposition of this matrix using householder reflections. pub fn hessenberg(self) -> Hessenberg where diff --git a/tests/linalg/udu.rs b/tests/linalg/udu.rs index 95cdc9c2..006cc95f 100644 --- a/tests/linalg/udu.rs +++ b/tests/linalg/udu.rs @@ -1,4 +1,4 @@ -use na::{Matrix3, UDU}; +use na::Matrix3; #[test] #[rustfmt::skip] @@ -8,7 +8,8 @@ fn udu_simple() { -1.0, 2.0, -1.0, 0.0, -1.0, 2.0); - let udu = UDU::new(m); + let udu = m.udu(); + // Rebuild let p = udu.u * udu.d_matrix() * udu.u.transpose(); @@ -24,7 +25,7 @@ fn udu_non_sym_panic() { 1.0, -2.0, 3.0, -2.0, 1.0, 0.0); - let udu = UDU::new(m); + let udu = m.udu(); // Rebuild let p = udu.u * udu.d_matrix() * udu.u.transpose(); @@ -39,7 +40,7 @@ mod quickcheck_tests { macro_rules! gen_tests( ($module: ident, $scalar: ty) => { mod $module { - use na::{UDU, DMatrix, Matrix4}; + use na::{DMatrix, Matrix4}; #[allow(unused_imports)] use crate::core::helper::{RandScalar, RandComplex}; @@ -48,7 +49,7 @@ mod quickcheck_tests { let n = std::cmp::max(1, std::cmp::min(n, 10)); let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0).hermitian_part(); - let udu = UDU::new(m.clone()); + let udu = m.clone().udu(); let p = &udu.u * &udu.d_matrix() * &udu.u.transpose(); relative_eq!(m, p, epsilon = 1.0e-7) @@ -57,7 +58,7 @@ mod quickcheck_tests { fn udu_static(m: Matrix4<$scalar>) -> bool { let m = m.map(|e| e.0).hermitian_part(); - let udu = UDU::new(m.clone()); + let udu = m.udu(); let p = udu.u * udu.d_matrix() * udu.u.transpose(); relative_eq!(m, p, epsilon = 1.0e-7) From 6699039fecbb9271539432d0dfbeac822d50c9ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Crozet=20S=C3=A9bastien?= Date: Thu, 25 Feb 2021 15:51:13 +0100 Subject: [PATCH 13/14] Fix rebase-induced compilation error. --- src/linalg/col_piv_qr.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/linalg/col_piv_qr.rs b/src/linalg/col_piv_qr.rs index 4fa97482..302dcd66 100644 --- a/src/linalg/col_piv_qr.rs +++ b/src/linalg/col_piv_qr.rs @@ -66,7 +66,8 @@ where let min_nrows_ncols = nrows.min(ncols); let mut p = PermutationSequence::identity_generic(min_nrows_ncols); - let mut diag = unsafe { MatrixMN::new_uninitialized_generic(min_nrows_ncols, U1) }; + let mut diag = + unsafe { crate::unimplemented_or_uninitialized_generic!(min_nrows_ncols, U1) }; if min_nrows_ncols.value() == 0 { return ColPivQR { From 7b6b3649f27770fe456add8debff443e87d4147f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Crozet=20S=C3=A9bastien?= Date: Thu, 25 Feb 2021 16:20:11 +0100 Subject: [PATCH 14/14] Run cargo fmt. --- src/linalg/decomposition.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 99f2bcfe..e167861a 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -2,7 +2,7 @@ use crate::storage::Storage; use crate::{ Allocator, Bidiagonal, Cholesky, ColPivQR, ComplexField, DefaultAllocator, Dim, DimDiff, DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, RealField, Schur, SymmetricEigen, - SymmetricTridiagonal, LU, QR, SVD, U1, UDU + SymmetricTridiagonal, LU, QR, SVD, U1, UDU, }; /// # Rectangular matrix decomposition