1
1
//! # Cholesky Decomposition
2
2
//!
3
- //! every positive definite matrix \\(A \in R^{n \times n}\\) can be factored as
3
+ //! every positive definite matrix \\(A \in R^{n \times n}\\) can be factored as
4
4
//!
5
5
//! \\[A = R^TR\\]
6
- //!
6
+ //!
7
7
//! where \\(R\\) is upper triangular matrix with positive diagonal elements
8
8
//!
9
9
//! Example:
12
12
//! use crate::smartcore::linalg::cholesky::*;
13
13
//!
14
14
//! let A = DenseMatrix::from_2d_array(&[
15
- //! &[25., 15., -5.],
16
- //! &[15., 18., 0.],
15
+ //! &[25., 15., -5.],
16
+ //! &[15., 18., 0.],
17
17
//! &[-5., 0., 11.]
18
18
//! ]);
19
19
//!
@@ -41,14 +41,14 @@ use crate::math::num::RealNumber;
41
41
/// Results of Cholesky decomposition.
42
42
pub struct Cholesky < T : RealNumber , M : BaseMatrix < T > > {
43
43
R : M ,
44
- t : PhantomData < T >
44
+ t : PhantomData < T > ,
45
45
}
46
46
47
47
impl < T : RealNumber , M : BaseMatrix < T > > Cholesky < T , M > {
48
48
pub ( crate ) fn new ( R : M ) -> Cholesky < T , M > {
49
49
Cholesky {
50
50
R : R ,
51
- t : PhantomData
51
+ t : PhantomData ,
52
52
}
53
53
}
54
54
@@ -65,10 +65,10 @@ impl<T: RealNumber, M: BaseMatrix<T>> Cholesky<T, M> {
65
65
}
66
66
}
67
67
R
68
- }
69
-
68
+ }
69
+
70
70
/// Get upper triangular matrix.
71
- pub fn U ( & self ) -> M {
71
+ pub fn U ( & self ) -> M {
72
72
let ( n, _) = self . R . shape ( ) ;
73
73
let mut R = M :: zeros ( n, n) ;
74
74
@@ -80,20 +80,20 @@ impl<T: RealNumber, M: BaseMatrix<T>> Cholesky<T, M> {
80
80
}
81
81
}
82
82
R
83
- }
83
+ }
84
84
85
85
/// Solves Ax = b
86
- pub ( crate ) fn solve ( & self , mut b : M ) -> Result < M , Failed > {
87
-
86
+ pub ( crate ) fn solve ( & self , mut b : M ) -> Result < M , Failed > {
88
87
let ( bn, m) = b. shape ( ) ;
89
88
let ( rn, _) = self . R . shape ( ) ;
90
89
91
90
if bn != rn {
92
- return Err ( Failed :: because ( FailedError :: SolutionFailed , & format ! (
93
- "Can't solve Ax = b for x. Number of rows in b != number of rows in R."
94
- ) ) ) ;
91
+ return Err ( Failed :: because (
92
+ FailedError :: SolutionFailed ,
93
+ & format ! ( "Can't solve Ax = b for x. Number of rows in b != number of rows in R." ) ,
94
+ ) ) ;
95
95
}
96
-
96
+
97
97
for k in 0 ..bn {
98
98
for j in 0 ..m {
99
99
for i in 0 ..k {
@@ -102,7 +102,7 @@ impl<T: RealNumber, M: BaseMatrix<T>> Cholesky<T, M> {
102
102
b. div_element_mut ( k, j, self . R . get ( k, k) ) ;
103
103
}
104
104
}
105
-
105
+
106
106
for k in ( 0 ..bn) . rev ( ) {
107
107
for j in 0 ..m {
108
108
for i in k + 1 ..bn {
@@ -128,11 +128,12 @@ pub trait CholeskyDecomposableMatrix<T: RealNumber>: BaseMatrix<T> {
128
128
let ( m, n) = self . shape ( ) ;
129
129
130
130
if m != n {
131
- return Err ( Failed :: because ( FailedError :: DecompositionFailed , & format ! (
132
- "Can't do Cholesky decomposition on a non-square matrix"
133
- ) ) ) ;
131
+ return Err ( Failed :: because (
132
+ FailedError :: DecompositionFailed ,
133
+ & format ! ( "Can't do Cholesky decomposition on a non-square matrix" ) ,
134
+ ) ) ;
134
135
}
135
-
136
+
136
137
for j in 0 ..n {
137
138
let mut d = T :: zero ( ) ;
138
139
for k in 0 ..j {
@@ -147,9 +148,10 @@ pub trait CholeskyDecomposableMatrix<T: RealNumber>: BaseMatrix<T> {
147
148
d = self . get ( j, j) - d;
148
149
149
150
if d < T :: zero ( ) {
150
- return Err ( Failed :: because ( FailedError :: DecompositionFailed , & format ! (
151
- "The matrix is not positive definite."
152
- ) ) ) ;
151
+ return Err ( Failed :: because (
152
+ FailedError :: DecompositionFailed ,
153
+ & format ! ( "The matrix is not positive definite." ) ,
154
+ ) ) ;
153
155
}
154
156
155
157
self . set ( j, j, d. sqrt ( ) ) ;
@@ -172,35 +174,33 @@ mod tests {
172
174
#[ test]
173
175
fn cholesky_decompose ( ) {
174
176
let a = DenseMatrix :: from_2d_array ( & [ & [ 25. , 15. , -5. ] , & [ 15. , 18. , 0. ] , & [ -5. , 0. , 11. ] ] ) ;
175
- let l = DenseMatrix :: from_2d_array ( & [
176
- & [ 5.0 , 0.0 , 0.0 ] ,
177
- & [ 3.0 , 3.0 , 0.0 ] ,
178
- & [ -1.0 , 1.0 , 3.0 ] ,
179
- ] ) ;
180
- let u = DenseMatrix :: from_2d_array ( & [
181
- & [ 5.0 , 3.0 , -1.0 ] ,
182
- & [ 0.0 , 3.0 , 1.0 ] ,
183
- & [ 0.0 , 0.0 , 3.0 ] ,
184
- ] ) ;
177
+ let l =
178
+ DenseMatrix :: from_2d_array ( & [ & [ 5.0 , 0.0 , 0.0 ] , & [ 3.0 , 3.0 , 0.0 ] , & [ -1.0 , 1.0 , 3.0 ] ] ) ;
179
+ let u =
180
+ DenseMatrix :: from_2d_array ( & [ & [ 5.0 , 3.0 , -1.0 ] , & [ 0.0 , 3.0 , 1.0 ] , & [ 0.0 , 0.0 , 3.0 ] ] ) ;
185
181
let cholesky = a. cholesky ( ) . unwrap ( ) ;
186
-
187
- assert ! ( cholesky. L ( ) . abs( ) . approximate_eq( & l. abs( ) , 1e-4 ) ) ;
188
- assert ! ( cholesky. U ( ) . abs( ) . approximate_eq( & u. abs( ) , 1e-4 ) ) ;
189
- assert ! ( cholesky. L ( ) . matmul( & cholesky. U ( ) ) . abs( ) . approximate_eq( & a. abs( ) , 1e-4 ) ) ;
182
+
183
+ assert ! ( cholesky. L ( ) . abs( ) . approximate_eq( & l. abs( ) , 1e-4 ) ) ;
184
+ assert ! ( cholesky. U ( ) . abs( ) . approximate_eq( & u. abs( ) , 1e-4 ) ) ;
185
+ assert ! ( cholesky
186
+ . L ( )
187
+ . matmul( & cholesky. U ( ) )
188
+ . abs( )
189
+ . approximate_eq( & a. abs( ) , 1e-4 ) ) ;
190
190
}
191
191
192
192
#[ test]
193
193
fn cholesky_solve_mut ( ) {
194
194
let a = DenseMatrix :: from_2d_array ( & [ & [ 25. , 15. , -5. ] , & [ 15. , 18. , 0. ] , & [ -5. , 0. , 11. ] ] ) ;
195
195
let b = DenseMatrix :: from_2d_array ( & [ & [ 40. , 51. , 28. ] ] ) ;
196
- let expected = DenseMatrix :: from_2d_array ( & [
197
- & [ 1.0 , 2.0 , 3.0 ]
198
- ] ) ;
199
-
196
+ let expected = DenseMatrix :: from_2d_array ( & [ & [ 1.0 , 2.0 , 3.0 ] ] ) ;
197
+
200
198
let cholesky = a. cholesky ( ) . unwrap ( ) ;
201
199
202
- assert ! ( cholesky. solve( b. transpose( ) ) . unwrap( ) . transpose( ) . approximate_eq( & expected, 1e-4 ) ) ;
203
-
200
+ assert ! ( cholesky
201
+ . solve( b. transpose( ) )
202
+ . unwrap( )
203
+ . transpose( )
204
+ . approximate_eq( & expected, 1e-4 ) ) ;
204
205
}
205
-
206
206
}
0 commit comments