From 6ac6e7f75e17197ff9655b69d474fac05d8f61b3 Mon Sep 17 00:00:00 2001 From: metric-space Date: Wed, 22 Dec 2021 00:12:27 -0500 Subject: [PATCH 1/9] First compiling commit for take-2 of polar-decomposition: Code inspired by this thread: https://github.com/dimforge/nalgebra/pull/656 Main person behind this is LucasCampos --- src/linalg/svd.rs | 29 +++++++++++++++++++++++++++++ tests/linalg/svd.rs | 11 +++++++++++ 2 files changed, 40 insertions(+) diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 5fd3debc..dabcb491 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -643,6 +643,35 @@ where } } +impl, C: Dim> SVD +where + DefaultAllocator: Allocator, C> + + Allocator> + + Allocator>, +{ + /// converts SVD results to a polar form + + pub fn to_polar(&self) -> Result<(OMatrix, OMatrix), &'static str> + where DefaultAllocator: Allocator //result + + Allocator, R> // adjoint + + Allocator> // mapped vals + + Allocator // square matrix & result + + Allocator, DimMinimum> // ? + , + { + match (&self.u, &self.v_t) { + (Some(u), Some(v_t)) => Ok(( + u * OMatrix::from_diagonal(&self.singular_values.map(|e| T::from_real(e))) + * u.adjoint(), + u * v_t, + )), + (None, None) => Err("SVD solve: U and V^t have not been computed."), + (None, _) => Err("SVD solve: U has not been computed."), + (_, None) => Err("SVD solve: V^t has not been computed."), + } + } +} + impl, C: Dim> SVD where DimMinimum: DimSub, // for Bidiagonal. diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index d76d1759..65a92ddb 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -441,3 +441,14 @@ fn svd_sorted() { epsilon = 1.0e-5 ); } + +#[test] +fn svd_polar_decomposition() { + + let m = DMatrix::::new_random(4, 4); + let svd = m.clone().svd(true, true); + let (p,u) = svd.to_polar().unwrap(); + + assert_relative_eq!(m, p*u, epsilon = 1.0e-5); + +} From ac94fbe831ab4756d691f1e1e749b2b27a67390d Mon Sep 17 00:00:00 2001 From: metric-space Date: Sun, 26 Dec 2021 21:01:05 -0500 Subject: [PATCH 2/9] Add polar decomposition method to main matrix decomposition interface Add one more test for decomposition of polar decomposition of rectangular matrix --- src/linalg/decomposition.rs | 37 +++++++++++++++++++++++++++++++++++-- src/linalg/svd.rs | 29 +++++++++++------------------ tests/linalg/svd.rs | 14 ++++++++++++-- 3 files changed, 58 insertions(+), 22 deletions(-) diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index 91ad03d9..d75cae9c 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, RealField, Schur, SymmetricEigen, - SymmetricTridiagonal, LU, QR, SVD, U1, UDU, + DimMin, DimMinimum, DimSub, FullPivLU, Hessenberg, Matrix, OMatrix, RealField, Schur, + SymmetricEigen, SymmetricTridiagonal, LU, QR, SVD, U1, UDU, }; /// # Rectangular matrix decomposition @@ -17,6 +17,7 @@ use crate::{ /// | LU with partial pivoting | `P⁻¹ * L * U` | `L` is lower-triangular with a diagonal filled with `1` and `U` is upper-triangular. `P` is a permutation matrix. | /// | LU with full pivoting | `P⁻¹ * L * U * Q⁻¹` | `L` is lower-triangular with a diagonal filled with `1` and `U` is upper-triangular. `P` and `Q` are permutation matrices. | /// | SVD | `U * Σ * Vᵀ` | `U` and `V` are two orthogonal matrices and `Σ` is a diagonal matrix containing the singular values. | +/// | Polar (Left Polar) | `P' * U` | `U` is semi-unitary/unitary and `P'` is a positive semi-definite Hermitian Matrix impl> Matrix { /// Computes the bidiagonalization using householder reflections. pub fn bidiagonalize(self) -> Bidiagonal @@ -186,6 +187,38 @@ impl> Matrix { { SVD::try_new_unordered(self.into_owned(), compute_u, compute_v, eps, max_niter) } + + /// Attempts to compute the Polar Decomposition of a `matrix + /// + /// # Arguments + /// + /// * `eps` − tolerance used to determine when a value converged to 0. + /// * `max_niter` − maximum total number of iterations performed by the algorithm + pub fn polar( + self, + eps: T::RealField, + max_niter: usize, + ) -> Option<(OMatrix, OMatrix)> + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator, R> + + Allocator> + + Allocator + + Allocator, DimMinimum> + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>>, + { + SVD::try_new_unordered(self.into_owned(), true, true, eps, max_niter) + .and_then(|svd| svd.to_polar()) + } } /// # Square matrix decomposition diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index dabcb491..0e5d7f6c 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -641,33 +641,26 @@ where } } } -} -impl, C: Dim> SVD -where - DefaultAllocator: Allocator, C> - + Allocator> - + Allocator>, -{ - /// converts SVD results to a polar form - - pub fn to_polar(&self) -> Result<(OMatrix, OMatrix), &'static str> - where DefaultAllocator: Allocator //result + /// converts SVD results to Polar decomposition form of the original Matrix + /// A = P'U + /// The polar decomposition used here is Left Polar Decomposition (or Reverse Polar Decomposition) + /// Returns None if the SVD hasn't been calculated + pub fn to_polar(&self) -> Option<(OMatrix, OMatrix)> + where + DefaultAllocator: Allocator //result + Allocator, R> // adjoint + Allocator> // mapped vals - + Allocator // square matrix & result - + Allocator, DimMinimum> // ? - , + + Allocator // result + + Allocator, DimMinimum>, // square matrix { match (&self.u, &self.v_t) { - (Some(u), Some(v_t)) => Ok(( + (Some(u), Some(v_t)) => Some(( u * OMatrix::from_diagonal(&self.singular_values.map(|e| T::from_real(e))) * u.adjoint(), u * v_t, )), - (None, None) => Err("SVD solve: U and V^t have not been computed."), - (None, _) => Err("SVD solve: U has not been computed."), - (_, None) => Err("SVD solve: V^t has not been computed."), + _ => None, } } } diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 65a92ddb..251156b5 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -443,12 +443,22 @@ fn svd_sorted() { } #[test] -fn svd_polar_decomposition() { +fn dynamic_square_matrix_polar_decomposition() { - let m = DMatrix::::new_random(4, 4); + let m = DMatrix::::new_random(10, 10); let svd = m.clone().svd(true, true); let (p,u) = svd.to_polar().unwrap(); assert_relative_eq!(m, p*u, epsilon = 1.0e-5); } + +#[test] +fn dynamic_rectangular_matrix_polar_decomposition() { + + let m = DMatrix::::new_random(7, 5); + let svd = m.clone().svd(true, true); + let (p,u) = svd.to_polar().unwrap(); + + assert_relative_eq!(m, p*u, epsilon = 1.0e-5); +} From dbaefed8d14070b1425e1bd4be56e4f9e6490806 Mon Sep 17 00:00:00 2001 From: metric-space Date: Sun, 26 Dec 2021 21:05:42 -0500 Subject: [PATCH 3/9] Fix doc typos --- src/linalg/decomposition.rs | 2 +- src/linalg/svd.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index d75cae9c..e952f18b 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -188,7 +188,7 @@ impl> Matrix { SVD::try_new_unordered(self.into_owned(), compute_u, compute_v, eps, max_niter) } - /// Attempts to compute the Polar Decomposition of a `matrix + /// Attempts to compute the Polar Decomposition of a `matrix` (indirectly uses SVD) /// /// # Arguments /// diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 0e5d7f6c..70ab68a8 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -643,7 +643,7 @@ where } /// converts SVD results to Polar decomposition form of the original Matrix - /// A = P'U + /// A = P' * U /// The polar decomposition used here is Left Polar Decomposition (or Reverse Polar Decomposition) /// Returns None if the SVD hasn't been calculated pub fn to_polar(&self) -> Option<(OMatrix, OMatrix)> From 43c1f8fb9d931cbb7d3df2aac8dedbe644e483ba Mon Sep 17 00:00:00 2001 From: metric-space Date: Mon, 27 Dec 2021 02:12:54 -0500 Subject: [PATCH 4/9] Increased strength of tests for polar decomposition --- tests/linalg/svd.rs | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 251156b5..13c54a8b 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -449,7 +449,11 @@ fn dynamic_square_matrix_polar_decomposition() { let svd = m.clone().svd(true, true); let (p,u) = svd.to_polar().unwrap(); - assert_relative_eq!(m, p*u, epsilon = 1.0e-5); + assert_relative_eq!(m, &p*&u, epsilon = 1.0e-5); + // unitary check + assert_eq!(true, u.is_orthogonal(1.0e-5)); + // hermitian check + assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); } @@ -460,5 +464,9 @@ fn dynamic_rectangular_matrix_polar_decomposition() { let svd = m.clone().svd(true, true); let (p,u) = svd.to_polar().unwrap(); - assert_relative_eq!(m, p*u, epsilon = 1.0e-5); + assert_relative_eq!(m, &p*&u, epsilon = 1.0e-5); + // semi-unitary check + assert_eq!(true, (u.is_orthogonal(1.0e-5))); + // hermitian check + assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); } From cc10b67dd16653830de406e917877d3cd28b6827 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Thu, 30 Dec 2021 22:15:04 +0100 Subject: [PATCH 5/9] Add Matrix::try_polar that returns Option and make Matrix::polar not return Option --- src/linalg/decomposition.rs | 32 ++++++++++++++++++++++++++++---- src/linalg/svd.rs | 6 +++--- 2 files changed, 31 insertions(+), 7 deletions(-) diff --git a/src/linalg/decomposition.rs b/src/linalg/decomposition.rs index e952f18b..c72babf3 100644 --- a/src/linalg/decomposition.rs +++ b/src/linalg/decomposition.rs @@ -188,13 +188,37 @@ impl> Matrix { SVD::try_new_unordered(self.into_owned(), compute_u, compute_v, eps, max_niter) } - /// Attempts to compute the Polar Decomposition of a `matrix` (indirectly uses SVD) + /// Computes the Polar Decomposition of a `matrix` (indirectly uses SVD). + pub fn polar(self) -> (OMatrix, OMatrix) + where + R: DimMin, + DimMinimum: DimSub, // for Bidiagonal. + DefaultAllocator: Allocator + + Allocator, R> + + Allocator> + + Allocator + + Allocator, DimMinimum> + + Allocator + + Allocator + + Allocator, U1>> + + Allocator, C> + + Allocator> + + Allocator> + + Allocator> + + Allocator, U1>>, + { + SVD::new_unordered(self.into_owned(), true, true) + .to_polar() + .unwrap() + } + + /// Attempts to compute the Polar Decomposition of a `matrix` (indirectly uses SVD). /// /// # Arguments /// - /// * `eps` − tolerance used to determine when a value converged to 0. - /// * `max_niter` − maximum total number of iterations performed by the algorithm - pub fn polar( + /// * `eps` − tolerance used to determine when a value converged to 0 when computing the SVD. + /// * `max_niter` − maximum total number of iterations performed by the SVD computation algorithm. + pub fn try_polar( self, eps: T::RealField, max_niter: usize, diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 70ab68a8..3f945a65 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -642,10 +642,10 @@ where } } - /// converts SVD results to Polar decomposition form of the original Matrix - /// A = P' * U + /// converts SVD results to Polar decomposition form of the original Matrix: `A = P' * U`. + /// /// The polar decomposition used here is Left Polar Decomposition (or Reverse Polar Decomposition) - /// Returns None if the SVD hasn't been calculated + /// Returns None if the singular vectors of the SVD haven't been calculated pub fn to_polar(&self) -> Option<(OMatrix, OMatrix)> where DefaultAllocator: Allocator //result From 8e0ca439c2431c2e713e8355419a88e2d724f473 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Thu, 30 Dec 2021 22:15:22 +0100 Subject: [PATCH 6/9] Use proptest for testing the polar decomposition --- tests/linalg/svd.rs | 48 ++++++++++++++++++--------------------------- 1 file changed, 19 insertions(+), 29 deletions(-) diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 13c54a8b..0e725f8e 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -153,6 +153,25 @@ mod proptest_tests { prop_assert!(relative_eq!(&m * &sol2, b2, epsilon = 1.0e-6)); } } + + #[test] + fn svd_polar_decomposition(m in dmatrix_($scalar)) { + let svd = m.clone().svd(true, true); + let (p, u) = svd.to_polar().unwrap(); + + assert_relative_eq!(m, &p* &u, epsilon = 1.0e-5); + // semi-unitary check + assert_eq!(true, (u.is_orthogonal(1.0e-5))); + // hermitian check + assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); + + /* + * Same thing, but using the method instead of calling the SVD explicitly. + */ + let (p2, u2) = m.clone().polar(); + assert_eq!(p, p2); + assert_eq!(u, u2); + } } } } @@ -441,32 +460,3 @@ fn svd_sorted() { epsilon = 1.0e-5 ); } - -#[test] -fn dynamic_square_matrix_polar_decomposition() { - - let m = DMatrix::::new_random(10, 10); - let svd = m.clone().svd(true, true); - let (p,u) = svd.to_polar().unwrap(); - - assert_relative_eq!(m, &p*&u, epsilon = 1.0e-5); - // unitary check - assert_eq!(true, u.is_orthogonal(1.0e-5)); - // hermitian check - assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); - -} - -#[test] -fn dynamic_rectangular_matrix_polar_decomposition() { - - let m = DMatrix::::new_random(7, 5); - let svd = m.clone().svd(true, true); - let (p,u) = svd.to_polar().unwrap(); - - assert_relative_eq!(m, &p*&u, epsilon = 1.0e-5); - // semi-unitary check - assert_eq!(true, (u.is_orthogonal(1.0e-5))); - // hermitian check - assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); -} From 67a82c2c882b590c4b0ea91c1ef526030aafecd6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Thu, 30 Dec 2021 22:28:55 +0100 Subject: [PATCH 7/9] Test: minor style fix --- tests/linalg/svd.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 0e725f8e..6cba0d5d 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -161,7 +161,7 @@ mod proptest_tests { assert_relative_eq!(m, &p* &u, epsilon = 1.0e-5); // semi-unitary check - assert_eq!(true, (u.is_orthogonal(1.0e-5))); + assert!(u.is_orthogonal(1.0e-5)); // hermitian check assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); From ae6fda7dc77ccfcc2fe2dbdfd72c9f83bc43ab31 Mon Sep 17 00:00:00 2001 From: metric-space Date: Thu, 30 Dec 2021 21:11:39 -0500 Subject: [PATCH 8/9] Change svd to svd_unordered for the method output to be equal Comment out unitary check for now --- tests/linalg/svd.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index 6cba0d5d..fe93466d 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -156,12 +156,12 @@ mod proptest_tests { #[test] fn svd_polar_decomposition(m in dmatrix_($scalar)) { - let svd = m.clone().svd(true, true); + let svd = m.clone().svd_unordered(true, true); let (p, u) = svd.to_polar().unwrap(); assert_relative_eq!(m, &p* &u, epsilon = 1.0e-5); // semi-unitary check - assert!(u.is_orthogonal(1.0e-5)); + //assert!(u.is_orthogonal(1.0e-5)); // hermitian check assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5); From 498d7e3d4cb167bd4d3a53c3b12aeefe67d0b6f7 Mon Sep 17 00:00:00 2001 From: metric-space Date: Thu, 30 Dec 2021 21:18:58 -0500 Subject: [PATCH 9/9] Semi-unitary test checks for if rows or cols are orthonomal --- tests/linalg/svd.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs index fe93466d..deb3d38d 100644 --- a/tests/linalg/svd.rs +++ b/tests/linalg/svd.rs @@ -161,7 +161,7 @@ mod proptest_tests { assert_relative_eq!(m, &p* &u, epsilon = 1.0e-5); // semi-unitary check - //assert!(u.is_orthogonal(1.0e-5)); + assert!(u.is_orthogonal(1.0e-5) || u.transpose().is_orthogonal(1.0e-5)); // hermitian check assert_relative_eq!(p, p.adjoint(), epsilon = 1.0e-5);