1: <?php
2: /**
3: * @package JAMA
4: *
5: * For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
6: * orthogonal matrix Q and an n-by-n upper triangular matrix R so that
7: * A = Q*R.
8: *
9: * The QR decompostion always exists, even if the matrix does not have
10: * full rank, so the constructor will never fail. The primary use of the
11: * QR decomposition is in the least squares solution of nonsquare systems
12: * of simultaneous linear equations. This will fail if isFullRank()
13: * returns false.
14: *
15: * @author Paul Meagher
16: * @license PHP v3.0
17: * @version 1.1
18: */
19: class PHPExcel_Shared_JAMA_QRDecomposition {
20:
21: const MatrixRankException = "Can only perform operation on full-rank matrix.";
22:
23: /**
24: * Array for internal storage of decomposition.
25: * @var array
26: */
27: private $QR = array();
28:
29: /**
30: * Row dimension.
31: * @var integer
32: */
33: private $m;
34:
35: /**
36: * Column dimension.
37: * @var integer
38: */
39: private $n;
40:
41: /**
42: * Array for internal storage of diagonal of R.
43: * @var array
44: */
45: private $Rdiag = array();
46:
47:
48: /**
49: * QR Decomposition computed by Householder reflections.
50: *
51: * @param matrix $A Rectangular matrix
52: * @return Structure to access R and the Householder vectors and compute Q.
53: */
54: public function __construct($A) {
55: if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
56: // Initialize.
57: $this->QR = $A->getArrayCopy();
58: $this->m = $A->getRowDimension();
59: $this->n = $A->getColumnDimension();
60: // Main loop.
61: for ($k = 0; $k < $this->n; ++$k) {
62: // Compute 2-norm of k-th column without under/overflow.
63: $nrm = 0.0;
64: for ($i = $k; $i < $this->m; ++$i) {
65: $nrm = hypo($nrm, $this->QR[$i][$k]);
66: }
67: if ($nrm != 0.0) {
68: // Form k-th Householder vector.
69: if ($this->QR[$k][$k] < 0) {
70: $nrm = -$nrm;
71: }
72: for ($i = $k; $i < $this->m; ++$i) {
73: $this->QR[$i][$k] /= $nrm;
74: }
75: $this->QR[$k][$k] += 1.0;
76: // Apply transformation to remaining columns.
77: for ($j = $k+1; $j < $this->n; ++$j) {
78: $s = 0.0;
79: for ($i = $k; $i < $this->m; ++$i) {
80: $s += $this->QR[$i][$k] * $this->QR[$i][$j];
81: }
82: $s = -$s/$this->QR[$k][$k];
83: for ($i = $k; $i < $this->m; ++$i) {
84: $this->QR[$i][$j] += $s * $this->QR[$i][$k];
85: }
86: }
87: }
88: $this->Rdiag[$k] = -$nrm;
89: }
90: } else {
91: throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
92: }
93: } // function __construct()
94:
95:
96: /**
97: * Is the matrix full rank?
98: *
99: * @return boolean true if R, and hence A, has full rank, else false.
100: */
101: public function isFullRank() {
102: for ($j = 0; $j < $this->n; ++$j) {
103: if ($this->Rdiag[$j] == 0) {
104: return false;
105: }
106: }
107: return true;
108: } // function isFullRank()
109:
110:
111: /**
112: * Return the Householder vectors
113: *
114: * @return Matrix Lower trapezoidal matrix whose columns define the reflections
115: */
116: public function getH() {
117: for ($i = 0; $i < $this->m; ++$i) {
118: for ($j = 0; $j < $this->n; ++$j) {
119: if ($i >= $j) {
120: $H[$i][$j] = $this->QR[$i][$j];
121: } else {
122: $H[$i][$j] = 0.0;
123: }
124: }
125: }
126: return new PHPExcel_Shared_JAMA_Matrix($H);
127: } // function getH()
128:
129:
130: /**
131: * Return the upper triangular factor
132: *
133: * @return Matrix upper triangular factor
134: */
135: public function getR() {
136: for ($i = 0; $i < $this->n; ++$i) {
137: for ($j = 0; $j < $this->n; ++$j) {
138: if ($i < $j) {
139: $R[$i][$j] = $this->QR[$i][$j];
140: } elseif ($i == $j) {
141: $R[$i][$j] = $this->Rdiag[$i];
142: } else {
143: $R[$i][$j] = 0.0;
144: }
145: }
146: }
147: return new PHPExcel_Shared_JAMA_Matrix($R);
148: } // function getR()
149:
150:
151: /**
152: * Generate and return the (economy-sized) orthogonal factor
153: *
154: * @return Matrix orthogonal factor
155: */
156: public function getQ() {
157: for ($k = $this->n-1; $k >= 0; --$k) {
158: for ($i = 0; $i < $this->m; ++$i) {
159: $Q[$i][$k] = 0.0;
160: }
161: $Q[$k][$k] = 1.0;
162: for ($j = $k; $j < $this->n; ++$j) {
163: if ($this->QR[$k][$k] != 0) {
164: $s = 0.0;
165: for ($i = $k; $i < $this->m; ++$i) {
166: $s += $this->QR[$i][$k] * $Q[$i][$j];
167: }
168: $s = -$s/$this->QR[$k][$k];
169: for ($i = $k; $i < $this->m; ++$i) {
170: $Q[$i][$j] += $s * $this->QR[$i][$k];
171: }
172: }
173: }
174: }
175: /*
176: for($i = 0; $i < count($Q); ++$i) {
177: for($j = 0; $j < count($Q); ++$j) {
178: if(! isset($Q[$i][$j]) ) {
179: $Q[$i][$j] = 0;
180: }
181: }
182: }
183: */
184: return new PHPExcel_Shared_JAMA_Matrix($Q);
185: } // function getQ()
186:
187:
188: /**
189: * Least squares solution of A*X = B
190: *
191: * @param Matrix $B A Matrix with as many rows as A and any number of columns.
192: * @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
193: */
194: public function solve($B) {
195: if ($B->getRowDimension() == $this->m) {
196: if ($this->isFullRank()) {
197: // Copy right hand side
198: $nx = $B->getColumnDimension();
199: $X = $B->getArrayCopy();
200: // Compute Y = transpose(Q)*B
201: for ($k = 0; $k < $this->n; ++$k) {
202: for ($j = 0; $j < $nx; ++$j) {
203: $s = 0.0;
204: for ($i = $k; $i < $this->m; ++$i) {
205: $s += $this->QR[$i][$k] * $X[$i][$j];
206: }
207: $s = -$s/$this->QR[$k][$k];
208: for ($i = $k; $i < $this->m; ++$i) {
209: $X[$i][$j] += $s * $this->QR[$i][$k];
210: }
211: }
212: }
213: // Solve R*X = Y;
214: for ($k = $this->n-1; $k >= 0; --$k) {
215: for ($j = 0; $j < $nx; ++$j) {
216: $X[$k][$j] /= $this->Rdiag[$k];
217: }
218: for ($i = 0; $i < $k; ++$i) {
219: for ($j = 0; $j < $nx; ++$j) {
220: $X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
221: }
222: }
223: }
224: $X = new PHPExcel_Shared_JAMA_Matrix($X);
225: return ($X->getMatrix(0, $this->n-1, 0, $nx));
226: } else {
227: throw new PHPExcel_Calculation_Exception(self::MatrixRankException);
228: }
229: } else {
230: throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
231: }
232: } // function solve()
233:
234: } // PHPExcel_Shared_JAMA_class QRDecomposition
235: