ILIAS  release_5-3 Revision v5.3.23-19-g915713cf615
QRDecomposition.php
Go to the documentation of this file.
1<?php
20
21 const MatrixRankException = "Can only perform operation on full-rank matrix.";
22
27 private $QR = array();
28
33 private $m;
34
39 private $n;
40
45 private $Rdiag = array();
46
47
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 {
92 }
93 } // function __construct()
94
95
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
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
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
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
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 {
228 }
229 } else {
231 }
232 } // function solve()
233
234} // PHPExcel_Shared_JAMA_class QRDecomposition
hypo($a, $b)
Definition: Maths.php:14
An exception for terminatinating execution or to throw for unit testing.
PHPExcel root directory.
Definition: Matrix.php:27
getR()
Return the upper triangular factor.
__construct($A)
QR Decomposition computed by Householder reflections.
getQ()
Generate and return the (economy-sized) orthogonal factor.
isFullRank()
Is the matrix full rank?
solve($B)
Least squares solution of A*X = B.
getH()
Return the Householder vectors.
$i
Definition: disco.tpl.php:19
const MatrixRankException
Definition: Error.php:57
$s
Definition: pwgen.php:45