From 8e8e8fbd3a885073e698a4dc27f6772db77f9c2e Mon Sep 17 00:00:00 2001 From: Andreas Borgen Longva Date: Fri, 17 Feb 2017 21:14:32 +0100 Subject: [PATCH] LU rewrite (#142) * Benchmarks for permutation-vector multiply * Implement is_{lower/upper}_triangular This commit adds a (for now private) testsupport module. Eventually this module may be made available to the public under a feature flag (e.g. `testsupport`). * Initial implementation of PartialPivLu For now it only wraps the existing lup_decomp function, without making any invasive changes. * Implement inverse, solve and det for PartialPivLu * Make `decomposition` a public module and export LU * Make utils module for benchmarks Moved reproducible_random_matrix into this new module (from svd) so that it can be used in other benchmarks. * Benchmarks for LU functionality * Rewrite LU decomposition with PermutationMatrix * Remove some allocation from lu inverse Also add some inline notes on performance potential for inverse() * Store L, U factors in a single matrix LU * Rewrite lu_forward_substitution for auto-vec * Documentation for Decomposition module * Documentation for PartialPivLu * Replace LU error with more informative one * Add reminder for the future on removing Any * Reformulate Factors doc * Deprecate lup_decomp Also switched internal use of lup_decomp with PartialPivLu. * Port existing lup_decomp test to PartialPivLu * Fix typo: compact -> compactly * Fix doc headers (### -> #) * Rephrase sentence in Decomposition docs * Restrict testsupport module to cfg(test) * Fix warnings Allow using deprecated lup_decomp in tests for lup_decomp * Remove redundant parity function --- benches/lib.rs | 7 +- benches/linalg/lu.rs | 140 +++++++++ benches/linalg/permutation.rs | 85 ++++++ benches/linalg/svd.rs | 13 +- benches/linalg/util.rs | 10 + src/lib.rs | 3 + src/matrix/decomposition/lu.rs | 490 +++++++++++++++++++++++++++++++- src/matrix/decomposition/mod.rs | 116 +++++++- src/matrix/impl_mat.rs | 90 ++---- src/matrix/mod.rs | 34 +-- src/testsupport/constraints.rs | 139 +++++++++ src/testsupport/mod.rs | 5 + tests/mat/mod.rs | 84 ++++++ 13 files changed, 1085 insertions(+), 131 deletions(-) create mode 100644 benches/linalg/lu.rs create mode 100644 benches/linalg/permutation.rs create mode 100644 benches/linalg/util.rs create mode 100644 src/testsupport/constraints.rs create mode 100644 src/testsupport/mod.rs diff --git a/benches/lib.rs b/benches/lib.rs index 2f2f551..0699c93 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -5,9 +5,12 @@ extern crate num as libnum; extern crate test; extern crate rand; -mod linalg { +pub mod linalg { mod iter; mod matrix; mod svd; + mod lu; mod norm; -} \ No newline at end of file + mod permutation; + pub mod util; +} diff --git a/benches/linalg/lu.rs b/benches/linalg/lu.rs new file mode 100644 index 0000000..68bdc39 --- /dev/null +++ b/benches/linalg/lu.rs @@ -0,0 +1,140 @@ +use rulinalg::matrix::{BaseMatrix, BaseMatrixMut, + Matrix}; +use rulinalg::matrix::decomposition::{PartialPivLu}; +use rulinalg::vector::Vector; + +use linalg::util::reproducible_random_matrix; +use libnum::{Zero, One}; + +use test::Bencher; + +fn nullify_upper_triangular_part(matrix: &mut Matrix) { + for i in 0 .. matrix.rows() { + for j in (i + 1) .. matrix.cols() { + matrix[[i, j]] = T::zero(); + } + } + + // TODO: assert is_lower_triangular +} + +fn nullify_lower_triangular_part(matrix: &mut Matrix) { + for i in 0 .. matrix.rows() { + for j in 0 .. i { + matrix[[i, j]] = T::zero(); + } + } + + // TODO: assert is_upper_triangular +} + +fn set_diagonal_to_one(matrix: &mut Matrix) { + use rulinalg::matrix::DiagOffset::Main; + for d in matrix.diag_iter_mut(Main) { + *d = T::one(); + } +} + +fn well_conditioned_matrix(n: usize) -> Matrix { + let mut l = reproducible_random_matrix(n, n); + let mut u = reproducible_random_matrix(n, n); + nullify_upper_triangular_part(&mut l); + nullify_lower_triangular_part(&mut u); + // By making the diagonal of both L and U consist + // exclusively of 1, the eigenvalues of the resulting + // matrix LU will also have eigenvalues of 1 + // (and thus be well-conditioned) + set_diagonal_to_one(&mut l); + set_diagonal_to_one(&mut u); + return l * u; +} + +#[bench] +fn partial_piv_lu_decompose_10x10(b: &mut Bencher) { + let n = 10; + let x = well_conditioned_matrix(n); + b.iter(|| { + // Assume that the cost of cloning x is roughly + // negligible in comparison with the cost of LU + PartialPivLu::decompose(x.clone()).expect("Matrix is well-conditioned") + }) +} + +#[bench] +fn partial_piv_lu_decompose_100x100(b: &mut Bencher) { + let n = 100; + let x = well_conditioned_matrix(n); + b.iter(|| { + // Assume that the cost of cloning x is roughly + // negligible in comparison with the cost of LU + PartialPivLu::decompose(x.clone()).expect("Matrix is well-conditioned") + }) +} + +#[bench] +fn partial_piv_lu_inverse_10x10(b: &mut Bencher) { + let n = 10; + let x = well_conditioned_matrix(n); + let lu = PartialPivLu::decompose(x) + .expect("Matrix is well-conditioned."); + b.iter(|| { + lu.inverse() + }) +} + +#[bench] +fn partial_piv_lu_inverse_100x100(b: &mut Bencher) { + let n = 100; + let x = well_conditioned_matrix(n); + let lu = PartialPivLu::decompose(x) + .expect("Matrix is well-conditioned."); + b.iter(|| { + lu.inverse() + }) +} + +#[bench] +fn partial_piv_lu_det_10x10(b: &mut Bencher) { + let n = 10; + let x = well_conditioned_matrix(n); + let lu = PartialPivLu::decompose(x) + .expect("Matrix is well-conditioned."); + b.iter(|| { + lu.det() + }) +} + +#[bench] +fn partial_piv_lu_det_100x100(b: &mut Bencher) { + let n = 100; + let x = well_conditioned_matrix(n); + let lu = PartialPivLu::decompose(x) + .expect("Matrix is well-conditioned."); + b.iter(|| { + lu.det() + }) +} + +#[bench] +fn partial_piv_lu_solve_10x10(b: &mut Bencher) { + let n = 10; + let x = well_conditioned_matrix(n); + let lu = PartialPivLu::decompose(x) + .expect("Matrix is well-conditioned."); + b.iter(|| { + let b = Vector::ones(n); + lu.solve(b).unwrap() + }) +} + +#[bench] +fn partial_piv_lu_solve_100x100(b: &mut Bencher) { + let n = 100; + let x = well_conditioned_matrix(n); + let lu = PartialPivLu::decompose(x) + .expect("Matrix is well-conditioned."); + b.iter(|| { + let b = Vector::ones(n); + lu.solve(b).unwrap() + }) +} diff --git a/benches/linalg/permutation.rs b/benches/linalg/permutation.rs new file mode 100644 index 0000000..00740fa --- /dev/null +++ b/benches/linalg/permutation.rs @@ -0,0 +1,85 @@ +use rulinalg::matrix::PermutationMatrix; +use rulinalg::vector::Vector; + +use test::Bencher; + +/// A perfect shuffle as defined at +/// https://en.wikipedia.org/wiki/Faro_shuffle +/// Note that it creates a permutation matrix of size +/// 2 * n and not n. +fn perfect_shuffle(n: usize) -> PermutationMatrix { + let array = (0 .. 2 * n).map(|k| + if (k + 1) % 2 == 0 { + n + (k + 1) / 2 - 1 + } else { + (k + 1) / 2 + } + ).collect::>(); + PermutationMatrix::from_array(array).unwrap() +} + +#[bench] +fn identity_permutation_mul_vector_100(b: &mut Bencher) { + let n = 100; + let p = PermutationMatrix::identity(n); + let x: Vector = Vector::zeros(n); + + b.iter(|| { + &p * &x + }) +} + +#[bench] +fn identity_permutation_as_matrix_mul_vector_100(b: &mut Bencher) { + let n = 100; + let p = PermutationMatrix::identity(n).as_matrix(); + let x: Vector = Vector::zeros(n); + + b.iter(|| { + &p * &x + }) +} + +#[bench] +fn identity_permutation_mul_vector_1000(b: &mut Bencher) { + let n = 1000; + let p = PermutationMatrix::identity(n); + let x: Vector = Vector::zeros(n); + + b.iter(|| { + &p * &x + }) +} + +#[bench] +fn identity_permutation_as_matrix_mul_vector_1000(b: &mut Bencher) { + let n = 1000; + let p = PermutationMatrix::identity(n).as_matrix(); + let x: Vector = Vector::zeros(n); + + b.iter(|| { + &p * &x + }) +} + +#[bench] +fn perfect_shuffle_permutation_mul_vector_1000(b: &mut Bencher) { + let n = 1000; + let p = perfect_shuffle(500); + let x: Vector = Vector::zeros(n); + + b.iter(|| { + &p * &x + }) +} + +#[bench] +fn perfect_shuffle_permutation_as_matrix_mul_vector_1000(b: &mut Bencher) { + let n = 1000; + let p = perfect_shuffle(500).as_matrix(); + let x: Vector = Vector::zeros(n); + + b.iter(|| { + &p * &x + }) +} diff --git a/benches/linalg/svd.rs b/benches/linalg/svd.rs index e7cc2dc..c08e2a7 100644 --- a/benches/linalg/svd.rs +++ b/benches/linalg/svd.rs @@ -1,14 +1,5 @@ use test::Bencher; -use rand; -use rand::{Rng, SeedableRng}; -use rulinalg::matrix::Matrix; - -fn reproducible_random_matrix(rows: usize, cols: usize) -> Matrix { - const STANDARD_SEED: [usize; 4] = [12, 2049, 4000, 33]; - let mut rng = rand::StdRng::from_seed(&STANDARD_SEED); - let elements: Vec<_> = rng.gen_iter::().take(rows * cols).collect(); - Matrix::new(rows, cols, elements) -} +use linalg::util::reproducible_random_matrix; #[bench] fn svd_10_10(b: &mut Bencher) { @@ -26,4 +17,4 @@ fn svd_100_100(b: &mut Bencher) { b.iter(|| mat.clone().svd() ) -} \ No newline at end of file +} diff --git a/benches/linalg/util.rs b/benches/linalg/util.rs new file mode 100644 index 0000000..5d3a06b --- /dev/null +++ b/benches/linalg/util.rs @@ -0,0 +1,10 @@ +use rulinalg::matrix::Matrix; +use rand; +use rand::{Rng, SeedableRng}; + +pub fn reproducible_random_matrix(rows: usize, cols: usize) -> Matrix { + const STANDARD_SEED: [usize; 4] = [12, 2049, 4000, 33]; + let mut rng = rand::StdRng::from_seed(&STANDARD_SEED); + let elements: Vec<_> = rng.gen_iter::().take(rows * cols).collect(); + Matrix::new(rows, cols, elements) +} diff --git a/src/lib.rs b/src/lib.rs index c3c1b86..54dbb94 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -100,6 +100,9 @@ pub mod vector; pub mod ulp; pub mod norm; +#[cfg(test)] +mod testsupport; + #[cfg(test)] #[macro_use] extern crate quickcheck; diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 7ac64ae..17fadd8 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -1,9 +1,361 @@ -use matrix::{Matrix, BaseMatrixMut}; +use matrix::{Matrix, BaseMatrix, BaseMatrixMut}; +use matrix::{back_substitution}; +use matrix::PermutationMatrix; +use vector::Vector; use error::{Error, ErrorKind}; use std::any::Any; -use libnum::Float; +use libnum::{Float, Zero, One}; + +use matrix::decomposition::Decomposition; + +/// Result of unpacking an instance of +/// [PartialPivLu](struct.PartialPivLu.html). +#[derive(Debug, Clone)] +pub struct LUP { + /// The lower triangular matrix in the decomposition. + pub l: Matrix, + /// The upper triangular matrix in the decomposition. + pub u: Matrix, + /// The permutation matrix in the decomposition. + pub p: PermutationMatrix +} + +/// LU decomposition with partial pivoting. +/// +/// For any square matrix A, there exist a permutation matrix +/// `P`, a lower triangular matrix `L` and an upper triangular +/// matrix `U` such that +/// +/// ```text +/// PA = LU. +/// ``` +/// +/// However, due to the way partial pivoting algorithms work, +/// LU decomposition with partial pivoting is in general +/// *only numerically stable for well-conditioned invertible matrices*. +/// +/// That said, partial pivoting is sufficient in the vast majority +/// of practical applications, and it is also the fastest of the +/// pivoting schemes in existence. +/// +/// +/// # Applications +/// +/// Given a matrix `x`, computing the LU(P) decomposition is simple: +/// +/// ``` +/// use rulinalg::matrix::decomposition::{PartialPivLu, LUP, Decomposition}; +/// use rulinalg::matrix::Matrix; +/// +/// let x = Matrix::::identity(4); +/// +/// // The matrix is consumed and its memory +/// // re-purposed for the decomposition +/// let lu = PartialPivLu::decompose(x).expect("Matrix is invertible."); +/// +/// // See below for applications +/// // ... +/// +/// // The factors L, U and P can be obtained by unpacking the +/// // decomposition, for example by destructuring as seen here +/// let LUP { l, u, p } = lu.unpack(); +/// +/// ``` +/// +/// ## Solving linear systems +/// +/// Arguably the most common use case of LU decomposition +/// is the computation of solutions to (multiple) linear systems +/// that share the same coefficient matrix. +/// +/// ``` +/// # #[macro_use] extern crate rulinalg; +/// # use rulinalg::matrix::decomposition::PartialPivLu; +/// # use rulinalg::matrix::Matrix; +/// # fn main() { +/// # let x = Matrix::identity(4); +/// # let lu = PartialPivLu::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); +/// +/// // We can efficiently solve multiple such systems +/// let c = vector![0.0, 0.0, 0.0, 0.0]; +/// let z = lu.solve(c).unwrap(); +/// assert_vector_eq!(z, vector![0.0, 0.0, 0.0, 0.0], comp = float); +/// # } +/// ``` +/// +/// ## Computing the inverse of a matrix +/// +/// The LU decomposition provides a convenient way to obtain +/// the inverse of the decomposed matrix. However, please keep +/// in mind that explicitly computing the inverse of a matrix +/// is *usually* a bad idea. In many cases, one might instead simply +/// solve multiple systems using `solve`. +/// +/// For example, a common misconception is that when one needs +/// to solve multiple linear systems `Ax = b` for different `b`, +/// one should pre-compute the inverse of the matrix for efficiency. +/// In fact, this is practically never a good idea! A far more efficient +/// and accurate method is to perform the LU decomposition once, and +/// then solve each system as shown in the examples of the previous +/// subsection. +/// +/// That said, there are definitely cases where an explicit inverse is +/// needed. In these cases, the inverse can easily be obtained +/// through the `inverse()` method. +/// +/// # Computing the determinant of a matrix +/// +/// Once the LU decomposition has been obtained, computing +/// the determinant of the decomposed matrix is a very cheap +/// operation. +/// +/// ``` +/// # #[macro_use] extern crate rulinalg; +/// # use rulinalg::matrix::decomposition::PartialPivLu; +/// # use rulinalg::matrix::Matrix; +/// # fn main() { +/// # let x = Matrix::::identity(4); +/// # let lu = PartialPivLu::decompose(x).unwrap(); +/// assert_eq!(lu.det(), 1.0); +/// # } +/// ``` +#[derive(Debug, Clone)] +pub struct PartialPivLu { + lu: Matrix, + p: PermutationMatrix +} + +impl Decomposition for PartialPivLu { + type Factors = LUP; + + fn unpack(self) -> LUP { + let l = unit_lower_triangular_part(&self.lu); + let u = nullify_lower_triangular_part(self.lu); + + LUP { + l: l, + u: u, + p: self.p + } + } +} + +impl PartialPivLu { + /// Performs the decomposition. + /// + /// # Panics + /// + /// The matrix must be square. + /// + /// # Errors + /// + /// An error will be returned if the matrix + /// is singular to working precision (badly conditioned). + pub fn decompose(matrix: Matrix) -> Result { + let n = matrix.cols; + assert!(matrix.rows == n, "Matrix must be square for LU decomposition."); + let mut lu = matrix; + let mut p = PermutationMatrix::identity(n); + + for index in 0..n { + let mut curr_max_idx = index; + let mut curr_max = lu[[curr_max_idx, curr_max_idx]]; + + for i in (curr_max_idx+1)..n { + if lu[[i, index]].abs() > curr_max.abs() { + curr_max = lu[[i, index]]; + curr_max_idx = i; + } + } + if curr_max.abs() < T::epsilon() { + return Err(Error::new(ErrorKind::DivByZero, + "The matrix is too ill-conditioned for + LU decomposition with partial pivoting.")); + } + + 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]]; + } + } + } + Ok(PartialPivLu { + lu: lu, + p: p.inverse() + }) + } +} + +// TODO: Remove Any bound (cannot for the time being, since +// back substitution uses Any bound) +impl PartialPivLu where T: Any + Float { + /// Solves the linear system `Ax = b`. + /// + /// Here, `A` is the decomposed matrix satisfying + /// `PA = 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::PartialPivLu; + /// # use rulinalg::matrix::Matrix; + /// # fn main() { + /// let x = Matrix::identity(4); + /// let lu = PartialPivLu::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."); + // Note that applying p here implicitly incurs a clone. + // TODO: Is it possible to avoid the clone somehow? + // To my knowledge, applying a permutation matrix + // in-place in O(n) time requires O(n) storage for bookkeeping. + // However, we might be able to get by with something like + // O(n log n) for the permutation as the forward/backward + // substitution algorithms are O(n^2), if this helps us + // avoid the memory overhead. + let b = lu_forward_substitution(&self.lu, &self.p * b); + back_substitution(&self.lu, b) + } + + /// Computes the inverse of the matrix which this LUP decomposition + /// represents. + /// + /// # Errors + /// The inversion might fail if the matrix is very ill-conditioned. + pub fn inverse(&self) -> Result, Error> { + let n = self.lu.rows(); + let mut inv = Matrix::zeros(n, n); + let mut e = Vector::zeros(n); + + // To compute the inverse of a matrix A, note that + // we can simply solve the system + // AX = I, + // where X is the inverse of A, and I is the identity + // matrix of appropriate dimension. + // + // Note that this is not optimal in terms of performance, + // and there is likely significant potential for improvement. + // + // A more performant technique is usually to compute the + // triangular inverse of each of the L and U triangular matrices, + // but this again requires efficient algorithms (blocked/recursive) + // to invert triangular matrices, which at this point + // we do not have available. + + // Solve for each column of the inverse matrix + 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. + pub fn det(&self) -> T { + // 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 determinant of P is equal to the + // determinant of P^T, so we don't have to invert it + let p_det = self.p.clone().det(); + p_det * u_det + } +} + +/// 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 +/// 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 { + 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; + + 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() + .iter() + .take(i) + .cloned() + .zip(x.iter().cloned()) + .fold(T::zero(), |sum, (l, x)| sum + l * x); + + x[i] = x[i] - adjustment; + } + x +} + +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); + + 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()); + } + + for _ in (i + 1) .. n { + data.push(T::zero()); + } + } + + Matrix::new(m, n, data) +} + +fn nullify_lower_triangular_part(mut matrix: Matrix) -> Matrix { + for (i, mut row) in matrix.row_iter_mut().enumerate() { + for element in row.iter_mut().take(i) { + *element = T::zero(); + } + } + matrix +} + impl Matrix where T: Any + Float { @@ -11,6 +363,10 @@ impl Matrix where T: Any + Float /// /// Returns L,U, and P respectively. /// + /// This function is deprecated. + /// Please see [PartialPivLu](decomposition/struct.PartialPivLu.html) + /// for a replacement. + /// /// # Examples /// /// ``` @@ -32,6 +388,7 @@ impl Matrix where T: Any + Float /// # Failures /// /// - Matrix cannot be LUP decomposed. + #[deprecated] pub fn lup_decomp(self) -> Result<(Matrix, Matrix, Matrix), Error> { let n = self.cols; assert!(self.rows == n, "Matrix must be square for LUP decomposition."); @@ -76,8 +433,15 @@ impl Matrix where T: Any + Float #[cfg(test)] mod tests { - use matrix::Matrix; + use matrix::{Matrix, PermutationMatrix}; + use testsupport::{is_lower_triangular, is_upper_triangular}; + + use super::{PartialPivLu, LUP}; + use matrix::decomposition::Decomposition; + + use libnum::Float; + #[allow(deprecated)] #[test] #[should_panic] fn test_non_square_lup_decomp() { @@ -86,6 +450,7 @@ mod tests { let _ = a.lup_decomp(); } + #[allow(deprecated)] #[test] fn test_lup_decomp() { use error::ErrorKind; @@ -101,4 +466,123 @@ mod tests { Ok(_) => panic!() } } + + #[test] + fn partial_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 LUP { l, u, p } = PartialPivLu::decompose(x.clone()) + .unwrap() + .unpack(); + let y = p.inverse() * &l * &u; + assert_matrix_eq!(x, y, comp = float); + assert!(is_lower_triangular(&l)); + assert!(is_upper_triangular(&u)); + } + + #[test] + pub fn partial_piv_lu_inverse_identity() { + let lu = PartialPivLu:: { + lu: Matrix::identity(3), + p: PermutationMatrix::identity(3) + }; + + let inv = lu.inverse().expect("Matrix is invertible."); + + assert_matrix_eq!(inv, Matrix::identity(3), comp = float); + } + + #[test] + pub fn partial_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 = PartialPivLu::decompose(x).unwrap(); + + assert_matrix_eq!(lu.inverse().unwrap(), inv, comp = float); + } + + #[test] + pub fn partial_piv_lu_det_identity() { + let lu = PartialPivLu:: { + lu: Matrix::identity(3), + p: PermutationMatrix::identity(3) + }; + + assert_eq!(lu.det(), 1.0); + } + + #[test] + pub fn partial_piv_lu_det_arbitrary_invertible_matrix() { + let x = matrix![ 5.0, 0.0, 0.0, 1.0; + 0.0, 2.0, 2.0, 1.0; + 15.0, 4.0, 7.0, 10.0; + 5.0, 2.0, 17.0, 32.0]; + + let lu = PartialPivLu::decompose(x).unwrap(); + + let expected_det = 149.99999999999997; + let diff = lu.det() - expected_det; + assert!(diff.abs() < 1e-12); + } + + #[test] + pub fn partial_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 = PartialPivLu::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 lu_forward_substitution() { + use super::lu_forward_substitution; + + { + let lu: Matrix = matrix![]; + let b = vector![]; + let x = lu_forward_substitution(&lu, b); + assert!(x.size() == 0); + } + + { + let lu = matrix![3.0]; + let b = vector![1.0]; + let x = lu_forward_substitution(&lu, b); + assert_eq!(x, vector![1.0]); + } + + { + let lu = matrix![3.0, 2.0; + 2.0, 2.0]; + let b = vector![1.0, 2.0]; + let x = lu_forward_substitution(&lu, b); + assert_eq!(x, vector![1.0, 0.0]); + } + } + } diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index 56b21a9..afa0e1c 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -1,14 +1,97 @@ -//! Matrix Decompositions +//! Decompositions for matrices. //! -//! References: -//! 1. [On Matrix Balancing and EigenVector computation] -//! (http://arxiv.org/pdf/1401.5766v1.pdf), James, Langou and Lowery +//! This module houses the decomposition API of `rulinalg`. +//! A decomposition - or factorization - of a matrix is an +//! ordered set of *factors* such that when multiplied reconstructs +//! the original matrix. The [Decomposition](trait.Decomposition.html) +//! trait encodes this property. //! -//! 2. [The QR algorithm for eigen decomposition] -//! (http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf) +//! # The decomposition API //! -//! 3. [Computation of the SVD] -//! (http://www.cs.utexas.edu/users/inderjit/public_papers/HLA_SVD.pdf) +//! Decompositions in `rulinalg` are in general modeled after +//! the following: +//! +//! 1. Given an appropriate matrix, an opaque decomposition object +//! may be computed which internally stores the factors +//! in an efficient and appropriate format. +//! 2. In general, the factors may not be immediately available +//! as distinct matrices after decomposition. If the user +//! desires the explicit matrix factors involved in the +//! decomposition, the user must `unpack` the decomposition. +//! 3. Before unpacking the decomposition, the decomposition +//! data structure in question may offer an API that provides +//! efficient implementations for some of the most common +//! applications of the decomposition. The user is encouraged +//! to use the decomposition-specific API rather than unpacking +//! the decompositions whenever possible. +//! +//! For a motivating example that explains the rationale behind +//! this design, let us consider the typical LU decomposition with +//! partial pivoting. In this case, given a square invertible matrix +//! `A`, one may find matrices `P`, `L` and `U` such that +//! `PA = LU`. Here `P` is a permutation matrix, `L` is a lower +//! triangular matrix and `U` is an upper triangular matrix. +//! +//! Once the decomposition has been obtained, one of its applications +//! is the efficient solution of multiple similar linear systems. +//! Consider that while computing the LU decomposition requires +//! O(n3) floating point operations, the solution to +//! the system `Ax = b` can be computed in O(n2) floating +//! point operations if the LU decomposition has already been obtained. +//! Since the right-hand side `b` has no bearing on the LU decomposition, +//! it follows that one can efficiently solve this system for any `b`. +//! +//! It turns out that the matrices `L` and `U` can be stored compactly +//! in the space of a single matrix. Indeed, this is how `PartialPivLu` +//! stores the LU decomposition internally. This allows `rulinalg` to +//! provide the user with efficient implementations of common applications +//! for the LU decomposition. However, the full matrix factors are easily +//! available to the user by unpacking the decomposition. +//! +//! # Available decompositions +//! +//! **The decompositions API is a work in progress.** +//! +//! Currently, only a portion of the available decompositions in `rulinalg` +//! are available through the decomposition API. Please see the +//! [Matrix](../struct.Matrix.html) API for the old decomposition +//! implementations that have yet not been implemented within +//! this framework. +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//! +//!
DecompositionApplicable toSupported features
PartialPivLuSquare, invertible matrices +//!
    +//!
  • Linear system solving
    • +//!
  • Matrix inverse
    • +//!
  • Determinant computation
    • +//!
+//!
+ +// References: +// +// 1. [On Matrix Balancing and EigenVector computation] +// (http://arxiv.org/pdf/1401.5766v1.pdf), James, Langou and Lowery +// +// 2. [The QR algorithm for eigen decomposition] +// (http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf) +// +// 3. [Computation of the SVD] +// (http://www.cs.utexas.edu/users/inderjit/public_papers/HLA_SVD.pdf) mod qr; mod cholesky; @@ -26,8 +109,25 @@ use vector::Vector; use utils; use error::{Error, ErrorKind}; +pub use self::lu::{PartialPivLu, LUP}; + use libnum::{Float}; +/// Base trait for decompositions. +/// +/// A matrix decomposition, or factorization, +/// is a procedure which takes a matrix `X` and returns +/// a set of `k` factors `X_1, X_2, ..., X_k` such that +/// `X = X_1 * X_2 * ... * X_k`. +pub trait Decomposition { + /// The type representing the ordered set of factors + /// that when multiplied yields the decomposed matrix. + type Factors; + + /// Extract the individual factors from this decomposition. + fn unpack(self) -> Self::Factors; +} + impl Matrix where T: Any + Float { diff --git a/src/matrix/impl_mat.rs b/src/matrix/impl_mat.rs index 50ee670..cf71364 100644 --- a/src/matrix/impl_mat.rs +++ b/src/matrix/impl_mat.rs @@ -2,11 +2,12 @@ use std::any::Any; use std::fmt; use libnum::{One, Zero, Float, FromPrimitive}; -use super::{Matrix, forward_substitution, back_substitution, parity}; +use super::{Matrix}; use super::{Axes}; use super::base::BaseMatrix; use error::{Error, ErrorKind}; use vector::Vector; +use matrix::decomposition::PartialPivLu; impl Matrix { /// Constructor for Matrix struct. @@ -314,6 +315,12 @@ impl Matrix { /// /// Requires a Vector `y` as input. /// + /// The method performs an LU decomposition internally, + /// consuming the matrix in the process. If solving + /// the same system for multiple right-hand sides + /// is desired, see + /// [PartialPivLu](decomposition/struct.PartialPivLu.html). + /// /// # Examples /// /// ``` @@ -341,14 +348,13 @@ impl Matrix { /// - The matrix cannot be decomposed into an LUP form to solve. /// - There is no valid solution as the matrix is singular. pub fn solve(self, y: Vector) -> Result, Error> { - let (l, u, p) = try!(self.lup_decomp()); - - let b = try!(forward_substitution(&l, p * y)); - back_substitution(&u, b) + PartialPivLu::decompose(self)?.solve(y) } /// Computes the inverse of the matrix. /// + /// Internally performs an LU decomposition. + /// /// # Examples /// /// ``` @@ -374,43 +380,7 @@ impl Matrix { /// - The matrix could not be LUP decomposed. /// - The matrix has zero determinant. pub fn inverse(self) -> Result, Error> { - let rows = self.rows; - let cols = self.cols; - assert!(rows == cols, "Matrix is not square."); - - let mut inv_t_data = Vec::::new(); - let (l, u, p) = try!(self.lup_decomp().map_err(|_| { - Error::new(ErrorKind::DecompFailure, - "Could not compute LUP factorization for inverse.") - })); - - let mut d = T::one(); - - unsafe { - for i in 0..l.cols { - d = d * *l.get_unchecked([i, i]); - d = d * *u.get_unchecked([i, i]); - } - } - - if d == T::zero() { - return Err(Error::new(ErrorKind::DecompFailure, - "Matrix is singular and cannot be inverted.")); - } - - for i in 0..rows { - let mut id_col = vec![T::zero(); cols]; - id_col[i] = T::one(); - - let b = forward_substitution(&l, &p * Vector::new(id_col)) - .expect("Matrix is singular AND has non-zero determinant!?"); - inv_t_data.append(&mut back_substitution(&u, b) - .expect("Matrix is singular AND has non-zero determinant!?") - .into_vec()); - - } - - Ok(Matrix::new(rows, cols, inv_t_data).transpose()) + PartialPivLu::decompose(self)?.inverse() } /// Computes the determinant of the matrix. @@ -436,18 +406,8 @@ impl Matrix { let n = self.cols; if self.is_diag() { - let mut d = T::one(); - - unsafe { - for i in 0..n { - d = d * *self.get_unchecked([i, i]); - } - } - - return d; - } - - if n == 2 { + self.diag().cloned().fold(T::one(), |d, entry| d * entry) + } else if n == 2 { (self[[0, 0]] * self[[1, 1]]) - (self[[0, 1]] * self[[1, 0]]) } else if n == 3 { (self[[0, 0]] * self[[1, 1]] * self[[2, 2]]) + @@ -457,26 +417,8 @@ impl Matrix { (self[[0, 1]] * self[[1, 0]] * self[[2, 2]]) - (self[[0, 2]] * self[[1, 1]] * self[[2, 0]]) } else { - let (l, u, p) = match self.lup_decomp() { - Ok(x) => x, - Err(ref e) if *e.kind() == ErrorKind::DivByZero => return T::zero(), - _ => { - panic!("Could not compute LUP decomposition."); - } - }; - - let mut d = T::one(); - - unsafe { - for i in 0..l.cols { - d = d * *l.get_unchecked([i, i]); - d = d * *u.get_unchecked([i, i]); - } - } - - let sgn = parity(&p); - - sgn * d + PartialPivLu::decompose(self).map(|lu| lu.det()) + .unwrap_or(T::zero()) } } } diff --git a/src/matrix/mod.rs b/src/matrix/mod.rs index e446584..c74e93a 100644 --- a/src/matrix/mod.rs +++ b/src/matrix/mod.rs @@ -12,11 +12,10 @@ use std::marker::PhantomData; use libnum::Float; use error::{Error, ErrorKind}; -use utils; use vector::Vector; +pub mod decomposition; mod base; -mod decomposition; mod deref; mod impl_mat; mod impl_ops; @@ -371,34 +370,3 @@ fn forward_substitution(m: &M, y: Vector) -> Result, Error> Ok(Vector::new(x)) } - -/// Computes the parity of a permutation matrix. -fn parity(m: &M) -> T - where T: Any + Float, - M: BaseMatrix -{ - let mut visited = vec![false; m.rows()]; - let mut sgn = T::one(); - - for k in 0..m.rows() { - if !visited[k] { - let mut next = k; - let mut len = 0; - - while !visited[next] { - len += 1; - visited[next] = true; - unsafe { - next = utils::find(&m.row_unchecked(next) - .raw_slice(), - T::one()); - } - } - - if len % 2 == 0 { - sgn = -sgn; - } - } - } - sgn -} diff --git a/src/testsupport/constraints.rs b/src/testsupport/constraints.rs new file mode 100644 index 0000000..4078645 --- /dev/null +++ b/src/testsupport/constraints.rs @@ -0,0 +1,139 @@ +use matrix::BaseMatrix; + +use libnum::Zero; + +use std::iter::Iterator; + +/// Returns true if the matrix is lower triangular, otherwise false. +/// This generalizes to rectangular matrices, in which case +/// it returns true if the matrix is lower trapezoidal. +pub fn is_lower_triangular(m: &M) -> bool + where T: Zero + PartialEq, M: BaseMatrix { + + m.row_iter() + .enumerate() + .all(|(i, row)| row.iter() + .skip(i + 1) + .all(|x| x == &T::zero())) +} + +/// Returns true if the matrix is upper triangular, otherwise false. +/// This generalizes to rectangular matrices, in which case +/// it returns true if the matrix is upper trapezoidal. +pub fn is_upper_triangular(m: &M) -> bool + where T: Zero + PartialEq, M: BaseMatrix { + + m.row_iter() + .enumerate() + .all(|(i, row)| row.iter() + .take(i) + .all(|x| x == &T::zero())) +} + +#[cfg(test)] +mod tests { + use super::is_lower_triangular; + use super::is_upper_triangular; + use matrix::Matrix; + + #[test] + fn is_lower_triangular_empty_matrix() { + let x: Matrix = matrix![]; + assert!(is_lower_triangular(&x)); + } + + #[test] + fn is_lower_triangular_1x1() { + let x = matrix![1]; + assert!(is_lower_triangular(&x)); + } + + #[test] + fn is_lower_triangular_square() { + { + let x = matrix![3, 0; + 4, 5]; + assert!(is_lower_triangular(&x)); + } + + { + let x = matrix![1, 0, 0; + 3, 3, 0; + 0, 4, 6]; + assert!(is_lower_triangular(&x)); + } + + { + // Diagonal is an important special case + let x = matrix![1, 0; + 0, 2]; + assert!(is_lower_triangular(&x)); + } + } + + #[test] + fn is_upper_triangular_empty_matrix() { + let x: Matrix = matrix![]; + assert!(is_upper_triangular(&x)); + } + + #[test] + fn is_upper_triangular_1x1() { + let x = matrix![1]; + assert!(is_upper_triangular(&x)); + } + + #[test] + fn is_upper_triangular_square() { + { + let x = matrix![3, 4; + 0, 5]; + assert!(is_upper_triangular(&x)); + } + + { + let x = matrix![1, 3, 0; + 0, 3, 4; + 0, 0, 6]; + assert!(is_upper_triangular(&x)); + } + + { + // Diagonal is an important special case + let x = matrix![1, 0; + 0, 2]; + assert!(is_upper_triangular(&x)); + } + } + + #[test] + fn is_upper_triangular_rectangular() { + { + let x = matrix![1, 2; + 0, 3; + 0, 0]; + assert!(is_upper_triangular(&x)); + } + + { + let x = matrix![1, 2, 3; + 0, 4, 5]; + assert!(is_upper_triangular(&x)); + } + } + + quickcheck! { + fn property_zero_is_lower_triangular(m: usize, n: usize) -> bool { + let x = Matrix::::zeros(m, n); + is_lower_triangular(&x) + } + } + + quickcheck! { + fn property_zero_is_upper_triangular(m: usize, n: usize) -> bool { + let x = Matrix::::zeros(m, n); + is_upper_triangular(&x) + } + } +} + diff --git a/src/testsupport/mod.rs b/src/testsupport/mod.rs new file mode 100644 index 0000000..3277d1f --- /dev/null +++ b/src/testsupport/mod.rs @@ -0,0 +1,5 @@ +//! Provides support functions for writing +//! tests involving data structures provided by rulinalg. +mod constraints; + +pub use self::constraints::{is_lower_triangular, is_upper_triangular}; diff --git a/tests/mat/mod.rs b/tests/mat/mod.rs index 1100e4c..b71b860 100644 --- a/tests/mat/mod.rs +++ b/tests/mat/mod.rs @@ -39,6 +39,7 @@ fn test_u_triangular_solve_errs() { assert!(a.solve_u_triangular(vector![1.0]).is_err()); } +#[allow(deprecated)] #[test] fn matrix_lup_decomp() { let a = matrix![1., 3., 5.; @@ -88,6 +89,89 @@ fn matrix_lup_decomp() { assert!(d.lup_decomp().is_ok()); } +#[test] +fn matrix_partial_piv_lu() { + use rulinalg::matrix::decomposition::{LUP, PartialPivLu}; + use rulinalg::matrix::decomposition::Decomposition; + // This is a port of the test for the old lup_decomp + // function, using the new PartialPivLu struct. + + { + // Note: this test only works with a _specific_ + // implementation of LU decomposition, because + // an LUP decomposition is not in general unique, + // so one cannot test against expected L, U and P + // matrices unless one knows exactly which ones the + // algorithm works. + let a = matrix![1., 3., 5.; + 2., 4., 7.; + 1., 1., 0.]; + + let LUP { l, u, p } = PartialPivLu::decompose(a) + .expect("Matrix is well-conditioned") + .unpack(); + + let l_true = matrix![1.0, 0.0, 0.0; + 0.5, 1.0, 0.0; + 0.5, -1.0, 1.0]; + let u_true = matrix![2.0, 4.0, 7.0; + 0.0, 1.0, 1.5; + 0.0, 0.0, -2.0]; + let p_true = matrix![0.0, 1.0, 0.0; + 1.0, 0.0, 0.0; + 0.0, 0.0, 1.0]; + + assert_matrix_eq!(l, l_true, comp = float); + assert_matrix_eq!(u, u_true, comp = float); + assert_matrix_eq!(p.as_matrix(), p_true, comp = float); + } + + { + let b = matrix![1., 2., 3., 4., 5.; + 3., 0., 4., 5., 6.; + 2., 1., 2., 3., 4.; + 0., 0., 0., 6., 5.; + 0., 0., 0., 5., 6.]; + + let LUP { l, u, p } = PartialPivLu::decompose(b.clone()) + .expect("Matrix is well-conditioned.") + .unpack(); + let k = p.inverse() * l * u; + + assert_matrix_eq!(k, b, comp = float); + } + + { + let c = matrix![-4.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0; + 1.0, -4.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0; + 0.0, 1.0, -4.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0; + 1.0, 0.0, 0.0, -4.0, 1.0, 0.0, 1.0, 0.0, 0.0; + 0.0, 1.0, 0.0, 1.0, -4.0, 1.0, 0.0, 1.0, 0.0; + 0.0, 0.0, 1.0, 0.0, 1.0, -4.0, 0.0, 0.0, 1.0; + 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, -4.0, 1.0, 0.0; + 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, -4.0, 1.0; + 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, -4.0]; + + let LUP { l, u, p } = PartialPivLu::decompose(c.clone()) + .unwrap() + .unpack(); + let c_reconstructed = p.inverse() * l * u; + assert_matrix_eq!(c_reconstructed, c, comp = float); + } + + { + let d = matrix![1.0, 1.0, 0.0, 0.0; + 0.0, 0.0, 1.0, 0.0; + -1.0, 0.0, 0.0, 0.0; + 0.0, 0.0, 0.0, 1.0]; + let LUP { l, u, p } = PartialPivLu::decompose(d.clone()) + .unwrap() + .unpack(); + let d_reconstructed = p.inverse() * l * u; + assert_matrix_eq!(d_reconstructed, d, comp = float); + } +} + #[test] fn cholesky() { let a = matrix![25., 15., -5.;