These results are based on our legacy PHP analysis, consider migrating to our new PHP analysis engine instead. Learn more
1 | <?php |
||
2 | |||
3 | declare(strict_types=1); |
||
4 | |||
5 | /** |
||
6 | * @package JAMA |
||
7 | * |
||
8 | * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n |
||
9 | * unit lower triangular matrix L, an n-by-n upper triangular matrix U, |
||
10 | * and a permutation vector piv of length m so that A(piv,:) = L*U. |
||
11 | * If m < n, then L is m-by-m and U is m-by-n. |
||
12 | * |
||
13 | * The LU decompostion with pivoting always exists, even if the matrix is |
||
14 | * singular, so the constructor will never fail. The primary use of the |
||
15 | * LU decomposition is in the solution of square systems of simultaneous |
||
16 | * linear equations. This will fail if isNonsingular() returns false. |
||
17 | * |
||
18 | * @author Paul Meagher |
||
19 | * @author Bartosz Matosiuk |
||
20 | * @author Michael Bommarito |
||
21 | * |
||
22 | * @version 1.1 |
||
23 | * |
||
24 | * @license PHP v3.0 |
||
25 | * |
||
26 | * Slightly changed to adapt the original code to PHP-ML library |
||
27 | * @date 2017/04/24 |
||
28 | * |
||
29 | * @author Mustafa Karabulut |
||
30 | */ |
||
31 | |||
32 | namespace Phpml\Math\LinearAlgebra; |
||
33 | |||
34 | use Phpml\Exception\MatrixException; |
||
35 | use Phpml\Math\Matrix; |
||
36 | |||
37 | class LUDecomposition |
||
38 | { |
||
39 | /** |
||
40 | * Decomposition storage |
||
41 | * |
||
42 | * @var array |
||
43 | */ |
||
44 | private $LU = []; |
||
45 | |||
46 | /** |
||
47 | * Row dimension. |
||
48 | * |
||
49 | * @var int |
||
50 | */ |
||
51 | private $m; |
||
52 | |||
53 | /** |
||
54 | * Column dimension. |
||
55 | * |
||
56 | * @var int |
||
57 | */ |
||
58 | private $n; |
||
59 | |||
60 | /** |
||
61 | * Pivot sign. |
||
62 | * |
||
63 | * @var int |
||
64 | */ |
||
65 | private $pivsign; |
||
66 | |||
67 | /** |
||
68 | * Internal storage of pivot vector. |
||
69 | * |
||
70 | * @var array |
||
71 | */ |
||
72 | private $piv = []; |
||
73 | |||
74 | /** |
||
75 | * Constructs Structure to access L, U and piv. |
||
76 | * |
||
77 | * @param Matrix $A Rectangular matrix |
||
78 | * |
||
79 | * @throws MatrixException |
||
80 | */ |
||
81 | public function __construct(Matrix $A) |
||
82 | { |
||
83 | if ($A->getRows() != $A->getColumns()) { |
||
84 | throw MatrixException::notSquareMatrix(); |
||
85 | } |
||
86 | |||
87 | // Use a "left-looking", dot-product, Crout/Doolittle algorithm. |
||
88 | $this->LU = $A->toArray(); |
||
89 | $this->m = $A->getRows(); |
||
90 | $this->n = $A->getColumns(); |
||
91 | for ($i = 0; $i < $this->m; ++$i) { |
||
92 | $this->piv[$i] = $i; |
||
93 | } |
||
94 | |||
95 | $this->pivsign = 1; |
||
96 | $LUcolj = []; |
||
97 | |||
98 | // Outer loop. |
||
99 | for ($j = 0; $j < $this->n; ++$j) { |
||
100 | // Make a copy of the j-th column to localize references. |
||
101 | View Code Duplication | for ($i = 0; $i < $this->m; ++$i) { |
|
0 ignored issues
–
show
|
|||
102 | $LUcolj[$i] = &$this->LU[$i][$j]; |
||
103 | } |
||
104 | |||
105 | // Apply previous transformations. |
||
106 | for ($i = 0; $i < $this->m; ++$i) { |
||
107 | $LUrowi = $this->LU[$i]; |
||
108 | // Most of the time is spent in the following dot product. |
||
109 | $kmax = min($i, $j); |
||
110 | $s = 0.0; |
||
111 | for ($k = 0; $k < $kmax; ++$k) { |
||
112 | $s += $LUrowi[$k] * $LUcolj[$k]; |
||
113 | } |
||
114 | |||
115 | $LUrowi[$j] = $LUcolj[$i] -= $s; |
||
116 | } |
||
117 | |||
118 | // Find pivot and exchange if necessary. |
||
119 | $p = $j; |
||
120 | for ($i = $j + 1; $i < $this->m; ++$i) { |
||
121 | if (abs($LUcolj[$i]) > abs($LUcolj[$p])) { |
||
122 | $p = $i; |
||
123 | } |
||
124 | } |
||
125 | |||
126 | if ($p != $j) { |
||
127 | View Code Duplication | for ($k = 0; $k < $this->n; ++$k) { |
|
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
|||
128 | $t = $this->LU[$p][$k]; |
||
129 | $this->LU[$p][$k] = $this->LU[$j][$k]; |
||
130 | $this->LU[$j][$k] = $t; |
||
131 | } |
||
132 | |||
133 | $k = $this->piv[$p]; |
||
134 | $this->piv[$p] = $this->piv[$j]; |
||
135 | $this->piv[$j] = $k; |
||
136 | $this->pivsign = $this->pivsign * -1; |
||
137 | } |
||
138 | |||
139 | // Compute multipliers. |
||
140 | if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) { |
||
141 | View Code Duplication | for ($i = $j + 1; $i < $this->m; ++$i) { |
|
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
|||
142 | $this->LU[$i][$j] /= $this->LU[$j][$j]; |
||
143 | } |
||
144 | } |
||
145 | } |
||
146 | } |
||
147 | |||
148 | /** |
||
149 | * Get lower triangular factor. |
||
150 | * |
||
151 | * @return Matrix Lower triangular factor |
||
152 | */ |
||
153 | public function getL(): Matrix |
||
154 | { |
||
155 | $L = []; |
||
156 | for ($i = 0; $i < $this->m; ++$i) { |
||
157 | for ($j = 0; $j < $this->n; ++$j) { |
||
158 | if ($i > $j) { |
||
159 | $L[$i][$j] = $this->LU[$i][$j]; |
||
160 | } elseif ($i == $j) { |
||
161 | $L[$i][$j] = 1.0; |
||
162 | } else { |
||
163 | $L[$i][$j] = 0.0; |
||
164 | } |
||
165 | } |
||
166 | } |
||
167 | |||
168 | return new Matrix($L); |
||
169 | } |
||
170 | |||
171 | /** |
||
172 | * Get upper triangular factor. |
||
173 | * |
||
174 | * @return Matrix Upper triangular factor |
||
175 | */ |
||
176 | public function getU(): Matrix |
||
177 | { |
||
178 | $U = []; |
||
179 | for ($i = 0; $i < $this->n; ++$i) { |
||
180 | for ($j = 0; $j < $this->n; ++$j) { |
||
181 | if ($i <= $j) { |
||
182 | $U[$i][$j] = $this->LU[$i][$j]; |
||
183 | } else { |
||
184 | $U[$i][$j] = 0.0; |
||
185 | } |
||
186 | } |
||
187 | } |
||
188 | |||
189 | return new Matrix($U); |
||
190 | } |
||
191 | |||
192 | /** |
||
193 | * Return pivot permutation vector. |
||
194 | * |
||
195 | * @return array Pivot vector |
||
196 | */ |
||
197 | public function getPivot(): array |
||
198 | { |
||
199 | return $this->piv; |
||
200 | } |
||
201 | |||
202 | /** |
||
203 | * Alias for getPivot |
||
204 | * |
||
205 | * @see getPivot |
||
206 | */ |
||
207 | public function getDoublePivot() |
||
208 | { |
||
209 | return $this->getPivot(); |
||
210 | } |
||
211 | |||
212 | /** |
||
213 | * Is the matrix nonsingular? |
||
214 | * |
||
215 | * @return bool true if U, and hence A, is nonsingular. |
||
216 | */ |
||
217 | public function isNonsingular(): bool |
||
218 | { |
||
219 | View Code Duplication | for ($j = 0; $j < $this->n; ++$j) { |
|
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
|||
220 | if ($this->LU[$j][$j] == 0) { |
||
221 | return false; |
||
222 | } |
||
223 | } |
||
224 | |||
225 | return true; |
||
226 | } |
||
227 | |||
228 | /** |
||
229 | * Count determinants |
||
230 | * |
||
231 | * @return float|int d matrix determinant |
||
232 | * |
||
233 | * @throws MatrixException |
||
234 | */ |
||
235 | public function det() |
||
236 | { |
||
237 | if ($this->m !== $this->n) { |
||
238 | throw MatrixException::notSquareMatrix(); |
||
239 | } |
||
240 | |||
241 | $d = $this->pivsign; |
||
242 | View Code Duplication | for ($j = 0; $j < $this->n; ++$j) { |
|
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
|||
243 | $d *= $this->LU[$j][$j]; |
||
244 | } |
||
245 | |||
246 | return $d; |
||
247 | } |
||
248 | |||
249 | /** |
||
250 | * Solve A*X = B |
||
251 | * |
||
252 | * @param Matrix $B A Matrix with as many rows as A and any number of columns. |
||
253 | * |
||
254 | * @return array X so that L*U*X = B(piv,:) |
||
255 | * |
||
256 | * @throws MatrixException |
||
257 | */ |
||
258 | public function solve(Matrix $B): array |
||
259 | { |
||
260 | if ($B->getRows() != $this->m) { |
||
261 | throw MatrixException::notSquareMatrix(); |
||
262 | } |
||
263 | |||
264 | if (!$this->isNonsingular()) { |
||
265 | throw MatrixException::singularMatrix(); |
||
266 | } |
||
267 | |||
268 | // Copy right hand side with pivoting |
||
269 | $nx = $B->getColumns(); |
||
270 | $X = $this->getSubMatrix($B->toArray(), $this->piv, 0, $nx - 1); |
||
271 | // Solve L*Y = B(piv,:) |
||
272 | for ($k = 0; $k < $this->n; ++$k) { |
||
273 | View Code Duplication | for ($i = $k + 1; $i < $this->n; ++$i) { |
|
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
|||
274 | for ($j = 0; $j < $nx; ++$j) { |
||
275 | $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k]; |
||
276 | } |
||
277 | } |
||
278 | } |
||
279 | |||
280 | // Solve U*X = Y; |
||
281 | for ($k = $this->n - 1; $k >= 0; --$k) { |
||
282 | for ($j = 0; $j < $nx; ++$j) { |
||
283 | $X[$k][$j] /= $this->LU[$k][$k]; |
||
284 | } |
||
285 | |||
286 | View Code Duplication | for ($i = 0; $i < $k; ++$i) { |
|
0 ignored issues
–
show
This code seems to be duplicated across your project.
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation. You can also find more detailed suggestions in the “Code” section of your repository.
Loading history...
|
|||
287 | for ($j = 0; $j < $nx; ++$j) { |
||
288 | $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k]; |
||
289 | } |
||
290 | } |
||
291 | } |
||
292 | |||
293 | return $X; |
||
294 | } |
||
295 | |||
296 | protected function getSubMatrix(array $matrix, array $RL, int $j0, int $jF): array |
||
297 | { |
||
298 | $m = count($RL); |
||
299 | $n = $jF - $j0; |
||
300 | $R = array_fill(0, $m, array_fill(0, $n + 1, 0.0)); |
||
301 | |||
302 | for ($i = 0; $i < $m; ++$i) { |
||
303 | for ($j = $j0; $j <= $jF; ++$j) { |
||
304 | $R[$i][$j - $j0] = $matrix[$RL[$i]][$j]; |
||
305 | } |
||
306 | } |
||
307 | |||
308 | return $R; |
||
309 | } |
||
310 | } |
||
311 |
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.
You can also find more detailed suggestions in the “Code” section of your repository.