diff --git a/benches/lib.rs b/benches/lib.rs index be53bfd..f313c99 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -11,6 +11,7 @@ pub mod linalg { mod matrix; mod svd; mod lu; + mod cholesky; mod norm; mod triangular; mod permutation; diff --git a/benches/linalg/cholesky.rs b/benches/linalg/cholesky.rs new file mode 100644 index 0000000..f7b6797 --- /dev/null +++ b/benches/linalg/cholesky.rs @@ -0,0 +1,69 @@ +use rulinalg::matrix::Matrix; +use rulinalg::matrix::decomposition::{Cholesky, Decomposition}; +use test::Bencher; + +#[bench] +fn cholesky_decompose_unpack_100x100(b: &mut Bencher) { + let n = 100; + let x = Matrix::::identity(n); + b.iter(|| { + // Assume that the cost of cloning x is roughly + // negligible in comparison with the cost of LU + Cholesky::decompose(x.clone()).expect("Matrix is invertible") + .unpack() + }) +} + +#[bench] +fn cholesky_decompose_unpack_500x500(b: &mut Bencher) { + let n = 500; + let x = Matrix::::identity(n); + b.iter(|| { + // Assume that the cost of cloning x is roughly + // negligible in comparison with the cost of LU + Cholesky::decompose(x.clone()).expect("Matrix is invertible") + .unpack() + }) +} + +#[bench] +fn cholesky_100x100(b: &mut Bencher) { + // Benchmark for legacy cholesky(). Remove when + // cholesky() has been removed. + let n = 100; + let x = Matrix::::identity(n); + b.iter(|| { + x.cholesky().expect("Matrix is invertible") + }) +} + +#[bench] +fn cholesky_500x500(b: &mut Bencher) { + // Benchmark for legacy cholesky(). Remove when + // cholesky() has been removed. + let n = 500; + let x = Matrix::::identity(n); + b.iter(|| { + x.cholesky().expect("Matrix is invertible") + }) +} + +#[bench] +fn cholesky_solve_1000x1000(b: &mut Bencher) { + let n = 1000; + let x = Matrix::identity(n); + let cholesky = Cholesky::decompose(x).unwrap(); + b.iter(|| { + cholesky.solve(vector![0.0; n]) + }); +} + +#[bench] +fn cholesky_solve_100x100(b: &mut Bencher) { + let n = 100; + let x = Matrix::identity(n); + let cholesky = Cholesky::decompose(x).unwrap(); + b.iter(|| { + cholesky.solve(vector![0.0; n]) + }); +} diff --git a/src/internal_utils.rs b/src/internal_utils.rs new file mode 100644 index 0000000..a1faaaf --- /dev/null +++ b/src/internal_utils.rs @@ -0,0 +1,52 @@ +use matrix::BaseMatrixMut; +use libnum::Zero; + +pub fn nullify_lower_triangular_part(matrix: &mut M) + where T: Zero, M: BaseMatrixMut { + for (i, mut row) in matrix.row_iter_mut().enumerate() { + for element in row.raw_slice_mut().iter_mut().take(i) { + *element = T::zero(); + } + } +} + +pub fn nullify_upper_triangular_part(matrix: &mut M) + where T: Zero, M: BaseMatrixMut { + for (i, mut row) in matrix.row_iter_mut().enumerate() { + for element in row.raw_slice_mut().iter_mut().skip(i + 1) { + *element = T::zero(); + } + } +} + +#[cfg(test)] +mod tests { + use super::nullify_lower_triangular_part; + use super::nullify_upper_triangular_part; + + #[test] + fn nullify_lower_triangular_part_examples() { + let mut x = matrix![1.0, 2.0, 3.0; + 4.0, 5.0, 6.0; + 7.0, 8.0, 9.0]; + nullify_lower_triangular_part(&mut x); + assert_matrix_eq!(x, matrix![ + 1.0, 2.0, 3.0; + 0.0, 5.0, 6.0; + 0.0, 0.0, 9.0 + ]); + } + + #[test] + fn nullify_upper_triangular_part_examples() { + let mut x = matrix![1.0, 2.0, 3.0; + 4.0, 5.0, 6.0; + 7.0, 8.0, 9.0]; + nullify_upper_triangular_part(&mut x); + assert_matrix_eq!(x, matrix![ + 1.0, 0.0, 0.0; + 4.0, 5.0, 0.0; + 7.0, 8.0, 9.0 + ]); + } +} diff --git a/src/lib.rs b/src/lib.rs index b955c07..da1d4fc 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -100,6 +100,8 @@ pub mod vector; pub mod ulp; pub mod norm; +mod internal_utils; + #[cfg(test)] mod testsupport; diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 81722fa..fa67686 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -1,9 +1,250 @@ use matrix::{Matrix, BaseMatrix}; use error::{Error, ErrorKind}; +use matrix::decomposition::Decomposition; +use matrix::forward_substitution; +use vector::Vector; +use utils::dot; use std::any::Any; -use libnum::Float; +use libnum::{Zero, Float}; + +/// Cholesky decomposition. +/// +/// Given a square, symmetric positive definite matrix A, +/// there exists an invertible lower triangular matrix L +/// such that +/// +/// A = L LT. +/// +/// This is called the Cholesky decomposition of A. +/// For not too ill-conditioned A, the computation +/// of the decomposition is very robust, and it takes about +/// half the effort of an LU decomposition with partial pivoting. +/// +/// # Applications +/// The Cholesky decomposition can be thought of as a specialized +/// LU decomposition for symmetric positive definite matrices, +/// and so its applications are similar to that of LU. +/// +/// The following example shows how to compute the Cholesky +/// decomposition of a given matrix. In this example, we also +/// unpack the decomposition to retrieve the L matrix, +/// but in many practical applications we are not so concerned +/// with the factor itself. Instead, we may wish to +/// solve linear systems or compute the determinant or the +/// inverse of a symmetric positive definite matrix. +/// In this case, see the next subsections. +/// +/// ``` +/// # #[macro_use] extern crate rulinalg; fn main() { +/// use rulinalg::matrix::decomposition::Cholesky; +/// +/// // Need to import Decomposition if we want to unpack +/// use rulinalg::matrix::decomposition::Decomposition; +/// +/// let x = matrix![ 1.0, 3.0, 1.0; +/// 3.0, 13.0, 11.0; +/// 1.0, 11.0, 21.0 ]; +/// let cholesky = Cholesky::decompose(x) +/// .expect("Matrix is SPD."); +/// +/// // Obtain the matrix factor L +/// let l = cholesky.unpack(); +/// +/// assert_matrix_eq!(l, matrix![1.0, 0.0, 0.0; +/// 3.0, 2.0, 0.0; +/// 1.0, 4.0, 2.0], comp = float); +/// # } +/// ``` +/// +/// ## Solving linear systems +/// After having decomposed the matrix, one may efficiently +/// solve linear systems for different right-hand sides. +/// +/// ``` +/// # #[macro_use] extern crate rulinalg; fn main() { +/// # use rulinalg::matrix::decomposition::Cholesky; +/// # let x = matrix![ 1.0, 3.0, 1.0; +/// # 3.0, 13.0, 11.0; +/// # 1.0, 11.0, 21.0 ]; +/// # let cholesky = Cholesky::decompose(x).unwrap(); +/// let b1 = vector![ 3.0, 2.0, 1.0]; +/// let b2 = vector![-2.0, 1.0, 0.0]; +/// let y1 = cholesky.solve(b1).expect("Matrix is invertible."); +/// let y2 = cholesky.solve(b2).expect("Matrix is invertible."); +/// assert_vector_eq!(y1, vector![ 23.25, -7.75, 3.0 ]); +/// assert_vector_eq!(y2, vector![-22.25, 7.75, -3.00 ]); +/// # } +/// ``` +/// +/// ## Computing the inverse of a matrix +/// +/// While computing the inverse explicitly is rarely +/// the best solution to any given problem, it is sometimes +/// necessary. In this case, it is easily accessible +/// through the `inverse()` method on `Cholesky`. +/// +/// # Computing the determinant of a matrix +/// +/// As with LU decomposition, the `Cholesky` decomposition +/// exposes a method `det` for computing the determinant +/// of the decomposed matrix. This is a very cheap operation. +#[derive(Clone, Debug)] +pub struct Cholesky { + l: Matrix +} + +impl Cholesky where T: 'static + Float { + /// Computes the Cholesky decomposition A = L LT + /// for the given square, symmetric positive definite matrix. + /// + /// Note that the implementation cannot reliably and efficiently + /// verify that the matrix truly is symmetric positive definite matrix, + /// so it is the responsibility of the user to make sure that this is + /// the case. In particular, if the input matrix is not SPD, + /// the returned decomposition may not be a valid decomposition + /// for the input matrix. + /// + /// # Errors + /// - A diagonal entry is effectively zero to working precision. + /// - A diagonal entry is negative. + /// + /// # Panics + /// + /// - The matrix must be square. + pub fn decompose(matrix: Matrix) -> Result { + assert!(matrix.rows() == matrix.cols(), + "Matrix must be square for Cholesky decomposition."); + let n = matrix.rows(); + + // The implementation here is based on the + // "Gaxpy-Rich Cholesky Factorization" + // from Chapter 4.2.5 in + // Matrix Computations, 4th Edition, + // (Golub and Van Loan). + + // We consume the matrix we're given, and overwrite its + // lower diagonal part with the L factor. However, + // we ignore the strictly upper triangular part of the matrix, + // because this saves us a few operations. + // When the decomposition is unpacked, we will completely zero + // the upper triangular part. + let mut a = matrix; + + for j in 0 .. n { + if j > 0 { + // This is essentially a GAXPY operation y = y - Bx + // where B is the [j .. n, 0 .. j] submatrix of A, + // x is the [ j, 0 .. j ] submatrix of A, + // and y is the [ j .. n, j ] submatrix of A + for k in j .. n { + let kj_dot = { + let j_row = a.row(j).raw_slice(); + let k_row = a.row(k).raw_slice(); + dot(&k_row[0 .. j], &j_row[0 .. j]) + }; + a[[k, j]] = a[[k, j]] - kj_dot; + } + } + + let diagonal = a[[j, j]]; + if diagonal.abs() < T::epsilon() { + return Err(Error::new(ErrorKind::DecompFailure, + "Matrix is singular to working precision.")); + } else if diagonal < T::zero() { + return Err(Error::new(ErrorKind::DecompFailure, + "Diagonal entries of matrix are not all positive.")); + } + + let divisor = diagonal.sqrt(); + for k in j .. n { + a[[k, j]] = a[[k, j]] / divisor; + } + } + + Ok(Cholesky { + l: a + }) + } + + /// Computes the determinant of the decomposed matrix. + /// + /// Note that the determinant of an empty matrix is considered + /// to be equal to 1. + pub fn det(&self) -> T { + let l_det = self.l.diag() + .cloned() + .fold(T::one(), |a, b| a * b); + l_det * l_det + } + + /// Solves the linear system Ax = b. + /// + /// Here A is the decomposed matrix and b is the + /// supplied vector. + /// + /// # Errors + /// If the matrix is sufficiently ill-conditioned, + /// it is possible that the solution cannot be obtained. + /// + /// # Panics + /// - The supplied right-hand side vector must be + /// dimensionally compatible with the supplied matrix. + pub fn solve(&self, b: Vector) -> Result, Error> { + assert!(self.l.rows() == b.size(), + "RHS vector and coefficient matrix must be + dimensionally compatible."); + // Solve Ly = b + let y = forward_substitution(&self.l, b)?; + // Solve L^T x = y + transpose_back_substitution(&self.l, y) + } + + /// Computes the inverse of the decomposed matrix. + /// + /// # Errors + /// If the matrix is sufficiently ill-conditioned, + /// it is possible that the inverse cannot be obtained. + pub fn inverse(&self) -> Result, Error> { + let n = self.l.rows(); + let mut inv = Matrix::zeros(n, n); + let mut e = Vector::zeros(n); + + // Note: this is essentially the same as + // PartialPivLu::inverse(), and consequently + // the data access patterns here can also be + // improved by way of using BLAS-3 calls. + // Please see that function's implementation + // for more details. + + // Solve for each column of the inverse matrix + for i in 0 .. n { + e[i] = T::one(); + let col = self.solve(e)?; + + for j in 0 .. n { + inv[[j, i]] = col[j]; + } + + e = col.apply(&|_| T::zero()); + } + + Ok(inv) + } +} + +impl Decomposition for Cholesky { + type Factors = Matrix; + + fn unpack(self) -> Matrix { + use internal_utils::nullify_upper_triangular_part; + let mut l = self.l; + nullify_upper_triangular_part(&mut l); + l + } +} + impl Matrix where T: Any + Float @@ -12,6 +253,11 @@ impl Matrix /// /// Returns the cholesky decomposition of a positive definite matrix. /// + /// *NOTE*: This function is deprecated, and will be removed in a + /// future release. Please see + /// [Cholesky](decomposition/struct.Cholesky.html) for its + /// replacement. + /// /// # Examples /// /// ``` @@ -33,6 +279,7 @@ impl Matrix /// # Failures /// /// - Matrix is not positive definite. + #[deprecated] pub fn cholesky(&self) -> Result, Error> { assert!(self.rows == self.cols, "Matrix must be square for Cholesky decomposition."); @@ -77,15 +324,274 @@ impl Matrix } } +/// Solves the square system L^T x = b, +/// where L is lower triangular +fn transpose_back_substitution(l: &Matrix, b: Vector) + -> Result, Error> where T: Float { + assert!(l.rows() == l.cols(), "Matrix L must be square."); + assert!(l.rows() == b.size(), "L and b must be dimensionally compatible."); + let n = l.rows(); + let mut x = b; + + for i in (0 .. n).rev() { + let row = l.row(i).raw_slice(); + let diagonal = l[[i, i]]; + if diagonal.abs() < T::epsilon() { + return Err(Error::new(ErrorKind::DivByZero, + "Matrix L is singular to working precision.")); + } + + x[i] = x[i] / diagonal; + + // Apply the BLAS-1 operation + // y <- y + α x + // where α = - x[i], + // y = x[0 .. i] + // and x = l[i, 0 .. i] + // TODO: Hopefully we'll have a more systematic way + // of applying optimized BLAS-like operations in the future. + // In this case, we should replace this loop with a call + // to the appropriate function. + for j in 0 .. i { + x[j] = x[j] - x[i] * row[j]; + } + } + + Ok(x) +} + #[cfg(test)] mod tests { use matrix::Matrix; + use matrix::decomposition::Decomposition; + use vector::Vector; + + use super::Cholesky; + use super::transpose_back_substitution; + + use libnum::Float; + use quickcheck::TestResult; #[test] #[should_panic] + #[allow(deprecated)] fn test_non_square_cholesky() { let a = Matrix::::ones(2, 3); let _ = a.cholesky(); } + + #[test] + fn cholesky_unpack_empty() { + let x: Matrix = matrix![]; + let l = Cholesky::decompose(x.clone()) + .unwrap() + .unpack(); + assert_matrix_eq!(l, x); + } + + #[test] + fn cholesky_unpack_1x1() { + let x = matrix![ 4.0 ]; + let expected = matrix![ 2.0 ]; + let l = Cholesky::decompose(x) + .unwrap() + .unpack(); + assert_matrix_eq!(l, expected, comp = float); + } + + #[test] + fn cholesky_unpack_2x2() { + { + let x = matrix![ 9.0, -6.0; + -6.0, 20.0]; + let expected = matrix![ 3.0, 0.0; + -2.0, 4.0]; + + let l = Cholesky::decompose(x) + .unwrap() + .unpack(); + assert_matrix_eq!(l, expected, comp = float); + } + } + + #[test] + fn cholesky_singular_fails() { + { + let x = matrix![0.0]; + assert!(Cholesky::decompose(x).is_err()); + } + + { + let x = matrix![0.0, 0.0; + 0.0, 1.0]; + assert!(Cholesky::decompose(x).is_err()); + } + + { + let x = matrix![1.0, 0.0; + 0.0, 0.0]; + assert!(Cholesky::decompose(x).is_err()); + } + + { + let x = matrix![1.0, 3.0, 5.0; + 3.0, 9.0, 15.0; + 5.0, 15.0, 65.0]; + assert!(Cholesky::decompose(x).is_err()); + } + } + + #[test] + fn cholesky_det_empty() { + let x: Matrix = matrix![]; + let cholesky = Cholesky::decompose(x).unwrap(); + assert_eq!(cholesky.det(), 1.0); + } + + #[test] + fn cholesky_det() { + { + let x = matrix![1.0]; + let cholesky = Cholesky::decompose(x).unwrap(); + let diff = cholesky.det() - 1.0; + assert!(diff.abs() < 1e-14); + } + + { + let x = matrix![1.0, 3.0, 5.0; + 3.0, 18.0, 33.0; + 5.0, 33.0, 65.0]; + let cholesky = Cholesky::decompose(x).unwrap(); + let diff = cholesky.det() - 36.0; + assert!(diff.abs() < 1e-14); + } + } + + #[test] + fn cholesky_solve_examples() { + { + let a: Matrix = matrix![]; + let b: Vector = vector![]; + let expected: Vector = vector![]; + let cholesky = Cholesky::decompose(a).unwrap(); + let x = cholesky.solve(b).unwrap(); + assert_eq!(x, expected); + } + + { + let a = matrix![ 1.0 ]; + let b = vector![ 4.0 ]; + let expected = vector![ 4.0 ]; + let cholesky = Cholesky::decompose(a).unwrap(); + let x = cholesky.solve(b).unwrap(); + assert_vector_eq!(x, expected, comp = float); + } + + { + let a = matrix![ 4.0, 6.0; + 6.0, 25.0]; + let b = vector![ 2.0, 4.0]; + let expected = vector![ 0.40625, 0.0625 ]; + let cholesky = Cholesky::decompose(a).unwrap(); + let x = cholesky.solve(b).unwrap(); + assert_vector_eq!(x, expected, comp = float); + } + } + + #[test] + fn cholesky_inverse_examples() { + { + let a: Matrix = matrix![]; + let expected: Matrix = matrix![]; + let cholesky = Cholesky::decompose(a).unwrap(); + assert_eq!(cholesky.inverse().unwrap(), expected); + } + + { + let a = matrix![ 2.0 ]; + let expected = matrix![ 0.5 ]; + let cholesky = Cholesky::decompose(a).unwrap(); + assert_matrix_eq!(cholesky.inverse().unwrap(), expected, + comp = float); + } + + { + let a = matrix![ 4.0, 6.0; + 6.0, 25.0]; + let expected = matrix![ 0.390625, -0.09375; + -0.093750 , 0.06250]; + let cholesky = Cholesky::decompose(a).unwrap(); + assert_matrix_eq!(cholesky.inverse().unwrap(), expected, + comp = float); + } + + { + let a = matrix![ 9.0, 6.0, 3.0; + 6.0, 20.0, 10.0; + 3.0, 10.0, 14.0]; + let expected = matrix![0.1388888888888889, -0.0416666666666667, 0.0 ; + -0.0416666666666667, 0.0902777777777778, -0.0555555555555556; + 0.0, -0.0555555555555556, 0.1111111111111111]; + let cholesky = Cholesky::decompose(a).unwrap(); + assert_matrix_eq!(cholesky.inverse().unwrap(), expected, + comp = float); + } + } + + quickcheck! { + fn property_cholesky_of_identity_is_identity(n: usize) -> TestResult { + if n > 30 { + return TestResult::discard(); + } + + let x = Matrix::::identity(n); + let l = Cholesky::decompose(x.clone()).map(|c| c.unpack()); + match l { + Ok(l) => { + assert_matrix_eq!(l, x, comp = float); + TestResult::passed() + }, + _ => TestResult::failed() + } + } + } + + #[test] + fn transpose_back_substitution_examples() { + { + let l: Matrix = matrix![]; + let b: Vector = vector![]; + let expected: Vector = vector![]; + let x = transpose_back_substitution(&l, b).unwrap(); + assert_vector_eq!(x, expected); + } + + { + let l = matrix![2.0]; + let b = vector![2.0]; + let expected = vector![1.0]; + let x = transpose_back_substitution(&l, b).unwrap(); + assert_vector_eq!(x, expected, comp = float); + } + + { + let l = matrix![2.0, 0.0; + 3.0, 4.0]; + let b = vector![2.0, 1.0]; + let expected = vector![0.625, 0.25 ]; + let x = transpose_back_substitution(&l, b).unwrap(); + assert_vector_eq!(x, expected, comp = float); + } + + { + let l = matrix![ 2.0, 0.0, 0.0; + 5.0, -1.0, 0.0; + -2.0, 0.0, 1.0]; + let b = vector![-1.0, 2.0, 3.0]; + let expected = vector![ 7.5, -2.0, 3.0 ]; + let x = transpose_back_substitution(&l, b).unwrap(); + assert_vector_eq!(x, expected, comp = float); + } + } } diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index a504516..7913935 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -136,8 +136,10 @@ 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 u = nullify_lower_triangular_part(self.lu); + let mut u = self.lu; + nullify_lower_triangular_part(&mut u); LUP { l: l, @@ -283,6 +285,9 @@ impl PartialPivLu where T: Any + Float { } /// Computes the determinant of the decomposed matrix. + /// + /// Note that the determinant of an empty matrix is considered + /// to be equal to 1. pub fn det(&self) -> T { // Recall that the determinant of a triangular matrix // is the product of its diagonal entries. Also, @@ -341,8 +346,10 @@ impl Decomposition for FullPivLu { type Factors = LUPQ; fn unpack(self) -> LUPQ { + use internal_utils::nullify_lower_triangular_part; let l = unit_lower_triangular_part(&self.lu); - let u = nullify_lower_triangular_part(self.lu); + let mut u = self.lu; + nullify_lower_triangular_part(&mut u); LUPQ { l: l, @@ -655,15 +662,6 @@ fn unit_lower_triangular_part(matrix: &M) -> Matrix Matrix::new(m, m, 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 { diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index 66efe8f..54aea23 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -62,14 +62,14 @@ //! //! //! Decomposition -//! Applicable to +//! Matrix requirements //! Supported features //! //! //! //! //! PartialPivLu -//! Square, invertible matrices +//! Square, invertible //! //!
    //!
  • Linear system solving
    • @@ -92,6 +92,18 @@ //! //! //! +//! +//! Cholesky +//! Square, symmetric positive definite +//! +//!
      +//!
    • Linear system solving
      • +//!
    • Matrix inverse
      • +//!
    • Determinant computation
      • +//!
    +//! +//! +//! //! //! @@ -123,6 +135,7 @@ use utils; use error::{Error, ErrorKind}; pub use self::lu::{PartialPivLu, LUP, FullPivLu, LUPQ}; +pub use self::cholesky::Cholesky; use libnum::{Float}; diff --git a/tests/mat/mod.rs b/tests/mat/mod.rs index 3cc45ac..239aa96 100644 --- a/tests/mat/mod.rs +++ b/tests/mat/mod.rs @@ -171,6 +171,7 @@ fn matrix_partial_piv_lu() { } #[test] +#[allow(deprecated)] fn cholesky() { let a = matrix![25., 15., -5.; 15., 18., 0.;