From 308d95386e1f775f437001331eecf2393d992668 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Crozet=20S=C3=A9bastien?= Date: Thu, 25 Feb 2021 12:06:04 +0100 Subject: [PATCH] Fix all tests and the ColPivQR::solve. --- src/linalg/{qrp.rs => col_piv_qr.rs} | 200 +++++++++++++------------ src/linalg/householder.rs | 2 +- src/linalg/mod.rs | 2 + src/linalg/qr.rs | 4 +- tests/linalg/{qrp.rs => col_piv_qr.rs} | 163 ++++++++++---------- tests/linalg/mod.rs | 1 + 6 files changed, 198 insertions(+), 174 deletions(-) rename src/linalg/{qrp.rs => col_piv_qr.rs} (60%) rename tests/linalg/{qrp.rs => col_piv_qr.rs} (52%) diff --git a/src/linalg/qrp.rs b/src/linalg/col_piv_qr.rs similarity index 60% rename from src/linalg/qrp.rs rename to src/linalg/col_piv_qr.rs index c7d2f27b..3f8ecda1 100644 --- a/src/linalg/qrp.rs +++ b/src/linalg/col_piv_qr.rs @@ -1,98 +1,94 @@ +use num::Zero; #[cfg(feature = "serde-serialize")] use serde::{Deserialize, Serialize}; -use num::Zero; -use alga::general::ComplexField; use crate::allocator::{Allocator, Reallocator}; use crate::base::{DefaultAllocator, Matrix, MatrixMN, MatrixN, Unit, VectorN}; use crate::constraint::{SameNumberOfRows, ShapeConstraint}; use crate::dimension::{Dim, DimMin, DimMinimum, U1}; use crate::storage::{Storage, StorageMut}; +use crate::ComplexField; use crate::geometry::Reflection; -use crate::linalg::householder; -use crate::linalg::PermutationSequence; +use crate::linalg::{householder, PermutationSequence}; -/// The QRP decomposition of a general matrix. +/// The QR decomposition (with column pivoting) of a general matrix. #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[cfg_attr( feature = "serde-serialize", - serde(bound( - serialize = "DefaultAllocator: Allocator + + serde(bound(serialize = "DefaultAllocator: Allocator + Allocator>, MatrixMN: Serialize, PermutationSequence>: Serialize, - VectorN>: Serialize" - )) + VectorN>: Serialize")) )] #[cfg_attr( feature = "serde-serialize", - serde(bound( - deserialize = "DefaultAllocator: Allocator + + serde(bound(deserialize = "DefaultAllocator: Allocator + Allocator>, MatrixMN: Deserialize<'de>, PermutationSequence>: Deserialize<'de>, - VectorN>: Deserialize<'de>" - )) + VectorN>: Deserialize<'de>")) )] #[derive(Clone, Debug)] -pub struct QRP, C: Dim> -where DefaultAllocator: Allocator + Allocator> + - Allocator<(usize, usize), DimMinimum>, +pub struct ColPivQR, C: Dim> +where + DefaultAllocator: Allocator + + Allocator> + + Allocator<(usize, usize), DimMinimum>, { - qrp: MatrixMN, + col_piv_qr: MatrixMN, p: PermutationSequence>, diag: VectorN>, } -impl, C: Dim> Copy for QRP +impl, C: Dim> Copy for ColPivQR where - DefaultAllocator: Allocator + Allocator> + - Allocator<(usize, usize), DimMinimum>, + DefaultAllocator: Allocator + + Allocator> + + Allocator<(usize, usize), DimMinimum>, MatrixMN: Copy, PermutationSequence>: Copy, VectorN>: Copy, -{} - -impl, C: Dim> QRP -where DefaultAllocator: Allocator + Allocator + Allocator> + Allocator<(usize, usize), DimMinimum> { - /// Computes the QRP decomposition using householder reflections. +} + +impl, C: Dim> ColPivQR +where + DefaultAllocator: Allocator + + Allocator + + Allocator> + + Allocator<(usize, usize), DimMinimum>, +{ + /// Computes the ColPivQR decomposition using householder reflections. pub fn new(mut matrix: MatrixMN) -> Self { let (nrows, ncols) = matrix.data.shape(); 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) }; - println!("diag: {:?}", &diag); if min_nrows_ncols.value() == 0 { - return QRP { - qrp: matrix, - p: p, - diag: diag, + return ColPivQR { + col_piv_qr: matrix, + p, + diag, }; } - for ite in 0..min_nrows_ncols.value() { - let mut col_norm = Vec::new(); - for column in matrix.column_iter() { - col_norm.push(column.norm_squared()); - } - - let piv = matrix.slice_range(ite.., ite..).icamax_full(); - let col_piv = piv.1 + ite; - matrix.swap_columns(ite, col_piv); - p.append_permutation(ite, col_piv); + for i in 0..min_nrows_ncols.value() { + let piv = matrix.slice_range(i.., i..).icamax_full(); + let col_piv = piv.1 + i; + matrix.swap_columns(i, col_piv); + p.append_permutation(i, col_piv); - householder::clear_column_unchecked(&mut matrix, &mut diag[ite], ite, 0, None); - println!("matrix: {:?}", &matrix.data); + householder::clear_column_unchecked(&mut matrix, &mut diag[i], i, 0, None); } - QRP { - qrp: matrix, - p: p, - diag: diag, + ColPivQR { + col_piv_qr: matrix, + p, + diag, } } @@ -102,8 +98,11 @@ where DefaultAllocator: Allocator + Allocator + Allocator, C>, { - let (nrows, ncols) = self.qrp.data.shape(); - let mut res = self.qrp.rows_generic(0, nrows.min(ncols)).upper_triangle(); + let (nrows, ncols) = self.col_piv_qr.data.shape(); + let mut res = self + .col_piv_qr + .rows_generic(0, nrows.min(ncols)) + .upper_triangle(); res.set_partial_diagonal(self.diag.iter().map(|e| N::from_real(e.modulus()))); res } @@ -116,8 +115,10 @@ where DefaultAllocator: Allocator + Allocator + Allocator, C>, { - let (nrows, ncols) = self.qrp.data.shape(); - let mut res = self.qrp.resize_generic(nrows.min(ncols), ncols, N::zero()); + let (nrows, ncols) = self.col_piv_qr.data.shape(); + let mut res = self + .col_piv_qr + .resize_generic(nrows.min(ncols), ncols, N::zero()); res.fill_lower_triangle(N::zero(), 1); res.set_partial_diagonal(self.diag.iter().map(|e| N::from_real(e.modulus()))); res @@ -125,8 +126,10 @@ where DefaultAllocator: Allocator + Allocator + Allocator MatrixMN> - where DefaultAllocator: Allocator> { - let (nrows, ncols) = self.qrp.data.shape(); + where + DefaultAllocator: Allocator>, + { + let (nrows, ncols) = self.col_piv_qr.data.shape(); // NOTE: we could build the identity matrix and call q_mul on it. // Instead we don't so that we take in account the matrix sparseness. @@ -134,8 +137,8 @@ where DefaultAllocator: Allocator + Allocator + Allocator + Allocator + Allocator: DimMin>, - DefaultAllocator: - Allocator> + Reallocator, C> + Allocator<(usize, usize), DimMinimum>, + DefaultAllocator: Allocator> + + Reallocator, C> + + Allocator<(usize, usize), DimMinimum>, { (self.q(), self.r(), self.p) } #[doc(hidden)] - pub fn qrp_internal(&self) -> &MatrixMN { - &self.qrp + pub fn col_piv_qr_internal(&self) -> &MatrixMN { + &self.col_piv_qr } /// Multiplies the provided matrix by the transpose of the `Q` matrix of this decomposition. pub fn q_tr_mul(&self, rhs: &mut Matrix) - // FIXME: do we need a static constraint on the number of rows of rhs? - where S2: StorageMut { + where + S2: StorageMut, + { let dim = self.diag.len(); for i in 0..dim { - let axis = self.qrp.slice_range(i.., i); + let axis = self.col_piv_qr.slice_range(i.., i); let refl = Reflection::new(Unit::new_unchecked(axis), N::zero()); let mut rhs_rows = rhs.rows_range_mut(i..); @@ -187,9 +192,10 @@ where DefaultAllocator: Allocator + Allocator + Allocator> QRP -where DefaultAllocator: Allocator + Allocator + - Allocator<(usize, usize), DimMinimum> +impl> ColPivQR +where + DefaultAllocator: + Allocator + Allocator + Allocator<(usize, usize), DimMinimum>, { /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined. /// @@ -222,20 +228,23 @@ where DefaultAllocator: Allocator + Allocator + ShapeConstraint: SameNumberOfRows, { assert_eq!( - self.qrp.nrows(), + self.col_piv_qr.nrows(), b.nrows(), - "QRP solve matrix dimension mismatch." + "ColPivQR solve matrix dimension mismatch." ); assert!( - self.qrp.is_square(), - "QRP solve: unable to solve a non-square system." + self.col_piv_qr.is_square(), + "ColPivQR solve: unable to solve a non-square system." ); self.q_tr_mul(b); - self.solve_upper_triangular_mut(b) + let solved = self.solve_upper_triangular_mut(b); + self.p.inv_permute_rows(b); + + solved } - // FIXME: duplicate code from the `solve` module. + // TODO: duplicate code from the `solve` module. fn solve_upper_triangular_mut( &self, b: &mut Matrix, @@ -244,7 +253,7 @@ where DefaultAllocator: Allocator + Allocator + S2: StorageMut, ShapeConstraint: SameNumberOfRows, { - let dim = self.qrp.nrows(); + let dim = self.col_piv_qr.nrows(); for k in 0..b.ncols() { let mut b = b.column_mut(k); @@ -263,7 +272,7 @@ where DefaultAllocator: Allocator + Allocator + } b.rows_range_mut(..i) - .axpy(-coeff, &self.qrp.slice_range(..i, i), N::one()); + .axpy(-coeff, &self.col_piv_qr.slice_range(..i, i), N::one()); } } @@ -275,12 +284,12 @@ where DefaultAllocator: Allocator + Allocator + /// Returns `None` if the decomposed matrix is not invertible. pub fn try_inverse(&self) -> Option> { assert!( - self.qrp.is_square(), - "QRP inverse: unable to compute the inverse of a non-square matrix." + self.col_piv_qr.is_square(), + "ColPivQR inverse: unable to compute the inverse of a non-square matrix." ); - // FIXME: is there a less naive method ? - let (nrows, ncols) = self.qrp.data.shape(); + // TODO: is there a less naive method ? + let (nrows, ncols) = self.col_piv_qr.data.shape(); let mut res = MatrixN::identity_generic(nrows, ncols); if self.solve_mut(&mut res) { @@ -293,8 +302,8 @@ where DefaultAllocator: Allocator + Allocator + /// Indicates if the decomposed matrix is invertible. pub fn is_invertible(&self) -> bool { assert!( - self.qrp.is_square(), - "QRP: unable to test the invertibility of a non-square matrix." + self.col_piv_qr.is_square(), + "ColPivQR: unable to test the invertibility of a non-square matrix." ); for i in 0..self.diag.len() { @@ -306,25 +315,32 @@ where DefaultAllocator: Allocator + Allocator + true } - /// Computes the determinant of the decomposed matrix. - pub fn determinant(&self) -> N { - let dim = self.qrp.nrows(); - assert!(self.qrp.is_square(), "QRP determinant: unable to compute the determinant of a non-square matrix."); + /// Computes the determinant of the decomposed matrix. + pub fn determinant(&self) -> N { + let dim = self.col_piv_qr.nrows(); + assert!( + self.col_piv_qr.is_square(), + "ColPivQR determinant: unable to compute the determinant of a non-square matrix." + ); - let mut res = N::one(); - for i in 0 .. dim { - res *= unsafe { *self.diag.vget_unchecked(i) }; - } + let mut res = N::one(); + for i in 0..dim { + res *= unsafe { *self.diag.vget_unchecked(i) }; + } - res * self.p.determinant() - } + res * self.p.determinant() + } } impl, C: Dim, S: Storage> Matrix -where DefaultAllocator: Allocator + Allocator + Allocator> + Allocator<(usize, usize), DimMinimum> +where + DefaultAllocator: Allocator + + Allocator + + Allocator> + + Allocator<(usize, usize), DimMinimum>, { - /// Computes the QRP decomposition of this matrix. - pub fn qrp(self) -> QRP { - QRP::new(self.into_owned()) + /// Computes the QR decomposition (with column pivoting) of this matrix. + pub fn col_piv_qr(self) -> ColPivQR { + ColPivQR::new(self.into_owned()) } } diff --git a/src/linalg/householder.rs b/src/linalg/householder.rs index 3692fe31..ebbffd30 100644 --- a/src/linalg/householder.rs +++ b/src/linalg/householder.rs @@ -34,7 +34,7 @@ pub fn reflection_axis_mut>( if !factor.is_zero() { column.unscale_mut(factor.sqrt()); - (signed_norm, true) + (-signed_norm, true) } else { // TODO: not sure why we don't have a - sign here. (signed_norm, false) diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index c602bfa7..81cf1f38 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -19,6 +19,7 @@ mod inverse; mod lu; mod permutation_sequence; mod qr; +mod col_piv_qr; mod schur; mod solve; mod svd; @@ -39,6 +40,7 @@ pub use self::hessenberg::*; pub use self::lu::*; pub use self::permutation_sequence::*; pub use self::qr::*; +pub use self::col_piv_qr::*; pub use self::schur::*; pub use self::svd::*; pub use self::symmetric_eigen::*; diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs index f404aa5a..fdf6b70a 100644 --- a/src/linalg/qr.rs +++ b/src/linalg/qr.rs @@ -60,8 +60,8 @@ where return QR { qr: matrix, diag }; } - for ite in 0..min_nrows_ncols.value() { - householder::clear_column_unchecked(&mut matrix, &mut diag[ite], ite, 0, None); + for i in 0..min_nrows_ncols.value() { + householder::clear_column_unchecked(&mut matrix, &mut diag[i], i, 0, None); } QR { qr: matrix, diag } diff --git a/tests/linalg/qrp.rs b/tests/linalg/col_piv_qr.rs similarity index 52% rename from tests/linalg/qrp.rs rename to tests/linalg/col_piv_qr.rs index 78b96b74..e9ffba86 100644 --- a/tests/linalg/qrp.rs +++ b/tests/linalg/col_piv_qr.rs @@ -3,19 +3,21 @@ use na::Matrix4; #[test] -fn qrp() { - let m = Matrix4::new ( - 1.0, -1.0, 2.0, 1.0, - -1.0, 3.0, -1.0, -1.0, - 3.0, -5.0, 5.0, 3.0, - 1.0, 2.0, 1.0, -2.0); - let qrp = m.qrp(); - assert!(relative_eq!(qrp.determinant(), 0.0, epsilon = 1.0e-7)); +fn col_piv_qr() { + let m = Matrix4::new( + 1.0, -1.0, 2.0, 1.0, -1.0, 3.0, -1.0, -1.0, 3.0, -5.0, 5.0, 3.0, 1.0, 2.0, 1.0, -2.0, + ); + let col_piv_qr = m.col_piv_qr(); + assert!(relative_eq!( + col_piv_qr.determinant(), + 0.0, + epsilon = 1.0e-7 + )); - let (q, r, p) = qrp.unpack(); + let (q, r, p) = col_piv_qr.unpack(); let mut qr = q * r; - p.inv_permute_columns(& mut qr); + p.inv_permute_columns(&mut qr); assert!(relative_eq!(m, qr, epsilon = 1.0e-7)); } @@ -29,85 +31,89 @@ mod quickcheck_tests { use std::cmp; #[allow(unused_imports)] use crate::core::helper::{RandScalar, RandComplex}; - + quickcheck! { - fn qrp(m: DMatrix<$scalar>) -> bool { + fn col_piv_qr(m: DMatrix<$scalar>) -> bool { let m = m.map(|e| e.0); - let qrp = m.clone().qrp(); - let q = qrp.q(); - let r = qrp.r(); - + let col_piv_qr = m.clone().col_piv_qr(); + let (q, r, p) = col_piv_qr.unpack(); + let mut qr = &q * &r; + p.inv_permute_columns(&mut qr); + println!("m: {}", m); - println!("qrp: {}", &q * &r); - - relative_eq!(m, &q * r, epsilon = 1.0e-7) && - q.is_orthogonal(1.0e-7) - } - - fn qrp_static_5_3(m: Matrix5x3<$scalar>) -> bool { - let m = m.map(|e| e.0); - let qrp = m.qrp(); - let q = qrp.q(); - let r = qrp.r(); - - relative_eq!(m, q * r, epsilon = 1.0e-7) && - q.is_orthogonal(1.0e-7) - } - - fn qrp_static_3_5(m: Matrix3x5<$scalar>) -> bool { - let m = m.map(|e| e.0); - let qrp = m.qrp(); - let q = qrp.q(); - let r = qrp.r(); + println!("col_piv_qr: {}", &q * &r); - relative_eq!(m, q * r, epsilon = 1.0e-7) && - q.is_orthogonal(1.0e-7) - } - - fn qrp_static_square(m: Matrix4<$scalar>) -> bool { - let m = m.map(|e| e.0); - let qrp = m.qrp(); - let q = qrp.q(); - let r = qrp.r(); - - println!("{}{}{}{}", q, r, q * r, m); - - relative_eq!(m, q * r, epsilon = 1.0e-7) && + relative_eq!(m, &qr, epsilon = 1.0e-7) && q.is_orthogonal(1.0e-7) } - - fn qrp_solve(n: usize, nb: usize) -> bool { + + fn col_piv_qr_static_5_3(m: Matrix5x3<$scalar>) -> bool { + let m = m.map(|e| e.0); + let col_piv_qr = m.col_piv_qr(); + let (q, r, p) = col_piv_qr.unpack(); + let mut qr = q * r; + p.inv_permute_columns(&mut qr); + + relative_eq!(m, qr, epsilon = 1.0e-7) && + q.is_orthogonal(1.0e-7) + } + + fn col_piv_qr_static_3_5(m: Matrix3x5<$scalar>) -> bool { + let m = m.map(|e| e.0); + let col_piv_qr = m.col_piv_qr(); + let (q, r, p) = col_piv_qr.unpack(); + let mut qr = q * r; + p.inv_permute_columns(&mut qr); + + relative_eq!(m, qr, epsilon = 1.0e-7) && + q.is_orthogonal(1.0e-7) + } + + fn col_piv_qr_static_square(m: Matrix4<$scalar>) -> bool { + let m = m.map(|e| e.0); + let col_piv_qr = m.col_piv_qr(); + let (q, r, p) = col_piv_qr.unpack(); + let mut qr = q * r; + p.inv_permute_columns(&mut qr); + + println!("{}{}{}{}", q, r, qr, m); + + relative_eq!(m, qr, epsilon = 1.0e-7) && + q.is_orthogonal(1.0e-7) + } + + fn col_piv_qr_solve(n: usize, nb: usize) -> bool { if n != 0 && nb != 0 { let n = cmp::min(n, 50); // To avoid slowing down the test too much. let nb = cmp::min(nb, 50); // To avoid slowing down the test too much. let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0); - - let qrp = m.clone().qrp(); + + let col_piv_qr = m.clone().col_piv_qr(); let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0); let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0); - - if qrp.is_invertible() { - let sol1 = qrp.solve(&b1).unwrap(); - let sol2 = qrp.solve(&b2).unwrap(); - + + if col_piv_qr.is_invertible() { + let sol1 = col_piv_qr.solve(&b1).unwrap(); + let sol2 = col_piv_qr.solve(&b2).unwrap(); + return relative_eq!(&m * sol1, b1, epsilon = 1.0e-6) && relative_eq!(&m * sol2, b2, epsilon = 1.0e-6) } } - + return true; } - - fn qrp_solve_static(m: Matrix4<$scalar>) -> bool { + + fn col_piv_qr_solve_static(m: Matrix4<$scalar>) -> bool { let m = m.map(|e| e.0); - let qrp = m.qrp(); + let col_piv_qr = m.col_piv_qr(); let b1 = Vector4::<$scalar>::new_random().map(|e| e.0); let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0); - - if qrp.is_invertible() { - let sol1 = qrp.solve(&b1).unwrap(); - let sol2 = qrp.solve(&b2).unwrap(); - + + if col_piv_qr.is_invertible() { + let sol1 = col_piv_qr.solve(&b1).unwrap(); + let sol2 = col_piv_qr.solve(&b2).unwrap(); + relative_eq!(m * sol1, b1, epsilon = 1.0e-6) && relative_eq!(m * sol2, b2, epsilon = 1.0e-6) } @@ -116,29 +122,29 @@ mod quickcheck_tests { } } - fn qrp_inverse(n: usize) -> bool { + fn col_piv_qr_inverse(n: usize) -> bool { let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much. let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0); - - if let Some(m1) = m.clone().qrp().try_inverse() { + + if let Some(m1) = m.clone().col_piv_qr().try_inverse() { let id1 = &m * &m1; let id2 = &m1 * &m; - + id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5) } else { true } } - - fn qrp_inverse_static(m: Matrix4<$scalar>) -> bool { + + fn col_piv_qr_inverse_static(m: Matrix4<$scalar>) -> bool { let m = m.map(|e| e.0); - let qrp = m.qrp(); - - if let Some(m1) = qrp.try_inverse() { + let col_piv_qr = m.col_piv_qr(); + + if let Some(m1) = col_piv_qr.try_inverse() { let id1 = &m * &m1; let id2 = &m1 * &m; - + id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5) } else { @@ -153,4 +159,3 @@ mod quickcheck_tests { gen_tests!(complex, RandComplex); gen_tests!(f64, RandScalar); } - diff --git a/tests/linalg/mod.rs b/tests/linalg/mod.rs index 7fc01396..66a0c06a 100644 --- a/tests/linalg/mod.rs +++ b/tests/linalg/mod.rs @@ -1,6 +1,7 @@ mod balancing; mod bidiagonal; mod cholesky; +mod col_piv_qr; mod convolution; mod eigen; mod exp;