From e1a1a0efd70fcb122bd3e35a016cc0aaa626ae47 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 11 Feb 2017 18:31:11 +0100 Subject: [PATCH 01/24] Add minimal internal_utils module --- src/internal_utils.rs | 52 +++++++++++++++++++++++++++++++++++++++++++ src/lib.rs | 2 ++ 2 files changed, 54 insertions(+) create mode 100644 src/internal_utils.rs diff --git a/src/internal_utils.rs b/src/internal_utils.rs new file mode 100644 index 0000000..5943985 --- /dev/null +++ b/src/internal_utils.rs @@ -0,0 +1,52 @@ +use matrix::BaseMatrixMut; +use libnum::Zero; + +pub fn nullify_lower_triangular_part<T, M>(matrix: &mut M) + where T: Zero, M: BaseMatrixMut<T> { + for (i, mut row) in matrix.row_iter_mut().enumerate() { + for element in row.iter_mut().take(i) { + *element = T::zero(); + } + } +} + +pub fn nullify_upper_triangular_part<T, M>(matrix: &mut M) + where T: Zero, M: BaseMatrixMut<T> { + for (i, mut row) in matrix.row_iter_mut().enumerate() { + for element in row.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; From 5f04dcb786991e853f1f7e6d8e62ac1850841648 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 11 Feb 2017 18:32:08 +0100 Subject: [PATCH 02/24] Basic implementation of Cholesky --- src/matrix/decomposition/cholesky.rs | 121 ++++++++++++++++++++++++++- src/matrix/decomposition/mod.rs | 1 + 2 files changed, 121 insertions(+), 1 deletion(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 81722fa..43a5781 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -1,9 +1,72 @@ use matrix::{Matrix, BaseMatrix}; use error::{Error, ErrorKind}; +use matrix::decomposition::Decomposition; use std::any::Any; -use libnum::Float; +use libnum::{Zero, Float}; + +/// TODO +#[derive(Clone, Debug)] +pub struct Cholesky<T> { + l: Matrix<T> +} + +impl<T> Cholesky<T> where T: Float { + /// TODO + fn decompose(matrix: Matrix<T>) -> Result<Self, Error> { + 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; + + // Resolve each submatrix (j .. n, j .. n) + for j in 0 .. n { + if j > 0 { + for k in j .. n { + for l in 0 .. j { + a[[k, j]] = a[[k, j]] - a[[k, l]] * a[[j, l]]; + } + } + } + + // TODO: Check for zero divisor + let divisor = a[[j, j]].sqrt(); + for k in j .. n { + a[[k, j]] = a[[k, j]] / divisor; + } + } + + Ok(Cholesky { + l: a + }) + } +} + +impl<T: Zero> Decomposition for Cholesky<T> { + type Factors = Matrix<T>; + + fn unpack(self) -> Matrix<T> { + use internal_utils::nullify_upper_triangular_part; + let mut l = self.l; + nullify_upper_triangular_part(&mut l); + l + } +} + impl<T> Matrix<T> where T: Any + Float @@ -80,6 +143,10 @@ impl<T> Matrix<T> #[cfg(test)] mod tests { use matrix::Matrix; + use matrix::decomposition::Decomposition; + use super::Cholesky; + + use quickcheck::TestResult; #[test] #[should_panic] @@ -88,4 +155,56 @@ mod tests { let _ = a.cholesky(); } + + #[test] + fn cholesky_unpack_empty() { + let x: Matrix<f64> = 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); + } + } + + quickcheck! { + fn property_cholesky_of_identity_is_identity(n: usize) -> TestResult { + if n > 30 { + return TestResult::discard(); + } + + let x = Matrix::<f64>::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() + } + } + } } diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index afa0e1c..adad954 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -110,6 +110,7 @@ use utils; use error::{Error, ErrorKind}; pub use self::lu::{PartialPivLu, LUP}; +pub use self::cholesky::Cholesky; use libnum::{Float}; From 754de355851b08eca7f6983a8ddc88a8ef2a6894 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 11 Feb 2017 18:36:46 +0100 Subject: [PATCH 03/24] Rewrite lu unpack to use nullify from internal_utils --- src/matrix/decomposition/lu.rs | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 17fadd8..ce81c7f 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -135,8 +135,10 @@ impl<T: Clone + One + Zero> Decomposition for PartialPivLu<T> { type Factors = LUP<T>; fn unpack(self) -> LUP<T> { + 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, @@ -347,15 +349,6 @@ fn unit_lower_triangular_part<T, M>(matrix: &M) -> Matrix<T> Matrix::new(m, n, data) } -fn nullify_lower_triangular_part<T: Zero>(mut matrix: Matrix<T>) -> Matrix<T> { - for (i, mut row) in matrix.row_iter_mut().enumerate() { - for element in row.iter_mut().take(i) { - *element = T::zero(); - } - } - matrix -} - impl<T> Matrix<T> where T: Any + Float { From f3807c78c4c34f0efb4a61737b410556b0de14ef Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 11 Feb 2017 21:04:22 +0100 Subject: [PATCH 04/24] Add Cholesky benchmarks --- benches/lib.rs | 1 + benches/linalg/cholesky.rs | 49 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 50 insertions(+) create mode 100644 benches/linalg/cholesky.rs 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..4ece138 --- /dev/null +++ b/benches/linalg/cholesky.rs @@ -0,0 +1,49 @@ +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::<f64>::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::<f64>::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::<f64>::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::<f64>::identity(n); + b.iter(|| { + x.cholesky().expect("Matrix is invertible") + }) +} From 38b08a373fea7a2bd4ad28092022d9e9fd05c7d0 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 11 Feb 2017 21:04:52 +0100 Subject: [PATCH 05/24] Take raw_slice_mut() in nullify --- src/internal_utils.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/internal_utils.rs b/src/internal_utils.rs index 5943985..a1faaaf 100644 --- a/src/internal_utils.rs +++ b/src/internal_utils.rs @@ -4,7 +4,7 @@ use libnum::Zero; pub fn nullify_lower_triangular_part<T, M>(matrix: &mut M) where T: Zero, M: BaseMatrixMut<T> { for (i, mut row) in matrix.row_iter_mut().enumerate() { - for element in row.iter_mut().take(i) { + for element in row.raw_slice_mut().iter_mut().take(i) { *element = T::zero(); } } @@ -13,7 +13,7 @@ pub fn nullify_lower_triangular_part<T, M>(matrix: &mut M) pub fn nullify_upper_triangular_part<T, M>(matrix: &mut M) where T: Zero, M: BaseMatrixMut<T> { for (i, mut row) in matrix.row_iter_mut().enumerate() { - for element in row.iter_mut().skip(i + 1) { + for element in row.raw_slice_mut().iter_mut().skip(i + 1) { *element = T::zero(); } } From 29163f4151eecb44f77fb866ff13a0cc6ec163c4 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 11 Feb 2017 21:08:48 +0100 Subject: [PATCH 06/24] Rewrite Cholesky in terms of dot() Benchmarks show significant performance gains. --- src/matrix/decomposition/cholesky.rs | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 43a5781..60929d9 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -1,6 +1,7 @@ use matrix::{Matrix, BaseMatrix}; use error::{Error, ErrorKind}; use matrix::decomposition::Decomposition; +use utils::dot; use std::any::Any; @@ -14,7 +15,7 @@ pub struct Cholesky<T> { impl<T> Cholesky<T> where T: Float { /// TODO - fn decompose(matrix: Matrix<T>) -> Result<Self, Error> { + pub fn decompose(matrix: Matrix<T>) -> Result<Self, Error> { assert!(matrix.rows() == matrix.cols(), "Matrix must be square for Cholesky decomposition."); let n = matrix.rows(); @@ -36,10 +37,17 @@ impl<T> Cholesky<T> where T: Float { // Resolve each submatrix (j .. n, j .. n) 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 { - for l in 0 .. j { - a[[k, j]] = a[[k, j]] - a[[k, l]] * a[[j, l]]; - } + 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; } } From 4520bd56d3e1415cea31a2314055202221fb0752 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Mon, 13 Feb 2017 16:12:39 +0100 Subject: [PATCH 07/24] Cholesky: validate diagonal entries --- src/matrix/decomposition/cholesky.rs | 39 ++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 60929d9..1566181 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -51,8 +51,16 @@ impl<T> Cholesky<T> where T: Float { } } - // TODO: Check for zero divisor - let divisor = a[[j, j]].sqrt(); + 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; } @@ -198,6 +206,33 @@ mod tests { } } + #[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()); + } + } + quickcheck! { fn property_cholesky_of_identity_is_identity(n: usize) -> TestResult { if n > 30 { From 91b075127cee57329cdf02e19c8ce50e16075953 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Mon, 13 Feb 2017 20:03:35 +0100 Subject: [PATCH 08/24] Implement det() for Cholesky --- src/matrix/decomposition/cholesky.rs | 35 ++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 1566181..178be9c 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -70,6 +70,14 @@ impl<T> Cholesky<T> where T: Float { l: a }) } + + /// TODO + pub fn det(&self) -> T { + let l_det = self.l.diag() + .cloned() + .fold(T::one(), |a, b| a * b); + l_det * l_det + } } impl<T: Zero> Decomposition for Cholesky<T> { @@ -162,6 +170,7 @@ mod tests { use matrix::decomposition::Decomposition; use super::Cholesky; + use libnum::Float; use quickcheck::TestResult; #[test] @@ -233,6 +242,32 @@ mod tests { } } + #[test] + fn cholesky_det_empty() { + let x: Matrix<f64> = 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); + } + } + quickcheck! { fn property_cholesky_of_identity_is_identity(n: usize) -> TestResult { if n > 30 { From 038f129a537ba74cf409318ced0e68126e4fcc4c Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Tue, 14 Feb 2017 15:27:25 +0100 Subject: [PATCH 09/24] Implement transpose_back_substitution --- src/matrix/decomposition/cholesky.rs | 72 ++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 178be9c..9c22f1d 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -1,6 +1,7 @@ use matrix::{Matrix, BaseMatrix}; use error::{Error, ErrorKind}; use matrix::decomposition::Decomposition; +use vector::Vector; use utils::dot; use std::any::Any; @@ -164,11 +165,44 @@ impl<T> Matrix<T> } } +/// Solves the square system L^T x = b, +/// where L is lower triangular +fn transpose_back_substitution<T>(l: &Matrix<T>, b: Vector<T>) + -> Result<Vector<T>, 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; + + // TODO: Make this implementation more cache efficient + // At the moment it is a simple naive (very cache inefficient) + // implementation for the sake of correctness. + for i in (0 .. n).rev() { + let mut inner_product = T::zero(); + for j in (i + 1) .. n { + inner_product = inner_product + l[[j, i]] * x[j]; + } + + 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] - inner_product) / diagonal; + } + + 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; @@ -285,4 +319,42 @@ mod tests { } } } + + #[test] + fn transpose_back_substitution_examples() { + { + let l: Matrix<f64> = matrix![]; + let b: Vector<f64> = vector![]; + let expected: Vector<f64> = 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); + } + } } From d24331376ad905917fb58a22bb55b609e7b68b0c Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Tue, 14 Feb 2017 15:58:54 +0100 Subject: [PATCH 10/24] Implement Cholesky::solve --- src/matrix/decomposition/cholesky.rs | 52 +++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 9c22f1d..4e9efa1 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -1,6 +1,7 @@ use matrix::{Matrix, BaseMatrix}; use error::{Error, ErrorKind}; use matrix::decomposition::Decomposition; +use matrix::forward_substitution; use vector::Vector; use utils::dot; @@ -14,7 +15,7 @@ pub struct Cholesky<T> { l: Matrix<T> } -impl<T> Cholesky<T> where T: Float { +impl<T> Cholesky<T> where T: 'static + Float { /// TODO pub fn decompose(matrix: Matrix<T>) -> Result<Self, Error> { assert!(matrix.rows() == matrix.cols(), @@ -79,6 +80,19 @@ impl<T> Cholesky<T> where T: Float { .fold(T::one(), |a, b| a * b); l_det * l_det } + + /// TODO + pub fn solve(&self, b: Vector<T>) -> Vector<T> { + 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) + .expect("Internal error: L should be invertible."); + // Solve L^T x = y + transpose_back_substitution(&self.l, y) + .expect("Internal error: L^T should be invertible.") + } } impl<T: Zero> Decomposition for Cholesky<T> { @@ -302,6 +316,42 @@ mod tests { } } + #[test] + fn cholesky_solve_examples() { + { + // TODO: Enable this test + // It is currently commented out because + // backward/forward substitution don't handle + // empty matrices. See PR #152 for more details. + + // let a: Matrix<f64> = matrix![]; + // let b: Vector<f64> = vector![]; + // let expected: Vector<f64> = vector![]; + // let cholesky = Cholesky::decompose(a).unwrap(); + // let x = cholesky.solve(b); + // 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); + 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); + assert_vector_eq!(x, expected, comp = float); + } + } + quickcheck! { fn property_cholesky_of_identity_is_identity(n: usize) -> TestResult { if n > 30 { From 8b398a1df20a9c6510acf15f8fa3401793acf125 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Tue, 14 Feb 2017 19:07:02 +0100 Subject: [PATCH 11/24] Add benchmarks for Cholesky::solve --- benches/linalg/cholesky.rs | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/benches/linalg/cholesky.rs b/benches/linalg/cholesky.rs index 4ece138..f7b6797 100644 --- a/benches/linalg/cholesky.rs +++ b/benches/linalg/cholesky.rs @@ -47,3 +47,23 @@ fn cholesky_500x500(b: &mut Bencher) { 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]) + }); +} From 5b47cd524aaa1b2c7e210bfd16a09e670356113c Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Mon, 20 Feb 2017 12:04:46 +0100 Subject: [PATCH 12/24] Rewrite transpose_back_sub.. for better data access pattern --- src/matrix/decomposition/cholesky.rs | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 4e9efa1..368e144 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -188,22 +188,28 @@ fn transpose_back_substitution<T>(l: &Matrix<T>, b: Vector<T>) let n = l.rows(); let mut x = b; - // TODO: Make this implementation more cache efficient - // At the moment it is a simple naive (very cache inefficient) - // implementation for the sake of correctness. for i in (0 .. n).rev() { - let mut inner_product = T::zero(); - for j in (i + 1) .. n { - inner_product = inner_product + l[[j, i]] * x[j]; - } - + 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] - inner_product) / diagonal; + 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) From 5c8a566c7be8289a0138215cffbc8d1de0bd846d Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Mon, 20 Feb 2017 12:11:52 +0100 Subject: [PATCH 13/24] Enable previously disabled test --- src/matrix/decomposition/cholesky.rs | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 368e144..5c351ea 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -325,17 +325,12 @@ mod tests { #[test] fn cholesky_solve_examples() { { - // TODO: Enable this test - // It is currently commented out because - // backward/forward substitution don't handle - // empty matrices. See PR #152 for more details. - - // let a: Matrix<f64> = matrix![]; - // let b: Vector<f64> = vector![]; - // let expected: Vector<f64> = vector![]; - // let cholesky = Cholesky::decompose(a).unwrap(); - // let x = cholesky.solve(b); - // assert_eq!(x, expected); + let a: Matrix<f64> = matrix![]; + let b: Vector<f64> = vector![]; + let expected: Vector<f64> = vector![]; + let cholesky = Cholesky::decompose(a).unwrap(); + let x = cholesky.solve(b); + assert_eq!(x, expected); } { From 2c45b7b7fa94b7ecf3125819e60fed67bcdbeca1 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Mon, 20 Feb 2017 21:36:04 +0100 Subject: [PATCH 14/24] Implement Cholesky::inverse() --- src/matrix/decomposition/cholesky.rs | 65 ++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 5c351ea..8e376d0 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -93,6 +93,34 @@ impl<T> Cholesky<T> where T: 'static + Float { transpose_back_substitution(&self.l, y) .expect("Internal error: L^T should be invertible.") } + + /// Computes the inverse of the decomposed matrix. + pub fn inverse(&self) -> Matrix<T> { + 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()); + } + + inv + } } impl<T: Zero> Decomposition for Cholesky<T> { @@ -353,6 +381,43 @@ mod tests { } } + #[test] + fn cholesky_inverse_examples() { + { + let a: Matrix<f64> = matrix![]; + let expected: Matrix<f64> = matrix![]; + let cholesky = Cholesky::decompose(a).unwrap(); + assert_eq!(cholesky.inverse(), expected); + } + + { + let a = matrix![ 2.0 ]; + let expected = matrix![ 0.5 ]; + let cholesky = Cholesky::decompose(a).unwrap(); + assert_matrix_eq!(cholesky.inverse(), 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(), 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(), expected, comp = float); + } + } + quickcheck! { fn property_cholesky_of_identity_is_identity(n: usize) -> TestResult { if n > 30 { From 5de2e509ea88a3ed0d24b06665188e84958f3c22 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Mon, 20 Feb 2017 22:20:09 +0100 Subject: [PATCH 15/24] Documentation for Cholesky --- src/matrix/decomposition/cholesky.rs | 104 +++++++++++++++++++++++++-- 1 file changed, 100 insertions(+), 4 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 8e376d0..bbaddb8 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -9,14 +9,103 @@ use std::any::Any; use libnum::{Zero, Float}; -/// TODO +/// Cholesky decomposition. +/// +/// Given a square, symmetric positive definite matrix A, +/// there exists an invertible lower triangular matrix L +/// such that +/// +/// A = L L<sup>T</sup>. +/// +/// 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]; +/// assert_vector_eq!(cholesky.solve(b1), vector![ 22.00, -7.25, 2.75 ]); +/// assert_vector_eq!(cholesky.solve(b2), 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<T> { l: Matrix<T> } impl<T> Cholesky<T> where T: 'static + Float { - /// TODO + /// Computes the Cholesky decomposition A = L L<sup>T</sup> + /// for the given square, symmetric positive definite matrix. + /// + /// # Errors + /// - A diagonal entry is very close to zero, which + /// corresponds to the matrix being effectively singular + /// to working precision. + /// - A diagonal entry is negative. + /// + /// # Panics + /// + /// - The matrix must be square. pub fn decompose(matrix: Matrix<T>) -> Result<Self, Error> { assert!(matrix.rows() == matrix.cols(), "Matrix must be square for Cholesky decomposition."); @@ -73,7 +162,7 @@ impl<T> Cholesky<T> where T: 'static + Float { }) } - /// TODO + /// Computes the determinant of the decomposed matrix. pub fn det(&self) -> T { let l_det = self.l.diag() .cloned() @@ -81,7 +170,14 @@ impl<T> Cholesky<T> where T: 'static + Float { l_det * l_det } - /// TODO + /// Solves the linear system Ax = b. + /// + /// Here A is the decomposed matrix and b is the + /// supplied vector. + /// + /// # Panics + /// - The supplied right-hand side vector must be + /// dimensionally compatible with the supplied matrix. pub fn solve(&self, b: Vector<T>) -> Vector<T> { assert!(self.l.rows() == b.size(), "RHS vector and coefficient matrix must be From a4eb96c7ea7b44a17d845bfaf7fd183848f7f633 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Mon, 20 Feb 2017 22:40:24 +0100 Subject: [PATCH 16/24] Fix errors in Cholesky doctests --- src/matrix/decomposition/cholesky.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index bbaddb8..68aa52f 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -49,7 +49,7 @@ use libnum::{Zero, Float}; /// let cholesky = Cholesky::decompose(x) /// .expect("Matrix is SPD."); /// -/// Obtain the matrix factor L +/// // Obtain the matrix factor L /// let l = cholesky.unpack(); /// /// assert_matrix_eq!(l, matrix![1.0, 0.0, 0.0; @@ -71,7 +71,7 @@ use libnum::{Zero, Float}; /// # let cholesky = Cholesky::decompose(x).unwrap(); /// let b1 = vector![ 3.0, 2.0, 1.0]; /// let b2 = vector![-2.0, 1.0, 0.0]; -/// assert_vector_eq!(cholesky.solve(b1), vector![ 22.00, -7.25, 2.75 ]); +/// assert_vector_eq!(cholesky.solve(b1), vector![ 23.25, -7.75, 3.0 ]); /// assert_vector_eq!(cholesky.solve(b2), vector![-22.25, 7.75, -3.00 ]); /// # } /// ``` From 9e28cdd95dfd6d529888df7f7465e98f836670aa Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 25 Feb 2017 14:43:38 +0100 Subject: [PATCH 17/24] Remove stray /// --- src/matrix/decomposition/cholesky.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 68aa52f..e7fdda9 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -34,7 +34,7 @@ use libnum::{Zero, Float}; /// 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./// +/// In this case, see the next subsections. /// /// ``` /// # #[macro_use] extern crate rulinalg; fn main() { From c222bb1045bc16fa845cec8f16c79d007c986c99 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 25 Feb 2017 14:46:19 +0100 Subject: [PATCH 18/24] Remove confusing comment --- src/matrix/decomposition/cholesky.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index e7fdda9..e4c980f 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -125,7 +125,6 @@ impl<T> Cholesky<T> where T: 'static + Float { // the upper triangular part. let mut a = matrix; - // Resolve each submatrix (j .. n, j .. n) for j in 0 .. n { if j > 0 { // This is essentially a GAXPY operation y = y - Bx From f38705626a7c3ed2ba3aabb1563196c37afeaa50 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 25 Feb 2017 14:55:47 +0100 Subject: [PATCH 19/24] Add Cholesky to decomposition module overview --- src/matrix/decomposition/mod.rs | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index adad954..044774e 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -62,14 +62,26 @@ //! <thead> //! <tr> //! <th>Decomposition</th> -//! <th>Applicable to</th> +//! <th>Matrix requirements</th> //! <th>Supported features</th> //! </tr> //! <tbody> //! //! <tr> //! <td><a href="struct.PartialPivLu.html">PartialPivLu</a></td> -//! <td>Square, invertible matrices</td> +//! <td>Square, invertible</td> +//! <td> +//! <ul> +//! <li>Linear system solving</li> +//! <li>Matrix inverse</li> +//! <li>Determinant computation</li> +//! </ul> +//! </td> +//! </tr> +//! +//! <tr> +//! <td><a href="struct.Cholesky.html">Cholesky</a></td> +//! <td>Square, symmetric positive definite</td> //! <td> //! <ul> //! <li>Linear system solving</li> From 5af62fa44da534ab31f822a0f1d1c11cc8426592 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 25 Feb 2017 15:04:37 +0100 Subject: [PATCH 20/24] Deprecate .cholesky() --- src/matrix/decomposition/cholesky.rs | 7 +++++++ tests/mat/mod.rs | 1 + 2 files changed, 8 insertions(+) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index e4c980f..c38e4b8 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -237,6 +237,11 @@ impl<T> Matrix<T> /// /// 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 /// /// ``` @@ -258,6 +263,7 @@ impl<T> Matrix<T> /// # Failures /// /// - Matrix is not positive definite. + #[deprecated] pub fn cholesky(&self) -> Result<Matrix<T>, Error> { assert!(self.rows == self.cols, "Matrix must be square for Cholesky decomposition."); @@ -352,6 +358,7 @@ mod tests { #[test] #[should_panic] + #[allow(deprecated)] fn test_non_square_cholesky() { let a = Matrix::<f64>::ones(2, 3); 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.; From 04c2d778cf30b8e407978d3967e9a555ba65e300 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Sat, 25 Feb 2017 15:11:55 +0100 Subject: [PATCH 21/24] Add warning to Cholesky::decompose --- src/matrix/decomposition/cholesky.rs | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index c38e4b8..24344d1 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -97,10 +97,15 @@ impl<T> Cholesky<T> where T: 'static + Float { /// Computes the Cholesky decomposition A = L L<sup>T</sup> /// 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 very close to zero, which - /// corresponds to the matrix being effectively singular - /// to working precision. + /// - A diagonal entry is effectively zero to working precision. /// - A diagonal entry is negative. /// /// # Panics From 1916c71d077a22a7b2ac8337928279cf94a40f0c Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Thu, 23 Mar 2017 14:53:57 +0100 Subject: [PATCH 22/24] Rewrite Cholesky solve/inverse to return Result --- src/matrix/decomposition/cholesky.rs | 43 +++++++++++++++++----------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 24344d1..5dc5218 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -71,8 +71,10 @@ use libnum::{Zero, Float}; /// # let cholesky = Cholesky::decompose(x).unwrap(); /// let b1 = vector![ 3.0, 2.0, 1.0]; /// let b2 = vector![-2.0, 1.0, 0.0]; -/// assert_vector_eq!(cholesky.solve(b1), vector![ 23.25, -7.75, 3.0 ]); -/// assert_vector_eq!(cholesky.solve(b2), vector![-22.25, 7.75, -3.00 ]); +/// 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 ]); /// # } /// ``` /// @@ -179,23 +181,29 @@ impl<T> Cholesky<T> where T: 'static + Float { /// 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<T>) -> Vector<T> { + pub fn solve(&self, b: Vector<T>) -> Result<Vector<T>, 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) - .expect("Internal error: L should be invertible."); + let y = forward_substitution(&self.l, b)?; // Solve L^T x = y transpose_back_substitution(&self.l, y) - .expect("Internal error: L^T should be invertible.") } /// Computes the inverse of the decomposed matrix. - pub fn inverse(&self) -> Matrix<T> { + /// + /// # Errors + /// If the matrix is sufficiently ill-conditioned, + /// it is possible that the inverse cannot be obtained. + pub fn inverse(&self) -> Result<Matrix<T>, Error> { let n = self.l.rows(); let mut inv = Matrix::zeros(n, n); let mut e = Vector::zeros(n); @@ -210,7 +218,7 @@ impl<T> Cholesky<T> where T: 'static + Float { // Solve for each column of the inverse matrix for i in 0 .. n { e[i] = T::one(); - let col = self.solve(e); + let col = self.solve(e)?; for j in 0 .. n { inv[[j, i]] = col[j]; @@ -219,7 +227,7 @@ impl<T> Cholesky<T> where T: 'static + Float { e = col.apply(&|_| T::zero()); } - inv + Ok(inv) } } @@ -464,7 +472,7 @@ mod tests { let b: Vector<f64> = vector![]; let expected: Vector<f64> = vector![]; let cholesky = Cholesky::decompose(a).unwrap(); - let x = cholesky.solve(b); + let x = cholesky.solve(b).unwrap(); assert_eq!(x, expected); } @@ -473,7 +481,7 @@ mod tests { let b = vector![ 4.0 ]; let expected = vector![ 4.0 ]; let cholesky = Cholesky::decompose(a).unwrap(); - let x = cholesky.solve(b); + let x = cholesky.solve(b).unwrap(); assert_vector_eq!(x, expected, comp = float); } @@ -483,7 +491,7 @@ mod tests { 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); + let x = cholesky.solve(b).unwrap(); assert_vector_eq!(x, expected, comp = float); } } @@ -494,14 +502,15 @@ mod tests { let a: Matrix<f64> = matrix![]; let expected: Matrix<f64> = matrix![]; let cholesky = Cholesky::decompose(a).unwrap(); - assert_eq!(cholesky.inverse(), expected); + 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(), expected, comp = float); + assert_matrix_eq!(cholesky.inverse().unwrap(), expected, + comp = float); } { @@ -510,7 +519,8 @@ mod tests { let expected = matrix![ 0.390625, -0.09375; -0.093750 , 0.06250]; let cholesky = Cholesky::decompose(a).unwrap(); - assert_matrix_eq!(cholesky.inverse(), expected, comp = float); + assert_matrix_eq!(cholesky.inverse().unwrap(), expected, + comp = float); } { @@ -521,7 +531,8 @@ mod tests { -0.0416666666666667, 0.0902777777777778, -0.0555555555555556; 0.0, -0.0555555555555556, 0.1111111111111111]; let cholesky = Cholesky::decompose(a).unwrap(); - assert_matrix_eq!(cholesky.inverse(), expected, comp = float); + assert_matrix_eq!(cholesky.inverse().unwrap(), expected, + comp = float); } } From 8dc150b84edb4c894b53d0f74a18bd281f44fe90 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Fri, 24 Mar 2017 13:47:40 +0100 Subject: [PATCH 23/24] Explain determinant of empty matrix to Cholesky/PartialPivLu --- src/matrix/decomposition/cholesky.rs | 3 +++ src/matrix/decomposition/lu.rs | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/matrix/decomposition/cholesky.rs b/src/matrix/decomposition/cholesky.rs index 5dc5218..fa67686 100644 --- a/src/matrix/decomposition/cholesky.rs +++ b/src/matrix/decomposition/cholesky.rs @@ -169,6 +169,9 @@ impl<T> Cholesky<T> where T: 'static + 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 { let l_det = self.l.diag() .cloned() diff --git a/src/matrix/decomposition/lu.rs b/src/matrix/decomposition/lu.rs index 5178103..7913935 100644 --- a/src/matrix/decomposition/lu.rs +++ b/src/matrix/decomposition/lu.rs @@ -285,6 +285,9 @@ impl<T> PartialPivLu<T> 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, From 9ec0d1c960f51106fb041a34b5a85b793444c508 Mon Sep 17 00:00:00 2001 From: Andreas Longva <andreas.b.longva@gmail.com> Date: Fri, 24 Mar 2017 14:09:46 +0100 Subject: [PATCH 24/24] Reorder Cholesky and FullPivLu in decomposition docs --- src/matrix/decomposition/mod.rs | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index 6231ed3..54aea23 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -80,26 +80,26 @@ //! </tr> //! //! <tr> -//! <td><a href="struct.Cholesky.html">Cholesky</a></td> -//! <td>Square, symmetric positive definite</td> +//! <td><a href="struct.FullPivLu.html">FullPivLu</a></td> +//! <td>Square matrices</td> //! <td> //! <ul> //! <li>Linear system solving</li> //! <li>Matrix inverse</li> //! <li>Determinant computation</li> +//! <li>Rank computation</li> //! </ul> //! </td> //! </tr> //! //! <tr> -//! <td><a href="struct.FullPivLu.html">FullPivLu</a></td> -//! <td>Square matrices</td> +//! <td><a href="struct.Cholesky.html">Cholesky</a></td> +//! <td>Square, symmetric positive definite</td> //! <td> //! <ul> //! <li>Linear system solving</li> //! <li>Matrix inverse</li> //! <li>Determinant computation</li> -//! <li>Rank computation</li> //! </ul> //! </td> //! </tr>