From 87eb03063f38b0e18cbca766205d407a6b2d4348 Mon Sep 17 00:00:00 2001 From: prismatica Date: Sat, 25 Mar 2017 10:57:45 -0400 Subject: [PATCH 1/6] Add rectangular matrix support for FullPivLu --- src/matrix/base/mod.rs | 8 +- src/matrix/decomposition/lu.rs | 278 ++++++++++++++++++++++++++++++--- src/vector/impl_vec.rs | 14 ++ 3 files changed, 281 insertions(+), 19 deletions(-) diff --git a/src/matrix/base/mod.rs b/src/matrix/base/mod.rs index 591ecae..c784412 100644 --- a/src/matrix/base/mod.rs +++ b/src/matrix/base/mod.rs @@ -51,11 +51,17 @@ pub trait BaseMatrix: Sized { /// Row stride in the matrix. fn row_stride(&self) -> usize; - /// Returns true if the matrix contais no elements + /// Returns true if the matrix contains no elements fn is_empty(&self) -> bool { self.rows() == 0 || self.cols() == 0 } + /// Returns true if the matrix has the same number of + /// rows as columns. + fn is_square(&self) -> bool { + self.rows() == self.cols() + } + /// Top left index of the matrix. fn as_ptr(&self) -> *const T; diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 7913935..750d8d3 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -5,10 +5,14 @@ use vector::Vector; use error::{Error, ErrorKind}; use std::any::Any; +use std::ops::Mul; use std::cmp; use libnum::{Float, Zero, One}; +use utils; +use norm::{Euclidean, VectorNorm}; + use matrix::decomposition::Decomposition; /// Result of unpacking an instance of @@ -383,10 +387,6 @@ impl FullPivLu { /// Performs the decomposition. pub fn decompose(matrix: Matrix) -> Result { - assert!( - matrix.rows() == matrix.cols(), - "Matrix must be square for LU decomposition."); - let mut lu = matrix; let nrows = lu.rows(); @@ -469,8 +469,43 @@ impl FullPivLu where T: Any + Float { assert!(b.size() == self.lu.rows(), "Right-hand side vector must have compatible size."); - let b = lu_forward_substitution(&self.lu, &self.p * b); - back_substitution(&self.lu, b).map(|x| &self.q * x) + // These two products are guarenteed to exist, since self.p + // is a permutation matrix, and the lower-triangular portion + // of self.lu is a composition of elementary row operations. + let x = lu_forward_substitution(&self.lu, &self.p * &b); + + // Here, though, the result might not exist. This happens + // if the columns of the matrix don't span the entire range + // (i.e. the matrix is rank-deficient or the matrix has more + // rows than columns). + // + // We get around this by only considering the top-left square + // invertible portion of U, solving that problem, and verifying + // the solution later. + let r = self.rank(); + let nonzero_u = self.lu.sub_slice([0, 0], r, r); + + let x = back_substitution(&nonzero_u, x.top(r)).unwrap(); + + // Now we pad x back up to its actual size and map through + // the last permutation matrix. + let x = Vector::from_fn( + self.lu.cols(), + |i| if i < r { unsafe{ *x.get_unchecked(i) }} else { T::zero() }); + + let x = &self.q * x; + + // We have a potential solution, but we need to verify that it actually + // solves the equation. + let norm = VectorNorm::norm(&Euclidean, &x); + let diff_norm = VectorNorm::norm(&Euclidean, &(self * &x - b)); + + if diff_norm < norm * self.epsilon() { + Ok(x) + } else { + Err(Error::new(ErrorKind::AlgebraFailure, + "No solution exists")) + } } /// Computes the inverse of the matrix which this LUP decomposition @@ -484,6 +519,10 @@ impl FullPivLu where T: Any + Float { let mut inv = Matrix::zeros(n, n); let mut e = Vector::zeros(n); + assert!( + self.lu.is_square(), + "Matrix must be square to take inverses."); + if !self.is_invertible() { return Err( Error::new( @@ -513,6 +552,10 @@ impl FullPivLu where T: Any + Float { /// # Panics /// If the underlying matrix is non-square. pub fn det(&self) -> T { + assert!( + self.lu.is_square(), + "Matrix must be square to take determinants."); + // Recall that the determinant of a triangular matrix // is the product of its diagonal entries. Also, // the determinant of L is implicitly 1. @@ -558,8 +601,7 @@ impl FullPivLu where T: Any + Float { /// Returns whether the matrix is invertible. /// - /// Empty matrices are considered to be invertible for - /// the sake of this function. + /// Empty matrices are invertible while rectangular matrices are not. /// /// # Examples /// @@ -580,16 +622,17 @@ impl FullPivLu where T: Any + Float { /// # } /// ``` pub fn is_invertible(&self) -> bool { - let diag_size = cmp::min(self.lu.rows(), self.lu.cols()); + if !self.lu.is_square() { + false - if diag_size > 0 { - let diag_last = diag_size - 1; - let last = - unsafe { self.lu.get_unchecked([diag_last, diag_last]) }; + } else if self.lu.is_empty() { + true - last.abs() > self.epsilon() } else { - true + let diag_last = self.lu.rows() - 1; + let last = self.lu.get([diag_last, diag_last]).unwrap(); + + last.abs() > self.epsilon() } } @@ -598,6 +641,79 @@ impl FullPivLu where T: Any + Float { } } +/// Multiplies an LU decomposed matrix by vector. +impl Mul> for FullPivLu where T: Float { + type Output = Vector; + + fn mul(self, m: Vector) -> Vector { + (&self) * (&m) + } +} + +/// Multiplies an LU decomposed matrix by vector. +impl<'a, T> Mul> for &'a FullPivLu where T: Float { + type Output = Vector; + + fn mul(self, m: Vector) -> Vector { + self * &m + } +} + +/// Multiplies an LU decomposed matrix by vector. +impl<'a, T> Mul<&'a Vector> for FullPivLu where T: Float { + type Output = Vector; + + fn mul(self, m: &Vector) -> Vector { + (&self) * m + } +} + +/// Multiplies an LU decomposed matrix by vector. +impl<'a, 'b, T> Mul<&'b Vector> for &'a FullPivLu where T: Float { + type Output = Vector; + + fn mul(self, v: &Vector) -> Vector { + assert!( + v.size() == self.lu.cols(), + "Matrix and Vector dimensions do not agree."); + + let r = self.lu.rows(); + let c = self.lu.cols(); + + let v = self.q.inverse() * v; + + let mut new_data = Vec::with_capacity(r); + + // First, do the U multiplication. + for i in 0..r { + let entry = + utils::dot( + &self.lu.data()[(i*(c+1)) .. ((i+1)*c)], + &v.data()[i..c]); + + new_data.push(entry); + } + + // Now, do L multiplication + for i in (1..r).rev() { + let row_begin = i * c; + + let entry = + utils::dot( + &self.lu.data()[row_begin .. row_begin + i], + &new_data[0..i]); + + new_data[i] = new_data[i] + entry; + } + + let mut v = Vector::new(new_data); + + self.p.inverse().permute_vector_in_place(&mut v); + + v + } +} + /// Performs Gaussian elimination in the lower-right hand corner starting at /// (index, index). fn gaussian_elimination(lu: &mut Matrix, index: usize) { @@ -622,7 +738,6 @@ fn gaussian_elimination(lu: &mut Matrix, index: usize) { /// /// This is equivalent to solving the system Lx = b. fn lu_forward_substitution(lu: &Matrix, b: Vector) -> Vector { - assert!(lu.rows() == lu.cols(), "LU matrix must be square."); assert!(b.size() == lu.rows(), "LU matrix and RHS vector must be compatible."); let mut x = b; @@ -935,14 +1050,125 @@ mod tests { } #[test] - #[should_panic] fn full_piv_lu_decompose_rectangular() { let x = matrix![ -3.0, 0.0, 4.0; -12.0, 5.0, 17.0; 15.0, 0.0, -18.0; -6.0, 0.0, 20.0]; - FullPivLu::decompose(x.clone()).unwrap(); + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + assert_eq!(lu.rank(), 3); + assert!(!lu.is_invertible()); + + let LUPQ { l, u, p, q } = lu.unpack(); + + let y = p.inverse() * &l * &u * q.inverse(); + + assert_matrix_eq!(x, y, comp = float); + assert!(is_lower_triangular(&l)); + assert!(is_upper_triangular(&u)); + } + + #[test] + fn full_piv_lu_decompose_rectangular2() { + let x = matrix![ -3.0, 0.0, 4.0; + -6.0, 0.0, 20.0]; + + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + assert_eq!(lu.rank(), 2); + assert!(!lu.is_invertible()); + + let LUPQ { l, u, p, q } = lu.unpack(); + + let y = p.inverse() * &l * &u * q.inverse(); + + assert_matrix_eq!(x, y, comp = float); + assert!(is_lower_triangular(&l)); + assert!(is_upper_triangular(&u)); + } + + #[test] + #[should_panic] + fn full_piv_lu_rectangular_det() { + let x = matrix![ -3.0, 0.0, 4.0; + -6.0, 0.0, 20.0]; + + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + lu.det(); + } + + #[test] + #[should_panic] + fn full_piv_lu_rectangular_inverse() { + let x = matrix![ -3.0, 0.0, 4.0; + -6.0, 0.0, 20.0]; + + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + lu.inverse().unwrap(); + } + + #[test] + pub fn full_piv_lu_solve_rectangular_matrix() { + let a = matrix![ 5.0, 0.0, 0.0, 1.0; + 2.0, 2.0, 2.0, 1.0; + 1.0, 6.0, 4.0, 5.0 ]; + let b = vector![9.0, 16.0, 45.0]; + + let lu = FullPivLu::decompose(a.clone()).unwrap(); + let x = lu.solve(b.clone()).unwrap(); + + // Note that we shouldn't just choose an expected + // value here, since there are an infinte number of + // solutions + let res = a * x; + + assert_vector_eq!(res, b, comp = ulp, tol = 100); + } + + #[test] + pub fn full_piv_lu_solve_rectangular_matrix2() { + let a = matrix![ 5.0, 1.0; + 2.0, 1.0; + 1.0, 5.0 ]; + let b = vector![7.0, 4.0, 11.0]; + + let lu = FullPivLu::decompose(a.clone()).unwrap(); + let x = lu.solve(b.clone()).unwrap(); + + let res = a * x; + + assert_vector_eq!(res, b, comp = ulp, tol = 100); + } + + #[test] + pub fn full_piv_lu_solve_rectangular_matrix3() { + let a = matrix![ 0.005, 0.001; + 0.002, 0.001; + 0.001, 0.005 ]; + let b = vector![0.0007, 0.0004, 0.0011]; + + let lu = FullPivLu::decompose(a.clone()).unwrap(); + let x = lu.solve(b.clone()).unwrap(); + + let res = a * x; + + assert_vector_eq!(res, b, comp = ulp, tol = 100); + } + + #[test] + pub fn full_piv_lu_solve_no_solution() { + let a = matrix![ 1.0, 0.0; + 0.0, 1.0; + 0.0, 0.0]; + let b = vector![0.0, 0.0, 1.0]; + + let lu = FullPivLu::decompose(a).unwrap(); + + assert!(lu.solve(b).is_err()) } #[test] @@ -991,6 +1217,22 @@ mod tests { assert!(lu.inverse().is_err()); } + #[test] + pub fn full_piv_lu_mul() { + let x = matrix![5.0, 0.0; + 4.0, 5.0; + 9.0, 5.0]; + + let lu = FullPivLu::decompose(x).unwrap(); + + let x = vector![1.0, 2.0]; + let expected = vector![5.0, 14.0, 19.0]; + + let y = &lu * &x; + + assert_vector_eq!(y, expected, comp = ulp, tol = 100); + } + #[test] pub fn full_piv_lu_empty_matrix() { use matrix::base::BaseMatrix; diff --git a/src/vector/impl_vec.rs b/src/vector/impl_vec.rs index 8cba075..9022d01 100644 --- a/src/vector/impl_vec.rs +++ b/src/vector/impl_vec.rs @@ -67,6 +67,20 @@ impl Vector { &self.data } + /// Returns the top `count` rows of the Vector. + pub fn top(&self, count: usize) -> Vector where T: Clone { + assert!( + count <= self.size, + "count: {}, self.size: {}", + count, + self.size); + + Vector { + size: count, + data: self.data[0..count].to_vec() + } + } + /// Returns a mutable slice of the underlying data. pub fn mut_data(&mut self) -> &mut [T] { &mut self.data From f7dd97f691608037e03211b698b13c6741e68e29 Mon Sep 17 00:00:00 2001 From: prismatica Date: Sat, 22 Apr 2017 17:11:59 -0400 Subject: [PATCH 2/6] Optimize memory for tall-matrix case for FullPivLu Previously, the L matrix would always be a mxm matrix, and the U matrix would be a mxn matrix. This is because the unpack method used the self.lu matrix exclusively to store the U matrix. In order to avoid wasting large amounts of memory, the internal storage is repurposed to store either the L or U matrix, whichever can fit into the mxn size. The other matrix is allocated as a square matrix. Also removed the Mul trait for FullPivLu, and moved that trait to a private internal function. --- src/matrix/decomposition/lu.rs | 186 +++++++++++++++++++++------------ 1 file changed, 120 insertions(+), 66 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 750d8d3..89e3eac 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -5,7 +5,6 @@ use vector::Vector; use error::{Error, ErrorKind}; use std::any::Any; -use std::ops::Mul; use std::cmp; use libnum::{Float, Zero, One}; @@ -13,6 +12,7 @@ use libnum::{Float, Zero, One}; use utils; use norm::{Euclidean, VectorNorm}; +use matrix::DiagOffset; use matrix::decomposition::Decomposition; /// Result of unpacking an instance of @@ -351,9 +351,37 @@ impl Decomposition for FullPivLu { fn unpack(self) -> LUPQ { use internal_utils::nullify_lower_triangular_part; - let l = unit_lower_triangular_part(&self.lu); - let mut u = self.lu; - nullify_lower_triangular_part(&mut u); + use internal_utils::nullify_upper_triangular_part; + + // We want the L matrix to be M x P, and the U matrix to be + // P x N, where P is the minimum of M and N. Since self.lu + // is already M x N, and we want to reuse that memory to store + // at least one of the matrices, we need to figure out which + // matrix to store in self.lu, and which to allocate anew. The + // nice thing here is that the new allocation will always be + // PxP, which is the smaller of the two matrices. + + let (l,u) = + if self.lu.rows() <= self.lu.cols() { + let l = unit_lower_triangular_part(&self.lu); + + let mut u = self.lu; + nullify_lower_triangular_part(&mut u); + + (l,u) + + } else { + let u = upper_triangular_part(&self.lu); + + let mut l = self.lu; + nullify_upper_triangular_part(&mut l); + + for val in l.diag_iter_mut(DiagOffset::Main) { + *val = T::one(); + } + + (l,u) + }; LUPQ { l: l, @@ -498,7 +526,10 @@ impl FullPivLu where T: Any + Float { // We have a potential solution, but we need to verify that it actually // solves the equation. let norm = VectorNorm::norm(&Euclidean, &x); - let diff_norm = VectorNorm::norm(&Euclidean, &(self * &x - b)); + + let rhs = self.map_vector(&x) - b; + + let diff_norm = VectorNorm::norm(&Euclidean, &rhs); if diff_norm < norm * self.epsilon() { Ok(x) @@ -639,53 +670,22 @@ impl FullPivLu where T: Any + Float { fn epsilon(&self) -> T { self.lu.get([0, 0]).unwrap_or(&T::one()).abs() * T::epsilon() } -} -/// Multiplies an LU decomposed matrix by vector. -impl Mul> for FullPivLu where T: Float { - type Output = Vector; - - fn mul(self, m: Vector) -> Vector { - (&self) * (&m) - } -} - -/// Multiplies an LU decomposed matrix by vector. -impl<'a, T> Mul> for &'a FullPivLu where T: Float { - type Output = Vector; - - fn mul(self, m: Vector) -> Vector { - self * &m - } -} - -/// Multiplies an LU decomposed matrix by vector. -impl<'a, T> Mul<&'a Vector> for FullPivLu where T: Float { - type Output = Vector; - - fn mul(self, m: &Vector) -> Vector { - (&self) * m - } -} - -/// Multiplies an LU decomposed matrix by vector. -impl<'a, 'b, T> Mul<&'b Vector> for &'a FullPivLu where T: Float { - type Output = Vector; - - fn mul(self, v: &Vector) -> Vector { + fn map_vector(&self, v: &Vector) -> Vector { assert!( v.size() == self.lu.cols(), "Matrix and Vector dimensions do not agree."); let r = self.lu.rows(); let c = self.lu.cols(); + let d = cmp::min(r, c); let v = self.q.inverse() * v; let mut new_data = Vec::with_capacity(r); // First, do the U multiplication. - for i in 0..r { + for i in 0..d { let entry = utils::dot( &self.lu.data()[(i*(c+1)) .. ((i+1)*c)], @@ -694,14 +694,22 @@ impl<'a, 'b, T> Mul<&'b Vector> for &'a FullPivLu where T: Float { new_data.push(entry); } + for _ in d..r { + new_data.push(T::zero()); + } + // Now, do L multiplication for i in (1..r).rev() { let row_begin = i * c; - let entry = - utils::dot( - &self.lu.data()[row_begin .. row_begin + i], - &new_data[0..i]); + let entry = { + let offset = cmp::min(d, i); + + let data_slice = &new_data[0..offset]; + let lu_slice = &self.lu.data()[row_begin .. row_begin + offset]; + + utils::dot(lu_slice, data_slice) + }; new_data[i] = new_data[i] + entry; } @@ -741,7 +749,9 @@ fn lu_forward_substitution(lu: &Matrix, b: Vector) -> Vector assert!(b.size() == lu.rows(), "LU matrix and RHS vector must be compatible."); let mut x = b; - for (i, row) in lu.row_iter().enumerate().skip(1) { + let diag_size = cmp::min(lu.rows(), lu.cols()); + + for (i, row) in lu.row_iter().enumerate().take(diag_size).skip(1) { // Note that at time of writing we need raw_slice here for // auto-vectorization to kick in let adjustment = row.raw_slice() @@ -777,6 +787,24 @@ fn unit_lower_triangular_part(matrix: &M) -> Matrix Matrix::new(m, m, data) } +fn upper_triangular_part(matrix: &M) -> Matrix + where T: Zero + One + Clone, M: BaseMatrix { + + let n = matrix.cols(); + let mut data = Vec::::with_capacity(n * n); + + for (i, row) in matrix.row_iter().take(n).enumerate() { + for _ in 0..i { + data.push(T::zero()); + } + + for element in row.iter().skip(i).cloned() { + data.push(element); + } + } + + Matrix::new(n, n, data) +} impl Matrix where T: Any + Float { @@ -1130,11 +1158,53 @@ mod tests { } #[test] - pub fn full_piv_lu_solve_rectangular_matrix2() { - let a = matrix![ 5.0, 1.0; - 2.0, 1.0; - 1.0, 5.0 ]; - let b = vector![7.0, 4.0, 11.0]; + pub fn full_piv_lu_solve_rectangular_matrix_tall() { + let a = matrix![ 1.0, 1.0; + 2.0, 1.0; + 3.0, 1.0; + 4.0, 1.0; + 5.0, 1.0; + 6.0, 1.0; + 7.0, 1.0; + 8.0, 1.0; + 9.0, 1.0; + 10.0, 1.0; + 11.0, 1.0; + 12.0, 1.0; + 13.0, 1.0; + 14.0, 1.0; + 15.0, 1.0; + 16.0, 1.0 ]; + let b = vector![ 2.0, + 3.0, + 4.0, + 5.0, + 6.0, + 7.0, + 8.0, + 9.0, + 10.0, + 11.0, + 12.0, + 13.0, + 14.0, + 15.0, + 16.0, + 17.0]; + + let lu = FullPivLu::decompose(a.clone()).unwrap(); + let x = lu.solve(b.clone()).unwrap(); + + let res = a * x; + + assert_vector_eq!(res, b, comp = ulp, tol = 100); + } + + #[test] + pub fn full_piv_lu_solve_rectangular_matrix_wide() { + let a = matrix![ 0.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0; + 1.0, 0.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]; + let b = vector![2.0, 3.0]; let lu = FullPivLu::decompose(a.clone()).unwrap(); let x = lu.solve(b.clone()).unwrap(); @@ -1145,7 +1215,7 @@ mod tests { } #[test] - pub fn full_piv_lu_solve_rectangular_matrix3() { + pub fn full_piv_lu_solve_rectangular_small_values() { let a = matrix![ 0.005, 0.001; 0.002, 0.001; 0.001, 0.005 ]; @@ -1217,22 +1287,6 @@ mod tests { assert!(lu.inverse().is_err()); } - #[test] - pub fn full_piv_lu_mul() { - let x = matrix![5.0, 0.0; - 4.0, 5.0; - 9.0, 5.0]; - - let lu = FullPivLu::decompose(x).unwrap(); - - let x = vector![1.0, 2.0]; - let expected = vector![5.0, 14.0, 19.0]; - - let y = &lu * &x; - - assert_vector_eq!(y, expected, comp = ulp, tol = 100); - } - #[test] pub fn full_piv_lu_empty_matrix() { use matrix::base::BaseMatrix; From 24f763c9a465590a3c43e8d1b586221feca4ad1e Mon Sep 17 00:00:00 2001 From: prismatica Date: Mon, 1 May 2017 18:33:27 -0400 Subject: [PATCH 3/6] Update verification steps to algorithm discussed in pull request. --- src/matrix/decomposition/lu.rs | 136 ++++++++++++--------------------- 1 file changed, 49 insertions(+), 87 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 89e3eac..869aa2a 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -6,15 +6,18 @@ use error::{Error, ErrorKind}; use std::any::Any; use std::cmp; +use std::fmt::Debug; use libnum::{Float, Zero, One}; use utils; -use norm::{Euclidean, VectorNorm}; use matrix::DiagOffset; use matrix::decomposition::Decomposition; +use internal_utils::nullify_lower_triangular_part; +use internal_utils::nullify_upper_triangular_part; + /// Result of unpacking an instance of /// [PartialPivLu](struct.PartialPivLu.html). #[derive(Debug, Clone)] @@ -140,7 +143,6 @@ impl Decomposition for PartialPivLu { type Factors = LUP; fn unpack(self) -> LUP { - use internal_utils::nullify_lower_triangular_part; let l = unit_lower_triangular_part(&self.lu); let mut u = self.lu; nullify_lower_triangular_part(&mut u); @@ -350,8 +352,6 @@ impl Decomposition for FullPivLu { type Factors = LUPQ; fn unpack(self) -> LUPQ { - use internal_utils::nullify_lower_triangular_part; - use internal_utils::nullify_upper_triangular_part; // We want the L matrix to be M x P, and the U matrix to be // P x N, where P is the minimum of M and N. Since self.lu @@ -460,7 +460,7 @@ impl FullPivLu { // TODO: Remove Any bound (cannot for the time being, since // back substitution uses Any bound) -impl FullPivLu where T: Any + Float { +impl FullPivLu where T: Any + Float + Debug { /// Solves the linear system `Ax = b`. /// @@ -497,19 +497,47 @@ impl FullPivLu where T: Any + Float { assert!(b.size() == self.lu.rows(), "Right-hand side vector must have compatible size."); - // These two products are guarenteed to exist, since self.p - // is a permutation matrix, and the lower-triangular portion - // of self.lu is a composition of elementary row operations. - let x = lu_forward_substitution(&self.lu, &self.p * &b); + let diag_size = cmp::min(self.lu.rows(), self.lu.cols()); - // Here, though, the result might not exist. This happens - // if the columns of the matrix don't span the entire range - // (i.e. the matrix is rank-deficient or the matrix has more - // rows than columns). - // - // We get around this by only considering the top-left square - // invertible portion of U, solving that problem, and verifying - // the solution later. + let b = &self.p * &b; + + let x = + lu_forward_substitution( + &self.lu.sub_slice([0, 0], diag_size, diag_size), + b.top(diag_size)); + + // Check that the bottom elements of X are all zero. + let eps = self.epsilon(); + + for v in x.iter().skip(diag_size) { + if v.abs() > eps { + return Err( + Error::new( + ErrorKind::AlgebraFailure, + "No solution exists")); + } + } + + // Get the lower portion of the L matrix and verify the solution + // found. + for (i, expected) in b.iter().enumerate().skip(diag_size) { + let row_begin = self.lu.row_stride() * i; + let row_end = row_begin + diag_size; + + let val = + utils::dot( + &self.lu.data()[row_begin .. row_end], + &x.data()[0 .. diag_size]); + + if (val - *expected).abs() > eps { + return Err( + Error::new( + ErrorKind::AlgebraFailure, + "No solution exists")); + } + } + + // Consider the top-left square invertible portion of U. let r = self.rank(); let nonzero_u = self.lu.sub_slice([0, 0], r, r); @@ -521,22 +549,7 @@ impl FullPivLu where T: Any + Float { self.lu.cols(), |i| if i < r { unsafe{ *x.get_unchecked(i) }} else { T::zero() }); - let x = &self.q * x; - - // We have a potential solution, but we need to verify that it actually - // solves the equation. - let norm = VectorNorm::norm(&Euclidean, &x); - - let rhs = self.map_vector(&x) - b; - - let diff_norm = VectorNorm::norm(&Euclidean, &rhs); - - if diff_norm < norm * self.epsilon() { - Ok(x) - } else { - Err(Error::new(ErrorKind::AlgebraFailure, - "No solution exists")) - } + Ok(&self.q * x) } /// Computes the inverse of the matrix which this LUP decomposition @@ -670,56 +683,6 @@ impl FullPivLu where T: Any + Float { fn epsilon(&self) -> T { self.lu.get([0, 0]).unwrap_or(&T::one()).abs() * T::epsilon() } - - fn map_vector(&self, v: &Vector) -> Vector { - assert!( - v.size() == self.lu.cols(), - "Matrix and Vector dimensions do not agree."); - - let r = self.lu.rows(); - let c = self.lu.cols(); - let d = cmp::min(r, c); - - let v = self.q.inverse() * v; - - let mut new_data = Vec::with_capacity(r); - - // First, do the U multiplication. - for i in 0..d { - let entry = - utils::dot( - &self.lu.data()[(i*(c+1)) .. ((i+1)*c)], - &v.data()[i..c]); - - new_data.push(entry); - } - - for _ in d..r { - new_data.push(T::zero()); - } - - // Now, do L multiplication - for i in (1..r).rev() { - let row_begin = i * c; - - let entry = { - let offset = cmp::min(d, i); - - let data_slice = &new_data[0..offset]; - let lu_slice = &self.lu.data()[row_begin .. row_begin + offset]; - - utils::dot(lu_slice, data_slice) - }; - - new_data[i] = new_data[i] + entry; - } - - let mut v = Vector::new(new_data); - - self.p.inverse().permute_vector_in_place(&mut v); - - v - } } /// Performs Gaussian elimination in the lower-right hand corner starting at @@ -745,13 +708,12 @@ fn gaussian_elimination(lu: &mut Matrix, index: usize) { /// to the strictly lower triangular part of L. /// /// This is equivalent to solving the system Lx = b. -fn lu_forward_substitution(lu: &Matrix, b: Vector) -> Vector { +fn lu_forward_substitution(lu: &M, b: Vector) -> Vector + where T: Float, M: BaseMatrix { assert!(b.size() == lu.rows(), "LU matrix and RHS vector must be compatible."); let mut x = b; - let diag_size = cmp::min(lu.rows(), lu.cols()); - - for (i, row) in lu.row_iter().enumerate().take(diag_size).skip(1) { + for (i, row) in lu.row_iter().enumerate().skip(1) { // Note that at time of writing we need raw_slice here for // auto-vectorization to kick in let adjustment = row.raw_slice() From 17c6409f0a3ec4d3eda6f0d11c8e5647a2d84b73 Mon Sep 17 00:00:00 2001 From: prismatica Date: Mon, 1 May 2017 18:36:06 -0400 Subject: [PATCH 4/6] Remove stray Debug bounds --- src/matrix/decomposition/lu.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 869aa2a..9e2e0c2 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -6,7 +6,6 @@ use error::{Error, ErrorKind}; use std::any::Any; use std::cmp; -use std::fmt::Debug; use libnum::{Float, Zero, One}; @@ -460,7 +459,7 @@ impl FullPivLu { // TODO: Remove Any bound (cannot for the time being, since // back substitution uses Any bound) -impl FullPivLu where T: Any + Float + Debug { +impl FullPivLu where T: Any + Float { /// Solves the linear system `Ax = b`. /// From ecebaf4fd81206458741eda5801d7c6410bdf32b Mon Sep 17 00:00:00 2001 From: prismatica Date: Mon, 1 May 2017 19:26:23 -0400 Subject: [PATCH 5/6] Add more tests for rectangular solve. --- src/matrix/decomposition/lu.rs | 50 ++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 9e2e0c2..d565c9f 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -1202,6 +1202,56 @@ mod tests { assert!(lu.solve(b).is_err()) } + #[test] + pub fn full_piv_lu_solve_no_solution_tall_matrix() { + let a = matrix![ 1.0, 0.0; + 2.0, 1.0; + 3.0, 0.0; + 4.0, 1.0; + 5.0, 0.0; + 6.0, 1.0; + 7.0, 0.0; + 8.0, 1.0; + 9.0, 0.0; + 10.0, 1.0; + 11.0, 0.0; + 12.0, 1.0; + 13.0, 0.0; + 14.0, 1.0; + 15.0, 0.0; + 16.0, 1.0 ]; + let b = vector![ 2.0, + 3.0, + 4.0, + 5.0, + 6.0, + 7.0, + 8.0, + 9.0, + 10.0, + 11.0, + 12.0, + 13.0, + 14.0, + 15.0, + 16.0, + 17.0]; + + let lu = FullPivLu::decompose(a.clone()).unwrap(); + assert!(lu.solve(b).is_err()); + } + + #[test] + pub fn full_piv_lu_solve_no_solution_matrix_wide() { + let a = matrix![ 1.0, 2.0, 3.0; + 1.0, 2.0, 3.0]; + let b = vector![2.0, 3.0]; + + let lu = FullPivLu::decompose(a.clone()).unwrap(); + + assert!(lu.solve(b).is_err()); + } + #[test] pub fn full_piv_lu_solve_arbitrary_matrix() { let x = matrix![ 5.0, 0.0, 0.0, 1.0; From ee6f44bc522b1b1488e4b7dd28344226d9c15d34 Mon Sep 17 00:00:00 2001 From: prismatica Date: Thu, 4 May 2017 19:24:11 -0400 Subject: [PATCH 6/6] Change to using rank instead of diag_size in solve This fixes a dumb bug in the previous commit. Instead of using diag_size, the verification step should use the rank. This applies to the full-pivot rectangular solve routine. --- src/matrix/decomposition/lu.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index d565c9f..c4e8703 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -497,6 +497,7 @@ impl FullPivLu where T: Any + Float { "Right-hand side vector must have compatible size."); let diag_size = cmp::min(self.lu.rows(), self.lu.cols()); + let rank = self.rank(); let b = &self.p * &b; @@ -508,7 +509,7 @@ impl FullPivLu where T: Any + Float { // Check that the bottom elements of X are all zero. let eps = self.epsilon(); - for v in x.iter().skip(diag_size) { + for v in x.iter().skip(rank) { if v.abs() > eps { return Err( Error::new( @@ -519,7 +520,7 @@ impl FullPivLu where T: Any + Float { // Get the lower portion of the L matrix and verify the solution // found. - for (i, expected) in b.iter().enumerate().skip(diag_size) { + for (i, expected) in b.iter().enumerate().skip(rank) { let row_begin = self.lu.row_stride() * i; let row_end = row_begin + diag_size; @@ -537,16 +538,15 @@ impl FullPivLu where T: Any + Float { } // Consider the top-left square invertible portion of U. - let r = self.rank(); - let nonzero_u = self.lu.sub_slice([0, 0], r, r); + let nonzero_u = self.lu.sub_slice([0, 0], rank, rank); - let x = back_substitution(&nonzero_u, x.top(r)).unwrap(); + let x = back_substitution(&nonzero_u, x.top(rank)).unwrap(); // Now we pad x back up to its actual size and map through // the last permutation matrix. let x = Vector::from_fn( self.lu.cols(), - |i| if i < r { unsafe{ *x.get_unchecked(i) }} else { T::zero() }); + |i| if i < rank { unsafe{ *x.get_unchecked(i) }} else { T::zero() }); Ok(&self.q * x) }