tablepress / libraries / vendor / markbaker / matrix / classes / Decomposition / QR.php

QR.php in TablePress – Tables in WordPress made easy 3.2.1, at libraries/vendor/markbaker/matrix/classes/Decomposition/QR.php

192 lines 4.3 KB

Raw Download Zip

1	<?php
2
3	namespace TablePress\Matrix\Decomposition;
4
5	use TablePress\Matrix\Exception;
6	use TablePress\Matrix\Matrix;
7
8	class 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
118	/**
119	* QR Decomposition computed by Householder reflections.
120	*/
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
165	/**
166	* Least squares solution of A*X = B.
167	*
168	* @param Matrix $B a Matrix with as many rows as A and any number of columns
169	*
170	* @throws Exception
171	*
172	* @return Matrix matrix that minimizes the two norm of QRX-B
173	*/
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	}
192