CholeskyDecomposition.php
3.3 KB
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
<?php
/**
* @package JAMA
*
* Cholesky decomposition class
*
* 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, the constructor
* returns a partial decomposition and sets an internal flag that may
* be queried by the isSPD() method.
*
* @author Paul Meagher
* @author Michael Bommarito
* @version 1.2
*/
class CholeskyDecomposition {
/**
* Decomposition storage
* @var array
* @access private
*/
private $L = array();
/**
* Matrix row and column dimension
* @var int
* @access private
*/
private $m;
/**
* Symmetric positive definite flag
* @var boolean
* @access private
*/
private $isspd = true;
/**
* CholeskyDecomposition
*
* Class constructor - decomposes symmetric positive definite matrix
* @param mixed Matrix square symmetric positive definite matrix
*/
public function __construct($A = null) {
if ($A instanceof Matrix) {
$this->L = $A->getArray();
$this->m = $A->getRowDimension();
for($i = 0; $i < $this->m; ++$i) {
for($j = $i; $j < $this->m; ++$j) {
for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
$sum -= $this->L[$i][$k] * $this->L[$j][$k];
}
if ($i == $j) {
if ($sum >= 0) {
$this->L[$i][$i] = sqrt($sum);
} else {
$this->isspd = false;
}
} else {
if ($this->L[$i][$i] != 0) {
$this->L[$j][$i] = $sum / $this->L[$i][$i];
}
}
}
for ($k = $i+1; $k < $this->m; ++$k) {
$this->L[$i][$k] = 0.0;
}
}
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
}
} // function __construct()
/**
* Is the matrix symmetric and positive definite?
*
* @return boolean
*/
public function isSPD() {
return $this->isspd;
} // function isSPD()
/**
* getL
*
* Return triangular factor.
* @return Matrix Lower triangular matrix
*/
public function getL() {
return new Matrix($this->L);
} // function getL()
/**
* Solve A*X = B
*
* @param $B Row-equal matrix
* @return Matrix L * L' * X = B
*/
public function solve($B = null) {
if ($B instanceof Matrix) {
if ($B->getRowDimension() == $this->m) {
if ($this->isspd) {
$X = $B->getArrayCopy();
$nx = $B->getColumnDimension();
for ($k = 0; $k < $this->m; ++$k) {
for ($i = $k + 1; $i < $this->m; ++$i) {
for ($j = 0; $j < $nx; ++$j) {
$X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
}
}
for ($j = 0; $j < $nx; ++$j) {
$X[$k][$j] /= $this->L[$k][$k];
}
}
for ($k = $this->m - 1; $k >= 0; --$k) {
for ($j = 0; $j < $nx; ++$j) {
$X[$k][$j] /= $this->L[$k][$k];
}
for ($i = 0; $i < $k; ++$i) {
for ($j = 0; $j < $nx; ++$j) {
$X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
}
}
}
return new Matrix($X, $this->m, $nx);
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));
}
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionException));
}
} else {
throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException));
}
} // function solve()
} // class CholeskyDecomposition