[Summary view]

   1  <?php
   2  
   3  namespace Matrix\Decomposition;
   4  
   5  use Matrix\Exception;
   6  use Matrix\Matrix;
   7  
   8  class LU
   9  {
  10      private $luMatrix;
  11      private $rows;
  12      private $columns;
  13  
  14      private $pivot = [];
  15  
  16      public function __construct(Matrix $matrix)
  17      {
  18          $this->luMatrix = $matrix->toArray();
  19          $this->rows = $matrix->rows;
  20          $this->columns = $matrix->columns;
  21  
  22          $this->buildPivot();
  23      }
  24  
  25      /**
  26       * Get lower triangular factor.
  27       *
  28       * @return Matrix Lower triangular factor
  29       */
  30      public function getL(): Matrix
  31      {
  32          $lower = [];
  33  
  34          $columns = min($this->rows, $this->columns);
  35          for ($row = 0; $row < $this->rows; ++$row) {
  36              for ($column = 0; $column < $columns; ++$column) {
  37                  if ($row > $column) {
  38                      $lower[$row][$column] = $this->luMatrix[$row][$column];
  39                  } elseif ($row === $column) {
  40                      $lower[$row][$column] = 1.0;
  41                  } else {
  42                      $lower[$row][$column] = 0.0;
  43                  }
  44              }
  45          }
  46  
  47          return new Matrix($lower);
  48      }
  49  
  50      /**
  51       * Get upper triangular factor.
  52       *
  53       * @return Matrix Upper triangular factor
  54       */
  55      public function getU(): Matrix
  56      {
  57          $upper = [];
  58  
  59          $rows = min($this->rows, $this->columns);
  60          for ($row = 0; $row < $rows; ++$row) {
  61              for ($column = 0; $column < $this->columns; ++$column) {
  62                  if ($row <= $column) {
  63                      $upper[$row][$column] = $this->luMatrix[$row][$column];
  64                  } else {
  65                      $upper[$row][$column] = 0.0;
  66                  }
  67              }
  68          }
  69  
  70          return new Matrix($upper);
  71      }
  72  
  73      /**
  74       * Return pivot permutation vector.
  75       *
  76       * @return Matrix Pivot matrix
  77       */
  78      public function getP(): Matrix
  79      {
  80          $pMatrix = [];
  81  
  82          $pivots = $this->pivot;
  83          $pivotCount = count($pivots);
  84          foreach ($pivots as $row => $pivot) {
  85              $pMatrix[$row] = array_fill(0, $pivotCount, 0);
  86              $pMatrix[$row][$pivot] = 1;
  87          }
  88  
  89          return new Matrix($pMatrix);
  90      }
  91  
  92      /**
  93       * Return pivot permutation vector.
  94       *
  95       * @return array Pivot vector
  96       */
  97      public function getPivot(): array
  98      {
  99          return $this->pivot;
 100      }
 101  
 102      /**
 103       *    Is the matrix nonsingular?
 104       *
 105       * @return bool true if U, and hence A, is nonsingular
 106       */
 107      public function isNonsingular(): bool
 108      {
 109          for ($diagonal = 0; $diagonal < $this->columns; ++$diagonal) {
 110              if ($this->luMatrix[$diagonal][$diagonal] === 0.0) {
 111                  return false;
 112              }
 113          }
 114  
 115          return true;
 116      }
 117  
 118      private function buildPivot(): void
 119      {
 120          for ($row = 0; $row < $this->rows; ++$row) {
 121              $this->pivot[$row] = $row;
 122          }
 123  
 124          for ($column = 0; $column < $this->columns; ++$column) {
 125              $luColumn = $this->localisedReferenceColumn($column);
 126  
 127              $this->applyTransformations($column, $luColumn);
 128  
 129              $pivot = $this->findPivot($column, $luColumn);
 130              if ($pivot !== $column) {
 131                  $this->pivotExchange($pivot, $column);
 132              }
 133  
 134              $this->computeMultipliers($column);
 135  
 136              unset($luColumn);
 137          }
 138      }
 139  
 140      private function localisedReferenceColumn($column): array
 141      {
 142          $luColumn = [];
 143  
 144          for ($row = 0; $row < $this->rows; ++$row) {
 145              $luColumn[$row] = &$this->luMatrix[$row][$column];
 146          }
 147  
 148          return $luColumn;
 149      }
 150  
 151      private function applyTransformations($column, array $luColumn): void
 152      {
 153          for ($row = 0; $row < $this->rows; ++$row) {
 154              $luRow = $this->luMatrix[$row];
 155              // Most of the time is spent in the following dot product.
 156              $kmax = min($row, $column);
 157              $sValue = 0.0;
 158              for ($kValue = 0; $kValue < $kmax; ++$kValue) {
 159                  $sValue += $luRow[$kValue] * $luColumn[$kValue];
 160              }
 161              $luRow[$column] = $luColumn[$row] -= $sValue;
 162          }
 163      }
 164  
 165      private function findPivot($column, array $luColumn): int
 166      {
 167          $pivot = $column;
 168          for ($row = $column + 1; $row < $this->rows; ++$row) {
 169              if (abs($luColumn[$row]) > abs($luColumn[$pivot])) {
 170                  $pivot = $row;
 171              }
 172          }
 173  
 174          return $pivot;
 175      }
 176  
 177      private function pivotExchange($pivot, $column): void
 178      {
 179          for ($kValue = 0; $kValue < $this->columns; ++$kValue) {
 180              $tValue = $this->luMatrix[$pivot][$kValue];
 181              $this->luMatrix[$pivot][$kValue] = $this->luMatrix[$column][$kValue];
 182              $this->luMatrix[$column][$kValue] = $tValue;
 183          }
 184  
 185          $lValue = $this->pivot[$pivot];
 186          $this->pivot[$pivot] = $this->pivot[$column];
 187          $this->pivot[$column] = $lValue;
 188      }
 189  
 190      private function computeMultipliers($diagonal): void
 191      {
 192          if (($diagonal < $this->rows) && ($this->luMatrix[$diagonal][$diagonal] != 0.0)) {
 193              for ($row = $diagonal + 1; $row < $this->rows; ++$row) {
 194                  $this->luMatrix[$row][$diagonal] /= $this->luMatrix[$diagonal][$diagonal];
 195              }
 196          }
 197      }
 198  
 199      private function pivotB(Matrix $B): array
 200      {
 201          $X = [];
 202          foreach ($this->pivot as $rowId) {
 203              $row = $B->getRows($rowId + 1)->toArray();
 204              $X[] = array_pop($row);
 205          }
 206  
 207          return $X;
 208      }
 209  
 210      /**
 211       * Solve A*X = B.
 212       *
 213       * @param Matrix $B a Matrix with as many rows as A and any number of columns
 214       *
 215       * @throws Exception
 216       *
 217       * @return Matrix X so that L*U*X = B(piv,:)
 218       */
 219      public function solve(Matrix $B): Matrix
 220      {
 221          if ($B->rows !== $this->rows) {
 222              throw new Exception('Matrix row dimensions are not equal');
 223          }
 224  
 225          if ($this->rows !== $this->columns) {
 226              throw new Exception('LU solve() only works on square matrices');
 227          }
 228  
 229          if (!$this->isNonsingular()) {
 230              throw new Exception('Can only perform operation on singular matrix');
 231          }
 232  
 233          // Copy right hand side with pivoting
 234          $nx = $B->columns;
 235          $X = $this->pivotB($B);
 236  
 237          // Solve L*Y = B(piv,:)
 238          for ($k = 0; $k < $this->columns; ++$k) {
 239              for ($i = $k + 1; $i < $this->columns; ++$i) {
 240                  for ($j = 0; $j < $nx; ++$j) {
 241                      $X[$i][$j] -= $X[$k][$j] * $this->luMatrix[$i][$k];
 242                  }
 243              }
 244          }
 245  
 246          // Solve U*X = Y;
 247          for ($k = $this->columns - 1; $k >= 0; --$k) {
 248              for ($j = 0; $j < $nx; ++$j) {
 249                  $X[$k][$j] /= $this->luMatrix[$k][$k];
 250              }
 251              for ($i = 0; $i < $k; ++$i) {
 252                  for ($j = 0; $j < $nx; ++$j) {
 253                      $X[$i][$j] -= $X[$k][$j] * $this->luMatrix[$i][$k];
 254                  }
 255              }
 256          }
 257  
 258          return new Matrix($X);
 259      }
 260  }