See Release Notes
Long Term Support Release
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 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 }
title
Description
Body
title
Description
Body
title
Description
Body
title
Body