18
18
* @author Michael Bommarito
19
19
* @version 1.2
20
20
*/
21
- class CholeskyDecomposition {
21
+ class CholeskyDecomposition
22
+ {
22
23
23
24
/**
24
25
* Decomposition storage
25
26
* @var array
26
27
* @access private
27
28
*/
28
- private $ L = array ();
29
-
30
- /**
31
- * Matrix row and column dimension
32
- * @var int
33
- * @access private
34
- */
35
- private $ m ;
36
-
37
- /**
38
- * Symmetric positive definite flag
39
- * @var boolean
40
- * @access private
41
- */
42
- private $ isspd = true ;
43
-
44
-
45
- /**
46
- * CholeskyDecomposition
47
- *
48
- * Class constructor - decomposes symmetric positive definite matrix
49
- * @param mixed Matrix square symmetric positive definite matrix
50
- */
51
- public function __construct ($ A = null ) {
52
- if ($ A instanceof Matrix) {
53
- $ this ->L = $ A ->getArray ();
54
- $ this ->m = $ A ->getRowDimension ();
55
-
56
- for ($ i = 0 ; $ i < $ this ->m ; ++$ i ) {
57
- for ($ j = $ i ; $ j < $ this ->m ; ++$ j ) {
58
- for ($ sum = $ this ->L [$ i ][$ j ], $ k = $ i - 1 ; $ k >= 0 ; --$ k ) {
59
- $ sum -= $ this ->L [$ i ][$ k ] * $ this ->L [$ j ][$ k ];
60
- }
61
- if ($ i == $ j ) {
62
- if ($ sum >= 0 ) {
63
- $ this ->L [$ i ][$ i ] = sqrt ($ sum );
64
- } else {
65
- $ this ->isspd = false ;
66
- }
67
- } else {
68
- if ($ this ->L [$ i ][$ i ] != 0 ) {
69
- $ this ->L [$ j ][$ i ] = $ sum / $ this ->L [$ i ][$ i ];
70
- }
71
- }
72
- }
73
-
74
- for ($ k = $ i +1 ; $ k < $ this ->m ; ++$ k ) {
75
- $ this ->L [$ i ][$ k ] = 0.0 ;
76
- }
77
- }
78
- } else {
79
- throw new Exception (JAMAError (ArgumentTypeException));
80
- }
81
- } // function __construct()
82
-
83
-
84
- /**
85
- * Is the matrix symmetric and positive definite?
86
- *
87
- * @return boolean
88
- */
89
- public function isSPD () {
90
- return $ this ->isspd ;
91
- } // function isSPD()
92
-
93
-
94
- /**
95
- * getL
96
- *
97
- * Return triangular factor.
98
- * @return Matrix Lower triangular matrix
99
- */
100
- public function getL () {
101
- return new Matrix ($ this ->L );
102
- } // function getL()
103
-
104
-
105
- /**
106
- * Solve A*X = B
107
- *
108
- * @param $B Row-equal matrix
109
- * @return Matrix L * L' * X = B
110
- */
111
- public function solve ($ B = null ) {
112
- if ($ B instanceof Matrix) {
113
- if ($ B ->getRowDimension () == $ this ->m ) {
114
- if ($ this ->isspd ) {
115
- $ X = $ B ->getArrayCopy ();
116
- $ nx = $ B ->getColumnDimension ();
117
-
118
- for ($ k = 0 ; $ k < $ this ->m ; ++$ k ) {
119
- for ($ i = $ k + 1 ; $ i < $ this ->m ; ++$ i ) {
120
- for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
121
- $ X [$ i ][$ j ] -= $ X [$ k ][$ j ] * $ this ->L [$ i ][$ k ];
122
- }
123
- }
124
- for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
125
- $ X [$ k ][$ j ] /= $ this ->L [$ k ][$ k ];
126
- }
127
- }
128
-
129
- for ($ k = $ this ->m - 1 ; $ k >= 0 ; --$ k ) {
130
- for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
131
- $ X [$ k ][$ j ] /= $ this ->L [$ k ][$ k ];
132
- }
133
- for ($ i = 0 ; $ i < $ k ; ++$ i ) {
134
- for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
135
- $ X [$ i ][$ j ] -= $ X [$ k ][$ j ] * $ this ->L [$ k ][$ i ];
136
- }
137
- }
138
- }
139
-
140
- return new Matrix ($ X , $ this ->m , $ nx );
141
- } else {
142
- throw new Exception (JAMAError (MatrixSPDException));
143
- }
144
- } else {
145
- throw new Exception (JAMAError (MatrixDimensionException));
146
- }
147
- } else {
148
- throw new Exception (JAMAError (ArgumentTypeException));
149
- }
150
- } // function solve()
151
-
152
- } // class CholeskyDecomposition
29
+ private $ L = [];
30
+
31
+ /**
32
+ * Matrix row and column dimension
33
+ * @var int
34
+ * @access private
35
+ */
36
+ private $ m ;
37
+
38
+ /**
39
+ * Symmetric positive definite flag
40
+ * @var boolean
41
+ * @access private
42
+ */
43
+ private $ isspd = true ;
44
+
45
+
46
+ /**
47
+ * CholeskyDecomposition
48
+ *
49
+ * Class constructor - decomposes symmetric positive definite matrix
50
+ * @param mixed Matrix square symmetric positive definite matrix
51
+ */
52
+ public function __construct ($ A = null )
53
+ {
54
+ if ($ A instanceof Matrix) {
55
+ $ this ->L = $ A ->getArray ();
56
+ $ this ->m = $ A ->getRowDimension ();
57
+
58
+ for ($ i = 0 ; $ i < $ this ->m ; ++$ i ) {
59
+ for ($ j = $ i ; $ j < $ this ->m ; ++$ j ) {
60
+ for ($ sum = $ this ->L [$ i ][$ j ], $ k = $ i - 1 ; $ k >= 0 ; --$ k ) {
61
+ $ sum -= $ this ->L [$ i ][$ k ] * $ this ->L [$ j ][$ k ];
62
+ }
63
+ if ($ i == $ j ) {
64
+ if ($ sum >= 0 ) {
65
+ $ this ->L [$ i ][$ i ] = sqrt ($ sum );
66
+ } else {
67
+ $ this ->isspd = false ;
68
+ }
69
+ } else {
70
+ if ($ this ->L [$ i ][$ i ] != 0 ) {
71
+ $ this ->L [$ j ][$ i ] = $ sum / $ this ->L [$ i ][$ i ];
72
+ }
73
+ }
74
+ }
75
+
76
+ for ($ k = $ i +1 ; $ k < $ this ->m ; ++$ k ) {
77
+ $ this ->L [$ i ][$ k ] = 0.0 ;
78
+ }
79
+ }
80
+ } else {
81
+ throw new Exception (JAMAError (ArgumentTypeException));
82
+ }
83
+ } // function __construct()
84
+
85
+
86
+ /**
87
+ * Is the matrix symmetric and positive definite?
88
+ *
89
+ * @return boolean
90
+ */
91
+ public function isSPD ()
92
+ {
93
+ return $ this ->isspd ;
94
+ } // function isSPD()
95
+
96
+
97
+ /**
98
+ * getL
99
+ *
100
+ * Return triangular factor.
101
+ * @return Matrix Lower triangular matrix
102
+ */
103
+ public function getL ()
104
+ {
105
+ return new Matrix ($ this ->L );
106
+ } // function getL()
107
+
108
+
109
+ /**
110
+ * Solve A*X = B
111
+ *
112
+ * @param $B Row-equal matrix
113
+ * @return Matrix L * L' * X = B
114
+ */
115
+ public function solve ($ B = null )
116
+ {
117
+ if ($ B instanceof Matrix) {
118
+ if ($ B ->getRowDimension () == $ this ->m ) {
119
+ if ($ this ->isspd ) {
120
+ $ X = $ B ->getArrayCopy ();
121
+ $ nx = $ B ->getColumnDimension ();
122
+
123
+ for ($ k = 0 ; $ k < $ this ->m ; ++$ k ) {
124
+ for ($ i = $ k + 1 ; $ i < $ this ->m ; ++$ i ) {
125
+ for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
126
+ $ X [$ i ][$ j ] -= $ X [$ k ][$ j ] * $ this ->L [$ i ][$ k ];
127
+ }
128
+ }
129
+ for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
130
+ $ X [$ k ][$ j ] /= $ this ->L [$ k ][$ k ];
131
+ }
132
+ }
133
+
134
+ for ($ k = $ this ->m - 1 ; $ k >= 0 ; --$ k ) {
135
+ for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
136
+ $ X [$ k ][$ j ] /= $ this ->L [$ k ][$ k ];
137
+ }
138
+ for ($ i = 0 ; $ i < $ k ; ++$ i ) {
139
+ for ($ j = 0 ; $ j < $ nx ; ++$ j ) {
140
+ $ X [$ i ][$ j ] -= $ X [$ k ][$ j ] * $ this ->L [$ k ][$ i ];
141
+ }
142
+ }
143
+ }
144
+
145
+ return new Matrix ($ X , $ this ->m , $ nx );
146
+ } else {
147
+ throw new Exception (JAMAError (MatrixSPDException));
148
+ }
149
+ } else {
150
+ throw new Exception (JAMAError (MatrixDimensionException));
151
+ }
152
+ } else {
153
+ throw new Exception (JAMAError (ArgumentTypeException));
154
+ }
155
+ } // function solve()
156
+ } // class CholeskyDecomposition
0 commit comments