[Summary view]

   1  <?php
   2  
   3  namespace Matrix\Decomposition;
   4  
   5  use Matrix\Exception;
   6  use Matrix\Matrix;
   7  
   8  class QR
   9  {
  10      private $qrMatrix;
  11      private $rows;
  12      private $columns;
  13  
  14      private $rDiagonal = [];
  15  
  16      public function __construct(Matrix $matrix)
  17      {
  18          $this->qrMatrix = $matrix->toArray();
  19          $this->rows = $matrix->rows;
  20          $this->columns = $matrix->columns;
  21  
  22          $this->decompose();
  23      }
  24  
  25      public function getHouseholdVectors(): Matrix
  26      {
  27          $householdVectors = [];
  28          for ($row = 0; $row < $this->rows; ++$row) {
  29              for ($column = 0; $column < $this->columns; ++$column) {
  30                  if ($row >= $column) {
  31                      $householdVectors[$row][$column] = $this->qrMatrix[$row][$column];
  32                  } else {
  33                      $householdVectors[$row][$column] = 0.0;
  34                  }
  35              }
  36          }
  37  
  38          return new Matrix($householdVectors);
  39      }
  40  
  41      public function getQ(): Matrix
  42      {
  43          $qGrid = [];
  44  
  45          $rowCount = $this->rows;
  46          for ($k = $this->columns - 1; $k >= 0; --$k) {
  47              for ($i = 0; $i < $this->rows; ++$i) {
  48                  $qGrid[$i][$k] = 0.0;
  49              }
  50              $qGrid[$k][$k] = 1.0;
  51              if ($this->columns > $this->rows) {
  52                  $qGrid = array_slice($qGrid, 0, $this->rows);
  53              }
  54  
  55              for ($j = $k; $j < $this->columns; ++$j) {
  56                  if (isset($this->qrMatrix[$k], $this->qrMatrix[$k][$k]) && $this->qrMatrix[$k][$k] != 0.0) {
  57                      $s = 0.0;
  58                      for ($i = $k; $i < $this->rows; ++$i) {
  59                          $s += $this->qrMatrix[$i][$k] * $qGrid[$i][$j];
  60                      }
  61                      $s = -$s / $this->qrMatrix[$k][$k];
  62                      for ($i = $k; $i < $this->rows; ++$i) {
  63                          $qGrid[$i][$j] += $s * $this->qrMatrix[$i][$k];
  64                      }
  65                  }
  66              }
  67          }
  68  
  69          array_walk(
  70              $qGrid,
  71              function (&$row) use ($rowCount) {
  72                  $row = array_reverse($row);
  73                  $row = array_slice($row, 0, $rowCount);
  74              }
  75          );
  76  
  77          return new Matrix($qGrid);
  78      }
  79  
  80      public function getR(): Matrix
  81      {
  82          $rGrid = [];
  83  
  84          for ($row = 0; $row < $this->columns; ++$row) {
  85              for ($column = 0; $column < $this->columns; ++$column) {
  86                  if ($row < $column) {
  87                      $rGrid[$row][$column] = $this->qrMatrix[$row][$column] ?? 0.0;
  88                  } elseif ($row === $column) {
  89                      $rGrid[$row][$column] = $this->rDiagonal[$row] ?? 0.0;
  90                  } else {
  91                      $rGrid[$row][$column] = 0.0;
  92                  }
  93              }
  94          }
  95  
  96          if ($this->columns > $this->rows) {
  97              $rGrid = array_slice($rGrid, 0, $this->rows);
  98          }
  99  
 100          return new Matrix($rGrid);
 101      }
 102  
 103      private function hypo($a, $b): float
 104      {
 105          if (abs($a) > abs($b)) {
 106              $r = $b / $a;
 107              $r = abs($a) * sqrt(1 + $r * $r);
 108          } elseif ($b != 0.0) {
 109              $r = $a / $b;
 110              $r = abs($b) * sqrt(1 + $r * $r);
 111          } else {
 112              $r = 0.0;
 113          }
 114  
 115          return $r;
 116      }
 117  
 118      /**
 119       * QR Decomposition computed by Householder reflections.
 120       */
 121      private function decompose(): void
 122      {
 123          for ($k = 0; $k < $this->columns; ++$k) {
 124              // Compute 2-norm of k-th column without under/overflow.
 125              $norm = 0.0;
 126              for ($i = $k; $i < $this->rows; ++$i) {
 127                  $norm = $this->hypo($norm, $this->qrMatrix[$i][$k]);
 128              }
 129              if ($norm != 0.0) {
 130                  // Form k-th Householder vector.
 131                  if ($this->qrMatrix[$k][$k] < 0.0) {
 132                      $norm = -$norm;
 133                  }
 134                  for ($i = $k; $i < $this->rows; ++$i) {
 135                      $this->qrMatrix[$i][$k] /= $norm;
 136                  }
 137                  $this->qrMatrix[$k][$k] += 1.0;
 138                  // Apply transformation to remaining columns.
 139                  for ($j = $k + 1; $j < $this->columns; ++$j) {
 140                      $s = 0.0;
 141                      for ($i = $k; $i < $this->rows; ++$i) {
 142                          $s += $this->qrMatrix[$i][$k] * $this->qrMatrix[$i][$j];
 143                      }
 144                      $s = -$s / $this->qrMatrix[$k][$k];
 145                      for ($i = $k; $i < $this->rows; ++$i) {
 146                          $this->qrMatrix[$i][$j] += $s * $this->qrMatrix[$i][$k];
 147                      }
 148                  }
 149              }
 150              $this->rDiagonal[$k] = -$norm;
 151          }
 152      }
 153  
 154      public function isFullRank(): bool
 155      {
 156          for ($j = 0; $j < $this->columns; ++$j) {
 157              if ($this->rDiagonal[$j] == 0.0) {
 158                  return false;
 159              }
 160          }
 161  
 162          return true;
 163      }
 164  
 165      /**
 166       * Least squares solution of A*X = B.
 167       *
 168       * @param Matrix $B a Matrix with as many rows as A and any number of columns
 169       *
 170       * @throws Exception
 171       *
 172       * @return Matrix matrix that minimizes the two norm of Q*R*X-B
 173       */
 174      public function solve(Matrix $B): Matrix
 175      {
 176          if ($B->rows !== $this->rows) {
 177              throw new Exception('Matrix row dimensions are not equal');
 178          }
 179  
 180          if (!$this->isFullRank()) {
 181              throw new Exception('Can only perform this operation on a full-rank matrix');
 182          }
 183  
 184          // Compute Y = transpose(Q)*B
 185          $Y = $this->getQ()->transpose()
 186              ->multiply($B);
 187          // Solve R*X = Y;
 188          return $this->getR()->inverse()
 189              ->multiply($Y);
 190      }
 191  }