1
|
|
|
<?php |
2
|
|
|
|
3
|
|
|
declare(strict_types=1); |
4
|
|
|
|
5
|
|
|
namespace Phpml\Math; |
6
|
|
|
|
7
|
|
|
use Phpml\Exception\InvalidArgumentException; |
8
|
|
|
use Phpml\Exception\MatrixException; |
9
|
|
|
use Phpml\Math\LinearAlgebra\LUDecomposition; |
10
|
|
|
|
11
|
|
|
class Matrix |
12
|
|
|
{ |
13
|
|
|
/** |
14
|
|
|
* @var array |
15
|
|
|
*/ |
16
|
|
|
private $matrix = []; |
17
|
|
|
|
18
|
|
|
/** |
19
|
|
|
* @var int |
20
|
|
|
*/ |
21
|
|
|
private $rows; |
22
|
|
|
|
23
|
|
|
/** |
24
|
|
|
* @var int |
25
|
|
|
*/ |
26
|
|
|
private $columns; |
27
|
|
|
|
28
|
|
|
/** |
29
|
|
|
* @var float |
30
|
|
|
*/ |
31
|
|
|
private $determinant; |
32
|
|
|
|
33
|
|
|
/** |
34
|
|
|
* @throws InvalidArgumentException |
35
|
|
|
*/ |
36
|
|
|
public function __construct(array $matrix, bool $validate = true) |
37
|
|
|
{ |
38
|
|
|
// When a row vector is given |
39
|
|
|
if (!is_array($matrix[0])) { |
40
|
|
|
$this->rows = 1; |
41
|
|
|
$this->columns = count($matrix); |
42
|
|
|
$matrix = [$matrix]; |
43
|
|
|
} else { |
44
|
|
|
$this->rows = count($matrix); |
45
|
|
|
$this->columns = count($matrix[0]); |
46
|
|
|
} |
47
|
|
|
|
48
|
|
|
if ($validate) { |
49
|
|
|
for ($i = 0; $i < $this->rows; ++$i) { |
50
|
|
|
if (count($matrix[$i]) !== $this->columns) { |
51
|
|
|
throw new InvalidArgumentException('Matrix dimensions did not match'); |
52
|
|
|
} |
53
|
|
|
} |
54
|
|
|
} |
55
|
|
|
|
56
|
|
|
$this->matrix = $matrix; |
57
|
|
|
} |
58
|
|
|
|
59
|
|
|
public static function fromFlatArray(array $array): self |
60
|
|
|
{ |
61
|
|
|
$matrix = []; |
62
|
|
|
foreach ($array as $value) { |
63
|
|
|
$matrix[] = [$value]; |
64
|
|
|
} |
65
|
|
|
|
66
|
|
|
return new self($matrix); |
67
|
|
|
} |
68
|
|
|
|
69
|
|
|
public function toArray(): array |
70
|
|
|
{ |
71
|
|
|
return $this->matrix; |
72
|
|
|
} |
73
|
|
|
|
74
|
|
|
public function toScalar(): float |
75
|
|
|
{ |
76
|
|
|
return $this->matrix[0][0]; |
77
|
|
|
} |
78
|
|
|
|
79
|
|
|
public function getRows(): int |
80
|
|
|
{ |
81
|
|
|
return $this->rows; |
82
|
|
|
} |
83
|
|
|
|
84
|
|
|
public function getColumns(): int |
85
|
|
|
{ |
86
|
|
|
return $this->columns; |
87
|
|
|
} |
88
|
|
|
|
89
|
|
|
/** |
90
|
|
|
* @throws MatrixException |
91
|
|
|
*/ |
92
|
|
|
public function getColumnValues(int $column): array |
93
|
|
|
{ |
94
|
|
|
if ($column >= $this->columns) { |
95
|
|
|
throw new MatrixException('Column out of range'); |
96
|
|
|
} |
97
|
|
|
|
98
|
|
|
return array_column($this->matrix, $column); |
99
|
|
|
} |
100
|
|
|
|
101
|
|
|
/** |
102
|
|
|
* @return float|int |
103
|
|
|
* |
104
|
|
|
* @throws MatrixException |
105
|
|
|
*/ |
106
|
|
|
public function getDeterminant() |
107
|
|
|
{ |
108
|
|
|
if ($this->determinant !== null) { |
109
|
|
|
return $this->determinant; |
110
|
|
|
} |
111
|
|
|
|
112
|
|
|
if (!$this->isSquare()) { |
113
|
|
|
throw new MatrixException('Matrix is not square matrix'); |
114
|
|
|
} |
115
|
|
|
|
116
|
|
|
$lu = new LUDecomposition($this); |
117
|
|
|
|
118
|
|
|
return $this->determinant = $lu->det(); |
119
|
|
|
} |
120
|
|
|
|
121
|
|
|
public function isSquare(): bool |
122
|
|
|
{ |
123
|
|
|
return $this->columns === $this->rows; |
124
|
|
|
} |
125
|
|
|
|
126
|
|
|
public function transpose(): self |
127
|
|
|
{ |
128
|
|
|
if ($this->rows === 1) { |
129
|
|
|
$matrix = array_map(static function ($el): array { |
130
|
|
|
return [$el]; |
131
|
|
|
}, $this->matrix[0]); |
132
|
|
|
} else { |
133
|
|
|
$matrix = array_map(null, ...$this->matrix); |
134
|
|
|
} |
135
|
|
|
|
136
|
|
|
return new self($matrix, false); |
137
|
|
|
} |
138
|
|
|
|
139
|
|
|
public function multiply(self $matrix): self |
140
|
|
|
{ |
141
|
|
|
if ($this->columns !== $matrix->getRows()) { |
142
|
|
|
throw new InvalidArgumentException('Inconsistent matrix supplied'); |
143
|
|
|
} |
144
|
|
|
|
145
|
|
|
$array1 = $this->toArray(); |
146
|
|
|
$array2 = $matrix->toArray(); |
147
|
|
|
$colCount = $matrix->columns; |
148
|
|
|
|
149
|
|
|
/* |
150
|
|
|
- To speed-up multiplication, we need to avoid use of array index operator [ ] as much as possible( See #255 for details) |
151
|
|
|
- A combination of "foreach" and "array_column" works much faster then accessing the array via index operator |
152
|
|
|
*/ |
153
|
|
|
$product = []; |
154
|
|
|
foreach ($array1 as $row => $rowData) { |
155
|
|
|
for ($col = 0; $col < $colCount; ++$col) { |
156
|
|
|
$columnData = array_column($array2, $col); |
157
|
|
|
$sum = 0; |
158
|
|
|
foreach ($rowData as $key => $valueData) { |
159
|
|
|
$sum += $valueData * $columnData[$key]; |
160
|
|
|
} |
161
|
|
|
|
162
|
|
|
$product[$row][$col] = $sum; |
163
|
|
|
} |
164
|
|
|
} |
165
|
|
|
|
166
|
|
|
return new self($product, false); |
167
|
|
|
} |
168
|
|
|
|
169
|
|
|
/** |
170
|
|
|
* @param float|int $value |
171
|
|
|
*/ |
172
|
|
|
public function divideByScalar($value): self |
173
|
|
|
{ |
174
|
|
|
$newMatrix = []; |
175
|
|
|
for ($i = 0; $i < $this->rows; ++$i) { |
176
|
|
|
for ($j = 0; $j < $this->columns; ++$j) { |
177
|
|
|
$newMatrix[$i][$j] = $this->matrix[$i][$j] / $value; |
178
|
|
|
} |
179
|
|
|
} |
180
|
|
|
|
181
|
|
|
return new self($newMatrix, false); |
182
|
|
|
} |
183
|
|
|
|
184
|
|
|
/** |
185
|
|
|
* @param float|int $value |
186
|
|
|
*/ |
187
|
|
|
public function multiplyByScalar($value): self |
188
|
|
|
{ |
189
|
|
|
$newMatrix = []; |
190
|
|
|
for ($i = 0; $i < $this->rows; ++$i) { |
191
|
|
|
for ($j = 0; $j < $this->columns; ++$j) { |
192
|
|
|
$newMatrix[$i][$j] = $this->matrix[$i][$j] * $value; |
193
|
|
|
} |
194
|
|
|
} |
195
|
|
|
|
196
|
|
|
return new self($newMatrix, false); |
197
|
|
|
} |
198
|
|
|
|
199
|
|
|
/** |
200
|
|
|
* Element-wise addition of the matrix with another one |
201
|
|
|
*/ |
202
|
|
|
public function add(self $other): self |
203
|
|
|
{ |
204
|
|
|
return $this->sum($other); |
205
|
|
|
} |
206
|
|
|
|
207
|
|
|
/** |
208
|
|
|
* Element-wise subtracting of another matrix from this one |
209
|
|
|
*/ |
210
|
|
|
public function subtract(self $other): self |
211
|
|
|
{ |
212
|
|
|
return $this->sum($other, -1); |
213
|
|
|
} |
214
|
|
|
|
215
|
|
|
public function inverse(): self |
216
|
|
|
{ |
217
|
|
|
if (!$this->isSquare()) { |
218
|
|
|
throw new MatrixException('Matrix is not square matrix'); |
219
|
|
|
} |
220
|
|
|
|
221
|
|
|
$LU = new LUDecomposition($this); |
222
|
|
|
$identity = $this->getIdentity(); |
223
|
|
|
$inverse = $LU->solve($identity); |
224
|
|
|
|
225
|
|
|
return new self($inverse, false); |
226
|
|
|
} |
227
|
|
|
|
228
|
|
|
public function crossOut(int $row, int $column): self |
229
|
|
|
{ |
230
|
|
|
$newMatrix = []; |
231
|
|
|
$r = 0; |
232
|
|
|
for ($i = 0; $i < $this->rows; ++$i) { |
233
|
|
|
$c = 0; |
234
|
|
|
if ($row != $i) { |
235
|
|
|
for ($j = 0; $j < $this->columns; ++$j) { |
236
|
|
|
if ($column != $j) { |
237
|
|
|
$newMatrix[$r][$c] = $this->matrix[$i][$j]; |
238
|
|
|
++$c; |
239
|
|
|
} |
240
|
|
|
} |
241
|
|
|
|
242
|
|
|
++$r; |
243
|
|
|
} |
244
|
|
|
} |
245
|
|
|
|
246
|
|
|
return new self($newMatrix, false); |
247
|
|
|
} |
248
|
|
|
|
249
|
|
|
public function isSingular(): bool |
250
|
|
|
{ |
251
|
|
|
return $this->getDeterminant() == 0; |
252
|
|
|
} |
253
|
|
|
|
254
|
|
|
/** |
255
|
|
|
* Frobenius norm (Hilbert–Schmidt norm, Euclidean norm) (‖A‖F) |
256
|
|
|
* Square root of the sum of the square of all elements. |
257
|
|
|
* |
258
|
|
|
* https://en.wikipedia.org/wiki/Matrix_norm#Frobenius_norm |
259
|
|
|
* |
260
|
|
|
* _____________ |
261
|
|
|
* /ᵐ ⁿ |
262
|
|
|
* ‖A‖F = √ Σ Σ |aᵢⱼ|² |
263
|
|
|
* ᵢ₌₁ ᵢ₌₁ |
264
|
|
|
*/ |
265
|
|
|
public function frobeniusNorm(): float |
266
|
|
|
{ |
267
|
|
|
$squareSum = 0; |
268
|
|
|
for ($i = 0; $i < $this->rows; ++$i) { |
269
|
|
|
for ($j = 0; $j < $this->columns; ++$j) { |
270
|
|
|
$squareSum += $this->matrix[$i][$j] ** 2; |
271
|
|
|
} |
272
|
|
|
} |
273
|
|
|
|
274
|
|
|
return $squareSum ** .5; |
275
|
|
|
} |
276
|
|
|
|
277
|
|
|
/** |
278
|
|
|
* Returns the transpose of given array |
279
|
|
|
*/ |
280
|
|
|
public static function transposeArray(array $array): array |
281
|
|
|
{ |
282
|
|
|
return (new self($array, false))->transpose()->toArray(); |
283
|
|
|
} |
284
|
|
|
|
285
|
|
|
/** |
286
|
|
|
* Returns the dot product of two arrays<br> |
287
|
|
|
* Matrix::dot(x, y) ==> x.y' |
288
|
|
|
*/ |
289
|
|
|
public static function dot(array $array1, array $array2): array |
290
|
|
|
{ |
291
|
|
|
$m1 = new self($array1, false); |
292
|
|
|
$m2 = new self($array2, false); |
293
|
|
|
|
294
|
|
|
return $m1->multiply($m2->transpose())->toArray()[0]; |
295
|
|
|
} |
296
|
|
|
|
297
|
|
|
/** |
298
|
|
|
* Element-wise addition or substraction depending on the given sign parameter |
299
|
|
|
*/ |
300
|
|
|
private function sum(self $other, int $sign = 1): self |
301
|
|
|
{ |
302
|
|
|
$a1 = $this->toArray(); |
303
|
|
|
$a2 = $other->toArray(); |
304
|
|
|
|
305
|
|
|
$newMatrix = []; |
306
|
|
|
for ($i = 0; $i < $this->rows; ++$i) { |
307
|
|
|
for ($k = 0; $k < $this->columns; ++$k) { |
308
|
|
|
$newMatrix[$i][$k] = $a1[$i][$k] + $sign * $a2[$i][$k]; |
309
|
|
|
} |
310
|
|
|
} |
311
|
|
|
|
312
|
|
|
return new self($newMatrix, false); |
313
|
|
|
} |
314
|
|
|
|
315
|
|
|
/** |
316
|
|
|
* Returns diagonal identity matrix of the same size of this matrix |
317
|
|
|
*/ |
318
|
|
|
private function getIdentity(): self |
319
|
|
|
{ |
320
|
|
|
$array = array_fill(0, $this->rows, array_fill(0, $this->columns, 0)); |
321
|
|
|
for ($i = 0; $i < $this->rows; ++$i) { |
322
|
|
|
$array[$i][$i] = 1; |
323
|
|
|
} |
324
|
|
|
|
325
|
|
|
return new self($array, false); |
326
|
|
|
} |
327
|
|
|
} |
328
|
|
|
|