ILIAS  release_5-4 Revision v5.4.26-12-gabc799a52e6
QRDecomposition.php
Go to the documentation of this file.
1<?php
2
4
5use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
6
23{
24 const MATRIX_RANK_EXCEPTION = 'Can only perform operation on full-rank matrix.';
25
31 private $QR = [];
32
38 private $m;
39
45 private $n;
46
52 private $Rdiag = [];
53
59 public function __construct(Matrix $A)
60 {
61 // Initialize.
62 $this->QR = $A->getArray();
63 $this->m = $A->getRowDimension();
64 $this->n = $A->getColumnDimension();
65 // Main loop.
66 for ($k = 0; $k < $this->n; ++$k) {
67 // Compute 2-norm of k-th column without under/overflow.
68 $nrm = 0.0;
69 for ($i = $k; $i < $this->m; ++$i) {
70 $nrm = hypo($nrm, $this->QR[$i][$k]);
71 }
72 if ($nrm != 0.0) {
73 // Form k-th Householder vector.
74 if ($this->QR[$k][$k] < 0) {
75 $nrm = -$nrm;
76 }
77 for ($i = $k; $i < $this->m; ++$i) {
78 $this->QR[$i][$k] /= $nrm;
79 }
80 $this->QR[$k][$k] += 1.0;
81 // Apply transformation to remaining columns.
82 for ($j = $k + 1; $j < $this->n; ++$j) {
83 $s = 0.0;
84 for ($i = $k; $i < $this->m; ++$i) {
85 $s += $this->QR[$i][$k] * $this->QR[$i][$j];
86 }
87 $s = -$s / $this->QR[$k][$k];
88 for ($i = $k; $i < $this->m; ++$i) {
89 $this->QR[$i][$j] += $s * $this->QR[$i][$k];
90 }
91 }
92 }
93 $this->Rdiag[$k] = -$nrm;
94 }
95 }
96
97 // function __construct()
98
104 public function isFullRank()
105 {
106 for ($j = 0; $j < $this->n; ++$j) {
107 if ($this->Rdiag[$j] == 0) {
108 return false;
109 }
110 }
111
112 return true;
113 }
114
115 // function isFullRank()
116
122 public function getH()
123 {
124 $H = [];
125 for ($i = 0; $i < $this->m; ++$i) {
126 for ($j = 0; $j < $this->n; ++$j) {
127 if ($i >= $j) {
128 $H[$i][$j] = $this->QR[$i][$j];
129 } else {
130 $H[$i][$j] = 0.0;
131 }
132 }
133 }
134
135 return new Matrix($H);
136 }
137
138 // function getH()
139
145 public function getR()
146 {
147 $R = [];
148 for ($i = 0; $i < $this->n; ++$i) {
149 for ($j = 0; $j < $this->n; ++$j) {
150 if ($i < $j) {
151 $R[$i][$j] = $this->QR[$i][$j];
152 } elseif ($i == $j) {
153 $R[$i][$j] = $this->Rdiag[$i];
154 } else {
155 $R[$i][$j] = 0.0;
156 }
157 }
158 }
159
160 return new Matrix($R);
161 }
162
163 // function getR()
164
170 public function getQ()
171 {
172 $Q = [];
173 for ($k = $this->n - 1; $k >= 0; --$k) {
174 for ($i = 0; $i < $this->m; ++$i) {
175 $Q[$i][$k] = 0.0;
176 }
177 $Q[$k][$k] = 1.0;
178 for ($j = $k; $j < $this->n; ++$j) {
179 if ($this->QR[$k][$k] != 0) {
180 $s = 0.0;
181 for ($i = $k; $i < $this->m; ++$i) {
182 $s += $this->QR[$i][$k] * $Q[$i][$j];
183 }
184 $s = -$s / $this->QR[$k][$k];
185 for ($i = $k; $i < $this->m; ++$i) {
186 $Q[$i][$j] += $s * $this->QR[$i][$k];
187 }
188 }
189 }
190 }
191
192 return new Matrix($Q);
193 }
194
195 // function getQ()
196
204 public function solve(Matrix $B)
205 {
206 if ($B->getRowDimension() == $this->m) {
207 if ($this->isFullRank()) {
208 // Copy right hand side
209 $nx = $B->getColumnDimension();
210 $X = $B->getArray();
211 // Compute Y = transpose(Q)*B
212 for ($k = 0; $k < $this->n; ++$k) {
213 for ($j = 0; $j < $nx; ++$j) {
214 $s = 0.0;
215 for ($i = $k; $i < $this->m; ++$i) {
216 $s += $this->QR[$i][$k] * $X[$i][$j];
217 }
218 $s = -$s / $this->QR[$k][$k];
219 for ($i = $k; $i < $this->m; ++$i) {
220 $X[$i][$j] += $s * $this->QR[$i][$k];
221 }
222 }
223 }
224 // Solve R*X = Y;
225 for ($k = $this->n - 1; $k >= 0; --$k) {
226 for ($j = 0; $j < $nx; ++$j) {
227 $X[$k][$j] /= $this->Rdiag[$k];
228 }
229 for ($i = 0; $i < $k; ++$i) {
230 for ($j = 0; $j < $nx; ++$j) {
231 $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
232 }
233 }
234 }
235 $X = new Matrix($X);
236
237 return $X->getMatrix(0, $this->n - 1, 0, $nx);
238 }
239
240 throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
241 }
242
243 throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
244 }
245}
hypo($a, $b)
Pythagorean Theorem:.
Definition: Maths.php:18
An exception for terminatinating execution or to throw for unit testing.
getColumnDimension()
getColumnDimension.
Definition: Matrix.php:137
For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n orthogonal matrix Q and an n-by...
__construct(Matrix $A)
QR Decomposition computed by Householder reflections.
solve(Matrix $B)
Least squares solution of A*X = B.
getQ()
Generate and return the (economy-sized) orthogonal factor.
$i
Definition: disco.tpl.php:19
$X
Definition: test.php:23
Class for the creating "special" Matrices.
Definition: Builder.php:11
$s
Definition: pwgen.php:45