Merge pull request #147 from dshizzle/master

Implemented Cholesky decomposition

src/lib.rs

@@ -150,7 +150,8 @@ pub use structs::{
 
 pub use linalg::{
     qr,
-    householder_matrix
+    householder_matrix,
+    cholesky
 };
 
 mod structs;

src/linalg/decompositions.rs

@@ -115,3 +115,51 @@ pub fn eigen_qr<N, V, VS, M>(m: &M, eps: &N, niter: usize) -> (M, V)
 
     (eigenvectors, eigenvalues.diag())
 }
+
+/// Cholesky decomposition G of a square symmetric positive definite matrix A, such that A = G * G^T
+///
+/// # Arguments
+/// * `m` - square symmetric positive definite matrix to decompose
+pub fn cholesky<N, V, VS, M>(m: &M) -> Result<M, &'static str>
+    where N:  BaseFloat,
+          VS: Indexable<usize, N> + Norm<N>,
+          M:  Indexable<(usize, usize), N> + SquareMat<N, V> + Add<M, Output = M> +
+              Sub<M, Output = M> + ColSlice<VS> +
+              ApproxEq<N> + Copy {
+
+    let mut out = m.clone().transpose();
+
+    if !ApproxEq::approx_eq(&out, &m) {
+        return Err("Cholesky: Input matrix is not symmetric");
+    }
+
+    for i in 0..out.nrows() {
+        for j in 0..(i+1) {
+
+            let mut sum: N = out[(i,j)];
+
+            for k in 0..j {
+                sum = sum - out[(i, k)] * out[(j, k)];
+            }
+
+            if i > j {
+                out[(i, j)] = sum / out[(j, j)];
+            }
+            else if sum > N::zero() {
+                out[(i,i)] = sum.sqrt();
+            }
+            else {
+                return Err("Cholesky: Input matrix is not positive definite to machine precision");
+            }
+        }
+    }
+
+    for i in 0..out.nrows() {
+        for j in i+1..out.ncols() {
+            out[(i,j)] = N::zero();
+        }
+
+    }
+
+    return Ok(out);
+}

src/linalg/mod.rs

@@ -1,4 +1,4 @@
 
-pub use self::decompositions::{qr, eigen_qr, householder_matrix};
+pub use self::decompositions::{qr, eigen_qr, householder_matrix, cholesky};
 
 mod decompositions;

tests/mat.rs

@@ -3,7 +3,7 @@ extern crate rand;
 
 use rand::random;
 use na::{Vec1, Vec3, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6, Rot2, Rot3, Persp3, PerspMat3, Ortho3,
-         OrthoMat3, DMat, DVec, Row, Col, BaseFloat};
+         OrthoMat3, DMat, DVec, Row, Col, BaseFloat, Diag};
 
 macro_rules! test_inv_mat_impl(
   ($t: ty) => (
@@ -41,6 +41,28 @@ macro_rules! test_qr_impl(
   );
 );
 
+macro_rules! test_cholesky_impl(
+  ($t: ty) => (
+    for _ in (0usize .. 10000) {
+      
+      // construct symmetric positive definite matrix
+      let mut randmat : $t = random();
+      let mut diagmat : $t = Diag::from_diag(&na::diag(&randmat));
+
+      diagmat = na::abs(&diagmat) + 1.0;
+      randmat = randmat * diagmat * na::transpose(&randmat);
+
+      let result = na::cholesky(&randmat);
+
+      assert!(result.is_ok());
+
+      let v = result.unwrap();
+      let recomp = v * na::transpose(&v);
+      assert!(na::approx_eq(&randmat,  &recomp));
+    }
+  );
+);
+
 // NOTE: deactivated untile we get a better convergence rate.
 // macro_rules! test_eigen_qr_impl(
 //     ($t: ty) => {
@@ -600,3 +622,70 @@ fn test_ortho() {
     assert!(na::approx_eq(&pm.znear(),  &24.0));
     assert!(na::approx_eq(&pm.zfar(),   &61.0));
 }
+
+#[test]
+fn test_cholesky_const() {
+    
+    let a : Mat3<f64> = Mat3::<f64>::new(1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 2.0, 3.0);
+    let g : Mat3<f64> = Mat3::<f64>::new(1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0);
+
+    let result = na::cholesky(&a);
+
+    assert!(result.is_ok());
+
+    let v = result.unwrap();
+    assert!(na::approx_eq(&v, &g));
+
+    let recomp = v * na::transpose(&v);
+    assert!(na::approx_eq(&recomp, &a));
+}
+
+#[test]
+fn test_cholesky_not_spd() {
+    
+    let a : Mat3<f64> = Mat3::<f64>::new(1.0, 2.0, 3.0, 3.0, 2.0, 1.0, 1.0, 1.0, 1.0);
+
+    let result = na::cholesky(&a);
+
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_cholesky_not_symmetric() {
+    
+    let a : Mat2<f64> = Mat2::<f64>::new(1.0, 1.0, -1.0, 1.0);
+
+    let result = na::cholesky(&a);
+
+    assert!(result.is_err());
+}
+
+#[test]
+fn test_cholesky_mat1() {
+    test_cholesky_impl!(Mat1<f64>);
+}
+
+#[test]
+fn test_cholesky_mat2() {
+    test_cholesky_impl!(Mat2<f64>);
+}
+
+#[test]
+fn test_cholesky_mat3() {
+    test_cholesky_impl!(Mat3<f64>);
+}
+
+#[test]
+fn test_cholesky_mat4() {
+    test_cholesky_impl!(Mat4<f64>);
+}
+
+#[test]
+fn test_cholesky_mat5() {
+    test_cholesky_impl!(Mat5<f64>);
+}
+
+#[test]
+fn test_cholesky_mat6() {
+    test_cholesky_impl!(Mat6<f64>);
+}