From 59c268ff31ec5f4b5c8b472d1f62760c1abf7ea0 Mon Sep 17 00:00:00 2001 From: prismatica Date: Sat, 25 Feb 2017 10:23:56 -0500 Subject: [PATCH 01/12] Added a "get" method to BaseMatrix. This resolves issue #151. --- src/matrix/base/mod.rs | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/matrix/base/mod.rs b/src/matrix/base/mod.rs index 6c19fb8..321401d 100644 --- a/src/matrix/base/mod.rs +++ b/src/matrix/base/mod.rs @@ -80,6 +80,35 @@ pub trait BaseMatrix: Sized { &*(self.as_ptr().offset((index[0] * self.row_stride() + index[1]) as isize)) } + /// Get a reference to a point in the matrix. + /// + /// # Examples + /// + /// ``` + /// # #[macro_use] extern crate rulinalg; fn main() { + /// use rulinalg::matrix::{Matrix, BaseMatrix}; + /// + /// let mat = matrix![0, 1; + /// 3, 4; + /// 6, 7]; + /// + /// assert_eq!(mat.get([0, 2]), None); + /// assert_eq!(mat.get([3, 0]), None); + /// + /// assert_eq!( *mat.get([0, 0]).unwrap(), 0) + /// # } + /// ``` + fn get(&self, index: [usize; 2]) -> Option<&T> { + let row_ind = index[0]; + let col_ind = index[1]; + + if row_ind >= self.rows() || col_ind >= self.cols() { + None + } else { + unsafe { Some(self.get_unchecked(index)) } + } + } + /// Returns the column of a matrix at the given index. /// `None` if the index is out of bounds. /// From 4244708ba972d658df43ce94d5dc8e7d4b11942c Mon Sep 17 00:00:00 2001 From: prismatica Date: Sat, 25 Feb 2017 10:41:06 -0500 Subject: [PATCH 02/12] Added a "get_mut" function to BaseMatrixMut. This is related to issue #151. --- src/matrix/base/mod.rs | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/src/matrix/base/mod.rs b/src/matrix/base/mod.rs index 321401d..79930d3 100644 --- a/src/matrix/base/mod.rs +++ b/src/matrix/base/mod.rs @@ -1173,6 +1173,37 @@ pub trait BaseMatrixMut: BaseMatrix { &mut *(self.as_mut_ptr().offset((index[0] * self.row_stride() + index[1]) as isize)) } + /// Get a mutable reference to a point in the matrix. + /// + /// # Examples + /// + /// ``` + /// # #[macro_use] extern crate rulinalg; fn main() { + /// use rulinalg::matrix::{Matrix, BaseMatrix, BaseMatrixMut}; + /// + /// let mut mat = matrix![0, 1; + /// 3, 4; + /// 6, 7]; + /// + /// assert_eq!(mat.get_mut([0, 2]), None); + /// assert_eq!(mat.get_mut([3, 0]), None); + /// + /// assert_eq!(*mat.get_mut([0, 0]).unwrap(), 0); + /// *mat.get_mut([0,0]).unwrap() = 2; + /// assert_eq!(*mat.get_mut([0, 0]).unwrap(), 2); + /// # } + /// ``` + fn get_mut(&mut self, index: [usize; 2]) -> Option<&mut T> { + let row_ind = index[0]; + let col_ind = index[1]; + + if row_ind >= self.rows() || col_ind >= self.cols() { + None + } else { + unsafe { Some(self.get_unchecked_mut(index)) } + } + } + /// Returns a mutable iterator over the matrix. /// /// # Examples From 9376d7d5763e8c7d333014ff1c91055cdc6e5a38 Mon Sep 17 00:00:00 2001 From: prismatica Date: Mon, 27 Feb 2017 16:04:03 -0500 Subject: [PATCH 03/12] Changed indentation and comments. As discussed, changed "a point" to "an element" in the comments for the "get" family of functions. Also fixed the indentation of the new "get" and "get_mut" functions by changing tabs to spaces. --- src/matrix/base/mod.rs | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/matrix/base/mod.rs b/src/matrix/base/mod.rs index 79930d3..591ecae 100644 --- a/src/matrix/base/mod.rs +++ b/src/matrix/base/mod.rs @@ -75,12 +75,12 @@ pub trait BaseMatrix: Sized { } } - /// Get a reference to a point in the matrix without bounds checking. + /// Get a reference to an element in the matrix without bounds checking. unsafe fn get_unchecked(&self, index: [usize; 2]) -> &T { &*(self.as_ptr().offset((index[0] * self.row_stride() + index[1]) as isize)) } - /// Get a reference to a point in the matrix. + /// Get a reference to an element in the matrix. /// /// # Examples /// @@ -99,13 +99,13 @@ pub trait BaseMatrix: Sized { /// # } /// ``` fn get(&self, index: [usize; 2]) -> Option<&T> { - let row_ind = index[0]; - let col_ind = index[1]; - + let row_ind = index[0]; + let col_ind = index[1]; + if row_ind >= self.rows() || col_ind >= self.cols() { - None + None } else { - unsafe { Some(self.get_unchecked(index)) } + unsafe { Some(self.get_unchecked(index)) } } } @@ -1168,12 +1168,12 @@ pub trait BaseMatrixMut: BaseMatrix { } } - /// Get a mutable reference to a point in the matrix without bounds checks. + /// Get a mutable reference to an element in the matrix without bounds checks. unsafe fn get_unchecked_mut(&mut self, index: [usize; 2]) -> &mut T { &mut *(self.as_mut_ptr().offset((index[0] * self.row_stride() + index[1]) as isize)) } - /// Get a mutable reference to a point in the matrix. + /// Get a mutable reference to an element in the matrix. /// /// # Examples /// @@ -1194,13 +1194,13 @@ pub trait BaseMatrixMut: BaseMatrix { /// # } /// ``` fn get_mut(&mut self, index: [usize; 2]) -> Option<&mut T> { - let row_ind = index[0]; - let col_ind = index[1]; + let row_ind = index[0]; + let col_ind = index[1]; if row_ind >= self.rows() || col_ind >= self.cols() { - None + None } else { - unsafe { Some(self.get_unchecked_mut(index)) } + unsafe { Some(self.get_unchecked_mut(index)) } } } From f0be95687cfb18a8e010aac25c9001c3956d6f29 Mon Sep 17 00:00:00 2001 From: prismatica Date: Thu, 9 Mar 2017 16:32:17 -0500 Subject: [PATCH 04/12] Added LU decomposition with complete pivoting, resolving #98 Added: * LUPQ, the result of unpacking a complete-pivot LU decomp. * FullPivLu, the parallel for PartialPivLu * A number of tests for full pivoting. --- src/matrix/decomposition/lu.rs | 428 ++++++++++++++++++++++++++++++-- src/matrix/decomposition/mod.rs | 2 +- 2 files changed, 414 insertions(+), 16 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 17fadd8..85fd35c 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -5,6 +5,7 @@ use vector::Vector; use error::{Error, ErrorKind}; use std::any::Any; +use std::cmp; use libnum::{Float, Zero, One}; @@ -181,13 +182,8 @@ impl PartialPivLu { lu.swap_rows(index, curr_max_idx); p.swap_rows(index, curr_max_idx); - for i in (index+1)..n { - let mult = lu[[i, index]] / curr_max; - lu[[i, index]] = mult; - for j in (index+1)..n { - lu[[i, j]] = lu[[i,j]] - mult*lu[[index, j]]; - } - } + + gaussian_elimination(&mut lu, index); } Ok(PartialPivLu { lu: lu, @@ -299,6 +295,260 @@ impl PartialPivLu where T: Any + Float { } } +/// Result of unpacking an instance of +/// [FullPivLu](struct.FullPivLu.html). +/// +/// PAQ = LU +#[derive(Debug, Clone)] +pub struct LUPQ { + /// The lower triangular matrix in the decomposition. + pub l: Matrix, + + /// The upper triangular matrix in the decomposition. + pub u: Matrix, + + /// The row-exchange permutation matrix in the decomposition. + pub p: PermutationMatrix, + + /// The column-exchange permutation matrix in the decomposition. + pub q: PermutationMatrix +} + +/// LU decomposition with complete pivoting. +/// +/// For any square matrix A, there exist two permutation matrices +/// `P` and `Q`, a lower triangular matrix `L` and an upper triangular +/// matrix `U` such that +/// +/// ```text +/// PAQ = LU. +/// ``` +/// +/// Unlike the LU decomposition computed with partial pivoting, this +/// decomposition is stable for singular matrices. It is also a rank- +/// revealing decomposition. +/// +/// See [PartialPivLu](decomposition/struct.PartialPivLu.html) for +/// applications of LU decompositions in general. +#[derive(Debug, Clone)] +pub struct FullPivLu { + lu: Matrix, + p: PermutationMatrix, + q: PermutationMatrix +} + +impl Decomposition for FullPivLu { + type Factors = LUPQ; + + fn unpack(self) -> LUPQ { + let l = unit_lower_triangular_part(&self.lu); + let u = nullify_lower_triangular_part(self.lu); + + LUPQ { + l: l, + u: u, + p: self.p, + q: self.q, + } + } +} + +impl FullPivLu { + fn select_pivot(mat: &Matrix, index: usize) -> (usize, usize, T) { + let mut piv_row = index; + let mut piv_col = index; + let mut piv_val = mat[[index,index]]; + + for row in index..mat.rows() { + for col in index..mat.cols() { + let val = mat[[row,col]]; + + if val.abs() > piv_val.abs() { + piv_val = val; + piv_row = row; + piv_col = col; + } + } + } + + (piv_row, piv_col, piv_val) + } + + /// Performs the decomposition. + pub fn decompose(matrix: Matrix) -> Result { + let mut lu = matrix.clone(); + + let nrows = matrix.rows(); + let ncols = matrix.cols(); + let diag_size = cmp::min(nrows, ncols); + + let mut p = PermutationMatrix::identity(nrows); + let mut q = PermutationMatrix::identity(ncols); + + for index in 0..diag_size { + // Select the current pivot. This is the largest value in + // the bottom right corner of the matrix, starting at + // (index, index). + let (piv_row, piv_col, piv_val) = FullPivLu::select_pivot(&lu, index); + + if piv_val.abs() < T::epsilon() { + break; + } + + lu.swap_rows(index, piv_row); + lu.swap_cols(index, piv_col); + + p.swap_rows(index, piv_row); + + // This is a little misleading, but even though + // we're calling swap_rows here, since q is applied on the + // right to A (i.e. P * A * Q), the result is a column swap of A. + q.swap_rows(index, piv_col); + + // We've swapped the pivot row and column so that the pivot + // ends up in the (index, index) position, so apply gaussian + // elimination to the bottom-right corner. + gaussian_elimination(&mut lu, index); + } + + Ok(FullPivLu { + lu: lu, + p: p.inverse(), + q: q.inverse() + }) + } +} + +// TODO: Remove Any bound (cannot for the time being, since +// back substitution uses Any bound) +impl FullPivLu where T: Any + Float { + + /// Solves the linear system `Ax = b`. + /// + /// Here, `A` is the decomposed matrix satisfying + /// `PAQ = LU`. Note that this method is particularly + /// well suited to solving multiple such linear systems + /// involving the same `A` but different `b`. + /// + /// # Errors + /// + /// If the matrix is very ill-conditioned, the function + /// might fail to obtain the solution to the system. + /// + /// # Panics + /// + /// The right-hand side vector `b` must have compatible size. + /// + /// # Examples + /// + /// ``` + /// # #[macro_use] extern crate rulinalg; + /// # use rulinalg::matrix::decomposition::FullPivLu; + /// # use rulinalg::matrix::Matrix; + /// # fn main() { + /// let x = Matrix::identity(4); + /// let lu = FullPivLu::decompose(x).unwrap(); + /// let b = vector![3.0, 4.0, 2.0, 1.0]; + /// let y = lu.solve(b) + /// .expect("Matrix is invertible."); + /// assert_vector_eq!(y, vector![3.0, 4.0, 2.0, 1.0], comp = float); + /// # } + /// ``` + pub fn solve(&self, b: Vector) -> Result, Error> { + 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) + } + + /// Computes the inverse of the matrix which this LUP decomposition + /// represents. + /// + /// # Errors + /// The inversion might fail if the matrix is very ill-conditioned. + /// The inversion fails if the matrix is not square, or is not invertible. + pub fn inverse(&self) -> Result, Error> { + let n = self.lu.rows(); + let mut inv = Matrix::zeros(n, n); + let mut e = Vector::zeros(n); + + if self.lu.cols() != n { + return Err( + Error::new( + ErrorKind::InvalidArg, + "Only square matrices have inverses")); + } + + if self.det() == T::zero() { + return Err( + Error::new( + ErrorKind::InvalidArg, + "Matrix not invertible")); + } + + for i in 0 .. n { + e[i] = T::one(); + + let col = try!(self.solve(e)); + + for j in 0 .. n { + inv[[j, i]] = col[j]; + } + + e = col.apply(&|_| T::zero()); + } + + Ok(inv) + } + + /// Computes the determinant of the decomposed matrix. + /// + /// # Panics + /// If the underlying matrix is non-square. + pub fn det(&self) -> T { + assert!( + self.lu.rows() == self.lu.cols(), + "Only square matrices have determinants"); + + // Recall that the determinant of a triangular matrix + // is the product of its diagonal entries. Also, + // the determinant of L is implicitly 1. + let u_det = self.lu.diag().fold(T::one(), |x, &y| x * y); + + // Note that the determinants of P and Q are equal to the + // determinant of P^T and Q^T, so we don't have to invert them + let p_det = self.p.clone().det(); + let q_det = self.q.clone().det(); + + p_det * u_det * q_det + } + + /// Computes the rank of the decomposed matrix. + pub fn rank(&self) -> usize { + self.lu.diag().fold( + 0 as usize, + |x, &y| if y.abs() > T::epsilon() { x + 1 } else { x } ) + } +} + +/// Performs Gaussian elimination in the lower-right hand corner starting at +/// (index, index). +fn gaussian_elimination(lu: &mut Matrix, index: usize) { + + let piv_val = lu[[index, index]]; + + for i in (index+1)..lu.rows() { + let mult = lu[[i, index]] / piv_val; + + lu[[i, index]] = mult; + + for j in (index+1)..lu.cols() { + lu[[i, j]] = lu[[i,j]] - mult*lu[[index, j]]; + } + } +} + /// Performs forward substitution using the LU matrix /// for which L has an implicit unit diagonal. That is, /// the strictly lower triangular part of LU corresponds @@ -327,24 +577,23 @@ fn lu_forward_substitution(lu: &Matrix, b: Vector) -> Vector fn unit_lower_triangular_part(matrix: &M) -> Matrix where T: Zero + One + Clone, M: BaseMatrix { - let (m, n) = (matrix.rows(), matrix.cols()); - let mut data = Vec::::with_capacity(m * n); + + let m = matrix.rows(); + let mut data = Vec::::with_capacity(m * m); for (i, row) in matrix.row_iter().enumerate() { for element in row.iter().take(i).cloned() { data.push(element); } - if i < n { - data.push(T::one()); - } + data.push(T::one()); - for _ in (i + 1) .. n { + for _ in (i + 1) .. m { data.push(T::zero()); } } - Matrix::new(m, n, data) + Matrix::new(m, m, data) } fn nullify_lower_triangular_part(mut matrix: Matrix) -> Matrix { @@ -436,7 +685,7 @@ mod tests { use matrix::{Matrix, PermutationMatrix}; use testsupport::{is_lower_triangular, is_upper_triangular}; - use super::{PartialPivLu, LUP}; + use super::{PartialPivLu, LUP, FullPivLu, LUPQ}; use matrix::decomposition::Decomposition; use libnum::Float; @@ -585,4 +834,153 @@ mod tests { } } + #[test] + fn full_piv_lu_decompose_arbitrary() { + // Since the LUP decomposition is not in general unique, + // we can not test against factors directly, but + // instead we must rely on the fact that the + // matrices P, L and U together construct the + // original matrix + let x = matrix![ -3.0, 0.0, 4.0, 1.0; + -12.0, 5.0, 17.0, 1.0; + 15.0, 0.0, -18.0, -5.0; + 6.0, 20.0, -10.0, -15.0 ]; + + let LUPQ { l, u, p, q } = FullPivLu::decompose(x.clone()) + .unwrap() + .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_singular() { + let x = matrix![ -3.0, 0.0, 4.0, 1.0; + -12.0, 5.0, 17.0, 1.0; + 15.0, 0.0, -18.0, -5.0; + -6.0, 0.0, 8.0, 2.0 ]; + + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + assert_eq!(lu.rank(), 3); + + 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_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]; + + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + assert_eq!(lu.rank(), 3); + + 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, 1.0, 20.0]; + + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + assert_eq!(lu.rank(), 2); + + 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_nonsquare_det() { + let x = matrix![ -3.0, 0.0, 4.0; + -6.0, 1.0, 20.0]; + + let lu = FullPivLu::decompose(x.clone()).unwrap(); + + lu.det(); + } + + #[test] + pub fn full_piv_lu_solve_arbitrary_matrix() { + let x = matrix![ 5.0, 0.0, 0.0, 1.0; + 2.0, 2.0, 2.0, 1.0; + 4.0, 5.0, 5.0, 5.0; + 1.0, 6.0, 4.0, 5.0 ]; + let b = vector![9.0, 16.0, 49.0, 45.0]; + let expected = vector![1.0, 2.0, 3.0, 4.0]; + + let lu = FullPivLu::decompose(x).unwrap(); + let y = lu.solve(b).unwrap(); + + // Need to up the tolerance to take into account + // numerical error. Ideally there'd be a more systematic + // way to test this. + assert_vector_eq!(y, expected, comp = ulp, tol = 100); + } + + #[test] + pub fn full_piv_lu_inverse_arbitrary_invertible_matrix() { + let x = matrix![5.0, 0.0, 0.0, 1.0; + 2.0, 2.0, 2.0, 1.0; + 4.0, 5.0, 5.0, 5.0; + 1.0, 6.0, 4.0, 5.0]; + + let inv = matrix![1.85185185185185203e-01, 1.85185185185185175e-01, -7.40740740740740561e-02, -1.02798428206033007e-17; + 1.66666666666666630e-01, 6.66666666666666519e-01, -6.66666666666666519e-01, 4.99999999999999833e-01; + -3.88888888888888840e-01, 1.11111111111111174e-01, 5.55555555555555358e-01, -4.99999999999999833e-01; + 7.40740740740740838e-02, -9.25925925925925819e-01, 3.70370370370370294e-01, 5.13992141030165006e-17]; + + let lu = FullPivLu::decompose(x).unwrap(); + + assert_matrix_eq!(lu.inverse().unwrap(), inv, comp = float); + } + + #[test] + pub fn full_piv_lu_inverse_rectangular() { + let x = matrix![5.0, 0.0, 1.0; + 2.0, 2.0, 1.0; + 4.0, 5.0, 5.0; + 1.0, 6.0, 5.0]; + + let lu = FullPivLu::decompose(x).unwrap(); + + assert!(lu.inverse().is_err()); + } + + #[test] + pub fn full_piv_lu_inverse_noninvertible() { + let x = matrix![5.0, 0.0, 1.0; + 4.0, 5.0, 5.0; + 9.0, 5.0, 6.0]; + + let lu = FullPivLu::decompose(x).unwrap(); + + assert!(lu.inverse().is_err()); + } } diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index afa0e1c..3142c34 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -109,7 +109,7 @@ use vector::Vector; use utils; use error::{Error, ErrorKind}; -pub use self::lu::{PartialPivLu, LUP}; +pub use self::lu::{PartialPivLu, LUP, FullPivLu, LUPQ}; use libnum::{Float}; From 7d0c67b0adbadc28832c355d1df58397c5678087 Mon Sep 17 00:00:00 2001 From: prismatica Date: Sun, 12 Mar 2017 11:27:16 -0400 Subject: [PATCH 05/12] Changed FullPivLu to only work with square matrices --- src/matrix/decomposition/lu.rs | 71 ++++------------------------------ 1 file changed, 7 insertions(+), 64 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 85fd35c..f9cad7c 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -376,6 +376,10 @@ 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.clone(); let nrows = matrix.rows(); @@ -473,14 +477,7 @@ impl FullPivLu where T: Any + Float { let mut inv = Matrix::zeros(n, n); let mut e = Vector::zeros(n); - if self.lu.cols() != n { - return Err( - Error::new( - ErrorKind::InvalidArg, - "Only square matrices have inverses")); - } - - if self.det() == T::zero() { + if self.rank() != n { return Err( Error::new( ErrorKind::InvalidArg, @@ -507,10 +504,6 @@ impl FullPivLu where T: Any + Float { /// # Panics /// If the underlying matrix is non-square. pub fn det(&self) -> T { - assert!( - self.lu.rows() == self.lu.cols(), - "Only square matrices have determinants"); - // Recall that the determinant of a triangular matrix // is the product of its diagonal entries. Also, // the determinant of L is implicitly 1. @@ -878,52 +871,14 @@ 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]; - let lu = FullPivLu::decompose(x.clone()).unwrap(); - - assert_eq!(lu.rank(), 3); - - 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, 1.0, 20.0]; - - let lu = FullPivLu::decompose(x.clone()).unwrap(); - - assert_eq!(lu.rank(), 2); - - 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_nonsquare_det() { - let x = matrix![ -3.0, 0.0, 4.0; - -6.0, 1.0, 20.0]; - - let lu = FullPivLu::decompose(x.clone()).unwrap(); - - lu.det(); + FullPivLu::decompose(x.clone()); } #[test] @@ -961,18 +916,6 @@ mod tests { assert_matrix_eq!(lu.inverse().unwrap(), inv, comp = float); } - #[test] - pub fn full_piv_lu_inverse_rectangular() { - let x = matrix![5.0, 0.0, 1.0; - 2.0, 2.0, 1.0; - 4.0, 5.0, 5.0; - 1.0, 6.0, 5.0]; - - let lu = FullPivLu::decompose(x).unwrap(); - - assert!(lu.inverse().is_err()); - } - #[test] pub fn full_piv_lu_inverse_noninvertible() { let x = matrix![5.0, 0.0, 1.0; From 2ffb026c1d230c487a8a908523f5b2b1fe4bfe6b Mon Sep 17 00:00:00 2001 From: prismatica Date: Sun, 12 Mar 2017 11:33:21 -0400 Subject: [PATCH 06/12] Added FullPivLu to the decomposition table --- src/matrix/decomposition/mod.rs | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index 3142c34..66efe8f 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -79,6 +79,19 @@ //! //! //! +//! +//! FullPivLu +//! Square matrices +//! +//!
    +//!
  • Linear system solving
    • +//!
  • Matrix inverse
    • +//!
  • Determinant computation
    • +//!
  • Rank computation
    • +//!
+//! +//! +//! //! //! From 572d45f277d3f79d679e3ae5887f97a63fc3eb3c Mon Sep 17 00:00:00 2001 From: prismatica Date: Sun, 12 Mar 2017 11:41:10 -0400 Subject: [PATCH 07/12] Fixed a warning I overlooked in a test --- src/matrix/decomposition/lu.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index f9cad7c..5ffbf9c 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -877,8 +877,8 @@ mod tests { -12.0, 5.0, 17.0; 15.0, 0.0, -18.0; -6.0, 0.0, 20.0]; - - FullPivLu::decompose(x.clone()); + + FullPivLu::decompose(x.clone()).unwrap(); } #[test] From aeeb09c6c8c80c532e8577c1c00dd4b479e7748c Mon Sep 17 00:00:00 2001 From: prismatica Date: Sat, 18 Mar 2017 11:41:21 -0400 Subject: [PATCH 08/12] Fix issues brought up in PR 1) Remove unneeded "clone" 2) Fix bad stopping condition 3) Remove unneeded singularity check that is handled in another portion of the code. --- src/matrix/decomposition/lu.rs | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 5ffbf9c..797bbe4 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -380,10 +380,10 @@ impl FullPivLu { matrix.rows() == matrix.cols(), "Matrix must be square for LU decomposition."); - let mut lu = matrix.clone(); + let mut lu = matrix; - let nrows = matrix.rows(); - let ncols = matrix.cols(); + let nrows = lu.rows(); + let ncols = lu.cols(); let diag_size = cmp::min(nrows, ncols); let mut p = PermutationMatrix::identity(nrows); @@ -395,7 +395,7 @@ impl FullPivLu { // (index, index). let (piv_row, piv_col, piv_val) = FullPivLu::select_pivot(&lu, index); - if piv_val.abs() < T::epsilon() { + if piv_val.abs() == T::zero() { break; } @@ -477,12 +477,8 @@ impl FullPivLu where T: Any + Float { let mut inv = Matrix::zeros(n, n); let mut e = Vector::zeros(n); - if self.rank() != n { - return Err( - Error::new( - ErrorKind::InvalidArg, - "Matrix not invertible")); - } + // "solve" will return an error if the matrix is + // singular, so no need to check for singularity here for i in 0 .. n { e[i] = T::one(); From ad9cdc46d6c140957826921854f71fad9724ac8f Mon Sep 17 00:00:00 2001 From: prismatica Date: Tue, 21 Mar 2017 18:21:12 -0400 Subject: [PATCH 09/12] Fix the tolerance for determining rank Instead of directly using T::epsilon to come up with a numerical rank, we now use lu[[0,0]].abs() * T::epsilon. --- src/matrix/decomposition/lu.rs | 67 +++++++++++++++++++++++++++++++--- 1 file changed, 62 insertions(+), 5 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 797bbe4..4ddcdba 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -477,8 +477,12 @@ impl FullPivLu where T: Any + Float { let mut inv = Matrix::zeros(n, n); let mut e = Vector::zeros(n); - // "solve" will return an error if the matrix is - // singular, so no need to check for singularity here + if self.is_singular() { + return Err( + Error::new( + ErrorKind::DivByZero, + "Singular matrix found while attempting inversion.")); + } for i in 0 .. n { e[i] = T::one(); @@ -514,10 +518,63 @@ impl FullPivLu where T: Any + Float { } /// Computes the rank of the decomposed matrix. + /// # Examples + /// + /// ``` + /// # #[macro_use] extern crate rulinalg; + /// # use rulinalg::matrix::decomposition::FullPivLu; + /// # use rulinalg::matrix::Matrix; + /// # fn main() { + /// let x = matrix![1.0, 2.0, 3.0; + /// 4.0, 5.0, 6.0; + /// 5.0, 7.0, 9.0]; + /// let lu = FullPivLu::decompose(x).unwrap(); + /// assert_eq!(lu.rank(), 2); + /// # } + /// ``` pub fn rank(&self) -> usize { - self.lu.diag().fold( - 0 as usize, - |x, &y| if y.abs() > T::epsilon() { x + 1 } else { x } ) + let eps = self.epsilon(); + let mut rank = 0; + + for d in self.lu.diag() { + if d.abs() > eps { + rank = rank + 1; + } else { + break; + } + } + + rank + } + + /// Returns whether the matrix is singular. + /// + /// # Examples + /// + /// ``` + /// # #[macro_use] extern crate rulinalg; + /// # use rulinalg::matrix::decomposition::FullPivLu; + /// # use rulinalg::matrix::Matrix; + /// # fn main() { + /// let x = Matrix::::identity(4); + /// let lu = FullPivLu::decompose(x).unwrap(); + /// assert!(!lu.is_singular()); + /// + /// let y = matrix![1.0, 2.0, 3.0; + /// 4.0, 5.0, 6.0; + /// 5.0, 7.0, 9.0]; + /// let lu = FullPivLu::decompose(y).unwrap(); + /// assert!(lu.is_singular()); + /// # } + /// ``` + pub fn is_singular(&self) -> bool { + let diag_size = cmp::min(self.lu.rows(), self.lu.cols()); + + self.rank() != diag_size + } + + fn epsilon(&self) -> T { + self.lu[[0, 0]].abs() * T::epsilon() } } From e340619e8ce946eb78dd9989b5fea489e9972302 Mon Sep 17 00:00:00 2001 From: prismatica Date: Wed, 22 Mar 2017 17:22:40 -0400 Subject: [PATCH 10/12] Change is_singular to is_invertible, other minor changes --- src/matrix/decomposition/lu.rs | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 4ddcdba..6328f61 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -477,11 +477,11 @@ impl FullPivLu where T: Any + Float { let mut inv = Matrix::zeros(n, n); let mut e = Vector::zeros(n); - if self.is_singular() { + if !self.is_invertible() { return Err( Error::new( ErrorKind::DivByZero, - "Singular matrix found while attempting inversion.")); + "Non-invertible matrix found while attempting inversion")); } for i in 0 .. n { @@ -547,7 +547,10 @@ impl FullPivLu where T: Any + Float { rank } - /// Returns whether the matrix is singular. + /// Returns whether the matrix is invertible. + /// + /// Empty matrices are considered to be invertible for + /// the sake of this function. /// /// # Examples /// @@ -558,23 +561,26 @@ impl FullPivLu where T: Any + Float { /// # fn main() { /// let x = Matrix::::identity(4); /// let lu = FullPivLu::decompose(x).unwrap(); - /// assert!(!lu.is_singular()); + /// assert!(lu.is_invertible()); /// /// let y = matrix![1.0, 2.0, 3.0; /// 4.0, 5.0, 6.0; /// 5.0, 7.0, 9.0]; /// let lu = FullPivLu::decompose(y).unwrap(); - /// assert!(lu.is_singular()); + /// assert!(!lu.is_invertible()); /// # } /// ``` - pub fn is_singular(&self) -> bool { + pub fn is_invertible(&self) -> bool { let diag_size = cmp::min(self.lu.rows(), self.lu.cols()); + let diag_last = diag_size - 1; - self.rank() != diag_size + self.lu.get([diag_last, diag_last]) + .map(|x| x.abs() > self.epsilon()) + .unwrap_or(true) } fn epsilon(&self) -> T { - self.lu[[0, 0]].abs() * T::epsilon() + self.lu.get([0, 0]).unwrap_or(&T::one()).abs() * T::epsilon() } } From 337405ef4ba3adda69698076f7e33de09ea2836b Mon Sep 17 00:00:00 2001 From: prismatica Date: Thu, 23 Mar 2017 18:42:07 -0400 Subject: [PATCH 11/12] Fix bad 'is_invertible' behavior for empty matrix --- src/matrix/decomposition/lu.rs | 46 ++++++++++++++++++++++++++++++---- 1 file changed, 41 insertions(+), 5 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 6328f61..52f0f55 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -471,7 +471,7 @@ impl FullPivLu where T: Any + Float { /// /// # Errors /// The inversion might fail if the matrix is very ill-conditioned. - /// The inversion fails if the matrix is not square, or is not invertible. + /// The inversion fails is not invertible. pub fn inverse(&self) -> Result, Error> { let n = self.lu.rows(); let mut inv = Matrix::zeros(n, n); @@ -501,6 +501,8 @@ impl FullPivLu where T: Any + Float { /// Computes the determinant of the decomposed matrix. /// + /// Empty matrices are considered to have a determinant of 1.0. + /// /// # Panics /// If the underlying matrix is non-square. pub fn det(&self) -> T { @@ -572,11 +574,16 @@ impl FullPivLu where T: Any + Float { /// ``` pub fn is_invertible(&self) -> bool { let diag_size = cmp::min(self.lu.rows(), self.lu.cols()); - let diag_last = diag_size - 1; - self.lu.get([diag_last, diag_last]) - .map(|x| x.abs() > self.epsilon()) - .unwrap_or(true) + if diag_size > 0 { + let diag_last = diag_size - 1; + let last = + unsafe { self.lu.get_unchecked([diag_last, diag_last]) }; + + last.abs() > self.epsilon() + } else { + true + } } fn epsilon(&self) -> T { @@ -985,4 +992,33 @@ mod tests { assert!(lu.inverse().is_err()); } + + #[test] + pub fn full_piv_lu_empty_matrix() { + use matrix::base::BaseMatrix; + + let x = Matrix::from_fn(0, 0, |_, _| 0.0); + assert_eq!(x.rows(), 0); + assert_eq!(x.cols(), 0); + + let lu = FullPivLu::decompose(x).unwrap(); + + assert!(lu.is_invertible()); + assert_eq!(lu.rank(), 0); + assert_eq!(lu.det(), 1.0); + + let inverse = lu.inverse().unwrap(); + assert_eq!(inverse.rows(), 0); + assert_eq!(inverse.cols(), 0); + + let LUPQ { l, u, p, q } = lu.unpack(); + assert_eq!(l.rows(), 0); + assert_eq!(l.cols(), 0); + + assert_eq!(u.rows(), 0); + assert_eq!(u.cols(), 0); + + assert_eq!(p.size(), 0); + assert_eq!(q.size(), 0); + } } From 5b546875ba0dc3a50450994a8077ec5c14a63176 Mon Sep 17 00:00:00 2001 From: prismatica Date: Thu, 23 Mar 2017 18:43:52 -0400 Subject: [PATCH 12/12] Fix a dumb typo in a comment --- src/matrix/decomposition/lu.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 52f0f55..a504516 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -471,7 +471,7 @@ impl FullPivLu where T: Any + Float { /// /// # Errors /// The inversion might fail if the matrix is very ill-conditioned. - /// The inversion fails is not invertible. + /// The inversion fails if the matrix is not invertible. pub fn inverse(&self) -> Result, Error> { let n = self.lu.rows(); let mut inv = Matrix::zeros(n, n);