[Summary view]

   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  }