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