1<?php
2
3declare(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
32namespace Phpml\Math\LinearAlgebra;
33
34use Phpml\Exception\MatrixException;
35use Phpml\Math\Matrix;
36
37class 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 new MatrixException('Matrix is not square matrix');
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            for ($i = 0; $i < $this->m; ++$i) {
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] ?? 0) > abs($LUcolj[$p] ?? 0)) {
122                    $p = $i;
123                }
124            }
125
126            if ($p != $j) {
127                for ($k = 0; $k < $this->n; ++$k) {
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 *= -1;
137            }
138
139            // Compute multipliers.
140            if (($j < $this->m) && ($this->LU[$j][$j] != 0.0)) {
141                for ($i = $j + 1; $i < $this->m; ++$i) {
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(): array
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        for ($j = 0; $j < $this->n; ++$j) {
220            if ($this->LU[$j][$j] == 0) {
221                return false;
222            }
223        }
224
225        return true;
226    }
227
228    public function det(): float
229    {
230        $d = $this->pivsign;
231        for ($j = 0; $j < $this->n; ++$j) {
232            $d *= $this->LU[$j][$j];
233        }
234
235        return (float) $d;
236    }
237
238    /**
239     * Solve A*X = B
240     *
241     * @param Matrix $B A Matrix with as many rows as A and any number of columns.
242     *
243     * @return array X so that L*U*X = B(piv,:)
244     *
245     * @throws MatrixException
246     */
247    public function solve(Matrix $B): array
248    {
249        if ($B->getRows() != $this->m) {
250            throw new MatrixException('Matrix is not square matrix');
251        }
252
253        if (!$this->isNonsingular()) {
254            throw new MatrixException('Matrix is singular');
255        }
256
257        // Copy right hand side with pivoting
258        $nx = $B->getColumns();
259        $X = $this->getSubMatrix($B->toArray(), $this->piv, 0, $nx - 1);
260        // Solve L*Y = B(piv,:)
261        for ($k = 0; $k < $this->n; ++$k) {
262            for ($i = $k + 1; $i < $this->n; ++$i) {
263                for ($j = 0; $j < $nx; ++$j) {
264                    $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
265                }
266            }
267        }
268
269        // Solve U*X = Y;
270        for ($k = $this->n - 1; $k >= 0; --$k) {
271            for ($j = 0; $j < $nx; ++$j) {
272                $X[$k][$j] /= $this->LU[$k][$k];
273            }
274
275            for ($i = 0; $i < $k; ++$i) {
276                for ($j = 0; $j < $nx; ++$j) {
277                    $X[$i][$j] -= $X[$k][$j] * $this->LU[$i][$k];
278                }
279            }
280        }
281
282        return $X;
283    }
284
285    protected function getSubMatrix(array $matrix, array $RL, int $j0, int $jF): array
286    {
287        $m = count($RL);
288        $n = $jF - $j0;
289        $R = array_fill(0, $m, array_fill(0, $n + 1, 0.0));
290
291        for ($i = 0; $i < $m; ++$i) {
292            for ($j = $j0; $j <= $jF; ++$j) {
293                $R[$i][$j - $j0] = $matrix[$RL[$i]][$j];
294            }
295        }
296
297        return $R;
298    }
299}
300