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>