1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

use num;
use num::Float;

use ApproxEq;
use Matrix;
use internalutil::{alloc_dirty_vec};

/// Cholesky Decomposition (for a real valued matrix).
///
/// Initial implementation based on JAMA.
///
/// For a symmetric, positive definite matrix A, the Cholesky decomposition
/// is an lower triangular matrix L so that A = L*L'.
///
/// If the matrix is not symmetric or positive definite, None is returned.
//
// Solve L one row at a time from row 1 to row n:
//   A = L * L'
//
//   a11 a12 .. .. a1n   L11   0  .. ..   0   L11 L21 .. .. Ln1
//   a21 a22 .. .. a2n   L21 L22  .. ..   0     0 L22 .. .. Ln2
//   ..  ..  .. .. ..  =  ..  ..  .. ..   0 *  ..  .. .. ..  ..
//   ..  ..  .. .. ..     ..  ..  .. ..   0    ..  .. .. ..  ..
//   an1 an2 .. .. ann   Ln1 Ln2  .. .. Lnn     0   0  0  0 Lnn
//
// Solve L[1][_]:
//   A[1][1] = L[1][1] * L[1][1]
//   => L[1][1] = sqrt(A[1][1]).
//   => L[1][2 .. n] = 0.
//
// Solve L[2][_]:
//   A[2][1] = L[2][1] * L[1][1]
//   => L[2][1] = A[2][1] / L[1][1].
//   A[2][2] = L[2][1] * L[2][1] + L[2][2] * L[2][2]
//   => L[2][2] = sqrt(A[2][2] - L[2][1] * L[2][1]).
//   => L[2][3 .. n] = 0.
//
// And in general: j = { 1 .. n }, i = { 1 .. n}
//   For j < i:
//     L[j][i] = 0.
//   For j = i:
//     L[j][i] = sqrt(A[j][i] - SUM { L[j][0 .. (j - 1)] * L'[(0 .. (j - 1)][j] })
//             = sqrt(A[j][j] - SUM { L[j][0 .. (j - 1)]^2 }
//   For j > i:
//     L[j][i] = (A[j][i] - SUM { L[i][0 .. (i - 1)] * L'[0 .. (i - 1)][j]}) / L[i][i].
//             = (A[j][i] - SUM { L[i][0 .. (i - 1)] * L[j][0 .. (i - 1)]) / L[i][i].
//
// As long as we follow the up->down, left->right order to compute the values, all the elements of L accessed on the right
// side will have been computed by the time they are needed.
//
pub struct CholeskyDecomposition<T> {
  l : Matrix<T>
}

impl<T : Float + ApproxEq<T>> CholeskyDecomposition<T> {
  pub fn new(m : &Matrix<T>) -> Option<CholeskyDecomposition<T>> {
    if m.rows() != m.cols() {
      return None
    }

    let n = m.rows();
    let mut data : Vec<T> = alloc_dirty_vec(n * n);

    for j in 0..n {
      // Solve row L[j].

      // d = SUM { L[j][0 .. (j - 1)]^2 }
      let mut d : T = num::zero();

      // Solve L[j][0 .. (j - 1)].
      for k in 0..j {
        // Solve L[j][k].

        // s = SUM { L[k][0 .. (k - 1)] * L'[0 .. (k - 1)][j] }
        //   = SUM { L[k][0 .. (k - 1)] * L[j][0 .. (k - 1) }
        let mut s : T = num::zero();
        for i in 0..k {
          s = s + data[k * n + i] * data[j * n + i];
        }

        // L[j][k] = (A[j][k] - SUM { L[k][0 .. (k - 1)] * L'[0 .. (k - 1)][j] }) / L[k][k].
        s = (m.get(j, k) - s) / data[k * n + k];
        data[j * n + k] = s;

        // Gather a sum of squres of L[j][0 .. (j - 1)] to d. Note: s = L[j][k].
        d = d + s * s;

        // Make sure input matrix is symmetric; Cholesky decomposition is not defined for non-symmetric matrixes.
        if m.get(k, j) != m.get(j, k) {
          return None
        }
      }

      // Solve L[j][j]. (Diagonals).
      // L[j][j] = sqrt(A[j][j] - SUM { L[j][0 .. (j - 1)]^2 }).
      d = m.get(j, j) - d;
      if d <= num::zero() {
        // A is not positive definite; Cholesky decomposition does not exists.
        return None
      }
      data[j * n + j] = d.sqrt();

      // Solve L[j][(j + 1) .. (n - 1)]. (Always zero as L is lower triangular).
      for k in (j + 1)..n {
        data[j * n + k] = num::zero();
      }
    }

    Some(CholeskyDecomposition { l : Matrix::new(n, n, data) })
  }

  #[inline]
  pub fn get_l<'lt>(&'lt self) -> &'lt Matrix<T> { &self.l }

  // Solve: Ax = b
  pub fn solve(&self, b : &Matrix<T>) -> Matrix<T> {
    let l = &self.l;
    assert!(l.rows() == b.rows());
    let n = l.rows();
    let mut xdata = b.get_data().clone();
    let nx = b.cols();

    // Solve L*Y = B
    for k in 0..n {
      for j in 0..nx {
        for i in 0..k {
          xdata[k * nx + j] = xdata[k * nx + j] - xdata[i * nx + j] * l.get_data()[k * n + i];
        }
        xdata[k * nx + j] = xdata[k * nx + j] / l.get_data()[k * n + k];
      }
    }

    // Solve L'*X = Y
    for k in (0..n).rev() {
      for j in 0..nx {
        for i in (k + 1)..n {
          xdata[k * nx + j] = xdata[k * nx + j] - xdata[i * nx + j] * l.get_data()[i * n + k];
        }
        xdata[k * nx + j] = xdata[k * nx + j] / l.get_data()[k * n + k];
      }
    }

    Matrix::new(n, nx, xdata)
  }
}

#[test]
fn cholesky_square_pos_def_test() {
  let a = m!(4.0, 12.0, -16.0; 12.0, 37.0, -43.0; -16.0, -43.0, 98.0);
  let c = CholeskyDecomposition::new(&a).unwrap();
  assert!(c.get_l() * c.get_l().t() == a);
  assert!(*c.get_l().get_data() == vec![2.0, 0.0, 0.0, 6.0, 1.0, 0.0, -8.0, 5.0, 3.0]);
}

#[test]
fn cholesky_not_pos_def_test() {
  let a = m!(4.0, 12.0, -16.0; 12.0, 37.0, 43.0; -16.0, 43.0, 98.0);
  assert!(CholeskyDecomposition::new(&a).is_none());
}

#[test]
fn cholesky_not_square_test() {
  let a = m!(4.0, 12.0, -16.0; 12.0, 37.0, 43.0);
  assert!(CholeskyDecomposition::new(&a).is_none());
}

#[test]
fn cholesky_solve_test() {
  let a = m!(2.0, 1.0, 0.0; 1.0, 1.0, 0.0; 0.0, 0.0, 1.0);
  let c = CholeskyDecomposition::new(&a).unwrap();
  let b = m!(1.0; 2.0; 3.0);
  assert!(c.solve(&b).approx_eq(&m!(-1.0; 3.0; 3.0)));
}

#[test]
#[should_panic]
fn cholesky_solve_test_incompatible() {
  let a = m!(2.0, 1.0, 0.0; 1.0, 1.0, 0.0; 0.0, 0.0, 1.0);
  let c = CholeskyDecomposition::new(&a).unwrap();
  let b = m!(1.0; 2.0; 3.0; 4.0);
  let _ = c.solve(&b);
}