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
use matrix::{Matrix, BaseMatrix};
use error::{Error, ErrorKind};
use std::any::Any;
use libnum::Float;
impl<T> Matrix<T>
where T: Any + Float
{
/// Cholesky decomposition
///
/// Returns the cholesky decomposition of a positive definite matrix.
///
/// # Examples
///
/// ```
/// # #[macro_use] extern crate rulinalg; fn main() {
/// use rulinalg::matrix::Matrix;
///
/// let m = matrix![1.0, 0.5, 0.5;
/// 0.5, 1.0, 0.5;
/// 0.5, 0.5, 1.0];
///
/// let l = m.cholesky();
/// # }
/// ```
///
/// # Panics
///
/// - The matrix is not square.
///
/// # Failures
///
/// - Matrix is not positive definite.
pub fn cholesky(&self) -> Result<Matrix<T>, Error> {
assert!(self.rows == self.cols,
"Matrix must be square for Cholesky decomposition.");
let mut new_data = Vec::<T>::with_capacity(self.rows() * self.cols());
for i in 0..self.rows() {
for j in 0..self.cols() {
if j > i {
new_data.push(T::zero());
continue;
}
let mut sum = T::zero();
for k in 0..j {
sum = sum + (new_data[i * self.cols() + k] * new_data[j * self.cols() + k]);
}
if j == i {
new_data.push((self[[i, i]] - sum).sqrt());
} else {
let p = (self[[i, j]] - sum) / new_data[j * self.cols + j];
if !p.is_finite() {
return Err(Error::new(ErrorKind::DecompFailure,
"Matrix is not positive definite."));
} else {
}
new_data.push(p);
}
}
}
Ok(Matrix {
rows: self.rows(),
cols: self.cols(),
data: new_data,
})
}
}
#[cfg(test)]
mod tests {
use matrix::Matrix;
#[test]
#[should_panic]
fn test_non_square_cholesky() {
let a = Matrix::<f64>::ones(2, 3);
let _ = a.cholesky();
}
}