ILIAS  release_5-4 Revision v5.4.26-12-gabc799a52e6
QR.php
Go to the documentation of this file.
1<?php
2
4
7
8class QR
9{
10 private $qrMatrix;
11 private $rows;
12 private $columns;
13
14 private $rDiagonal = [];
15
16 public function __construct(Matrix $matrix)
17 {
18 $this->qrMatrix = $matrix->toArray();
19 $this->rows = $matrix->rows;
20 $this->columns = $matrix->columns;
21
22 $this->decompose();
23 }
24
25 public function getHouseholdVectors(): Matrix
26 {
27 $householdVectors = [];
28 for ($row = 0; $row < $this->rows; ++$row) {
29 for ($column = 0; $column < $this->columns; ++$column) {
30 if ($row >= $column) {
31 $householdVectors[$row][$column] = $this->qrMatrix[$row][$column];
32 } else {
33 $householdVectors[$row][$column] = 0.0;
34 }
35 }
36 }
37
38 return new Matrix($householdVectors);
39 }
40
41 public function getQ(): Matrix
42 {
43 $qGrid = [];
44
45 $rowCount = $this->rows;
46 for ($k = $this->columns - 1; $k >= 0; --$k) {
47 for ($i = 0; $i < $this->rows; ++$i) {
48 $qGrid[$i][$k] = 0.0;
49 }
50 $qGrid[$k][$k] = 1.0;
51 if ($this->columns > $this->rows) {
52 $qGrid = array_slice($qGrid, 0, $this->rows);
53 }
54
55 for ($j = $k; $j < $this->columns; ++$j) {
56 if (isset($this->qrMatrix[$k], $this->qrMatrix[$k][$k]) && $this->qrMatrix[$k][$k] != 0.0) {
57 $s = 0.0;
58 for ($i = $k; $i < $this->rows; ++$i) {
59 $s += $this->qrMatrix[$i][$k] * $qGrid[$i][$j];
60 }
61 $s = -$s / $this->qrMatrix[$k][$k];
62 for ($i = $k; $i < $this->rows; ++$i) {
63 $qGrid[$i][$j] += $s * $this->qrMatrix[$i][$k];
64 }
65 }
66 }
67 }
68
69 array_walk(
70 $qGrid,
71 function (&$row) use ($rowCount) {
72 $row = array_reverse($row);
73 $row = array_slice($row, 0, $rowCount);
74 }
75 );
76
77 return new Matrix($qGrid);
78 }
79
80 public function getR(): Matrix
81 {
82 $rGrid = [];
83
84 for ($row = 0; $row < $this->columns; ++$row) {
85 for ($column = 0; $column < $this->columns; ++$column) {
86 if ($row < $column) {
87 $rGrid[$row][$column] = $this->qrMatrix[$row][$column] ?? 0.0;
88 } elseif ($row === $column) {
89 $rGrid[$row][$column] = $this->rDiagonal[$row] ?? 0.0;
90 } else {
91 $rGrid[$row][$column] = 0.0;
92 }
93 }
94 }
95
96 if ($this->columns > $this->rows) {
97 $rGrid = array_slice($rGrid, 0, $this->rows);
98 }
99
100 return new Matrix($rGrid);
101 }
102
103 private function hypo($a, $b): float
104 {
105 if (abs($a) > abs($b)) {
106 $r = $b / $a;
107 $r = abs($a) * sqrt(1 + $r * $r);
108 } elseif ($b != 0.0) {
109 $r = $a / $b;
110 $r = abs($b) * sqrt(1 + $r * $r);
111 } else {
112 $r = 0.0;
113 }
114
115 return $r;
116 }
117
121 private function decompose(): void
122 {
123 for ($k = 0; $k < $this->columns; ++$k) {
124 // Compute 2-norm of k-th column without under/overflow.
125 $norm = 0.0;
126 for ($i = $k; $i < $this->rows; ++$i) {
127 $norm = $this->hypo($norm, $this->qrMatrix[$i][$k]);
128 }
129 if ($norm != 0.0) {
130 // Form k-th Householder vector.
131 if ($this->qrMatrix[$k][$k] < 0.0) {
132 $norm = -$norm;
133 }
134 for ($i = $k; $i < $this->rows; ++$i) {
135 $this->qrMatrix[$i][$k] /= $norm;
136 }
137 $this->qrMatrix[$k][$k] += 1.0;
138 // Apply transformation to remaining columns.
139 for ($j = $k + 1; $j < $this->columns; ++$j) {
140 $s = 0.0;
141 for ($i = $k; $i < $this->rows; ++$i) {
142 $s += $this->qrMatrix[$i][$k] * $this->qrMatrix[$i][$j];
143 }
144 $s = -$s / $this->qrMatrix[$k][$k];
145 for ($i = $k; $i < $this->rows; ++$i) {
146 $this->qrMatrix[$i][$j] += $s * $this->qrMatrix[$i][$k];
147 }
148 }
149 }
150 $this->rDiagonal[$k] = -$norm;
151 }
152 }
153
154 public function isFullRank(): bool
155 {
156 for ($j = 0; $j < $this->columns; ++$j) {
157 if ($this->rDiagonal[$j] == 0.0) {
158 return false;
159 }
160 }
161
162 return true;
163 }
164
174 public function solve(Matrix $B): Matrix
175 {
176 if ($B->rows !== $this->rows) {
177 throw new Exception('Matrix row dimensions are not equal');
178 }
179
180 if (!$this->isFullRank()) {
181 throw new Exception('Can only perform this operation on a full-rank matrix');
182 }
183
184 // Compute Y = transpose(Q)*B
185 $Y = $this->getQ()->transpose()
186 ->multiply($B);
187 // Solve R*X = Y;
188 return $this->getR()->inverse()
189 ->multiply($Y);
190 }
191}
An exception for terminatinating execution or to throw for unit testing.
decompose()
QR Decomposition computed by Householder reflections.
Definition: QR.php:121
solve(Matrix $B)
Least squares solution of A*X = B.
Definition: QR.php:174
__construct(Matrix $matrix)
Definition: QR.php:16
rows()
Returns a Generator that will yield each row of the matrix in turn as a vector matrix or the value of...
Definition: Matrix.php:282
$i
Definition: disco.tpl.php:19
$r
Definition: example_031.php:79
$matrix
Definition: test.php:18
$row
Class for the creating "special" Matrices.
Definition: Builder.php:11
$s
Definition: pwgen.php:45