[Summary view]

Differences Between: [Versions 310 and 401] [Versions 311 and 401] [Versions 39 and 401]

   1  <?php
   2  
   3  namespace PhpOffice\PhpSpreadsheet\Shared\JAMA;
   4  
   5  use PhpOffice\PhpSpreadsheet\Calculation\Exception as CalculationException;
   6  
   7  /**
   8   *    For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
   9   *    orthogonal matrix Q and an n-by-n upper triangular matrix R so that
  10   *    A = Q*R.
  11   *
  12   *    The QR decompostion always exists, even if the matrix does not have
  13   *    full rank, so the constructor will never fail.  The primary use of the
  14   *    QR decomposition is in the least squares solution of nonsquare systems
  15   *    of simultaneous linear equations.  This will fail if isFullRank()
  16   *    returns false.
  17   *
  18   * @author  Paul Meagher
  19   *
  20   * @version 1.1
  21   */
  22  class QRDecomposition
  23  {
  24      const MATRIX_RANK_EXCEPTION = 'Can only perform operation on full-rank matrix.';
  25  
  26      /**
  27       * Array for internal storage of decomposition.
  28       *
  29       * @var array
  30       */
  31      private $QR = [];
  32  
  33      /**
  34       * Row dimension.
  35       *
  36       * @var int
  37       */
  38      private $m;
  39  
  40      /**
  41       * Column dimension.
  42       *
  43       * @var int
  44       */
  45      private $n;
  46  
  47      /**
  48       * Array for internal storage of diagonal of R.
  49       *
  50       * @var array
  51       */
  52      private $Rdiag = [];
  53  
  54      /**
  55       * QR Decomposition computed by Householder reflections.
  56       *
  57       * @param Matrix $A Rectangular matrix
  58       */
  59      public function __construct(Matrix $A)
  60      {
  61          // Initialize.
  62          $this->QR = $A->getArray();
  63          $this->m = $A->getRowDimension();
  64          $this->n = $A->getColumnDimension();
  65          // Main loop.
  66          for ($k = 0; $k < $this->n; ++$k) {
  67              // Compute 2-norm of k-th column without under/overflow.
  68              $nrm = 0.0;
  69              for ($i = $k; $i < $this->m; ++$i) {
  70                  $nrm = hypo($nrm, $this->QR[$i][$k]);
  71              }
  72              if ($nrm != 0.0) {
  73                  // Form k-th Householder vector.
  74                  if ($this->QR[$k][$k] < 0) {
  75                      $nrm = -$nrm;
  76                  }
  77                  for ($i = $k; $i < $this->m; ++$i) {
  78                      $this->QR[$i][$k] /= $nrm;
  79                  }
  80                  $this->QR[$k][$k] += 1.0;
  81                  // Apply transformation to remaining columns.
  82                  for ($j = $k + 1; $j < $this->n; ++$j) {
  83                      $s = 0.0;
  84                      for ($i = $k; $i < $this->m; ++$i) {
  85                          $s += $this->QR[$i][$k] * $this->QR[$i][$j];
  86                      }
  87                      $s = -$s / $this->QR[$k][$k];
  88                      for ($i = $k; $i < $this->m; ++$i) {
  89                          $this->QR[$i][$j] += $s * $this->QR[$i][$k];
  90                      }
  91                  }
  92              }
  93              $this->Rdiag[$k] = -$nrm;
  94          }
  95      }
  96  
  97      //    function __construct()
  98  
  99      /**
 100       *    Is the matrix full rank?
 101       *
 102       * @return bool true if R, and hence A, has full rank, else false
 103       */
 104      public function isFullRank()
 105      {
 106          for ($j = 0; $j < $this->n; ++$j) {
 107              if ($this->Rdiag[$j] == 0) {
 108                  return false;
 109              }
 110          }
 111  
 112          return true;
 113      }
 114  
 115      //    function isFullRank()
 116  
 117      /**
 118       * Return the Householder vectors.
 119       *
 120       * @return Matrix Lower trapezoidal matrix whose columns define the reflections
 121       */
 122      public function getH()
 123      {
 124          $H = [];
 125          for ($i = 0; $i < $this->m; ++$i) {
 126              for ($j = 0; $j < $this->n; ++$j) {
 127                  if ($i >= $j) {
 128                      $H[$i][$j] = $this->QR[$i][$j];
 129                  } else {
 130                      $H[$i][$j] = 0.0;
 131                  }
 132              }
 133          }
 134  
 135          return new Matrix($H);
 136      }
 137  
 138      //    function getH()
 139  
 140      /**
 141       * Return the upper triangular factor.
 142       *
 143       * @return Matrix upper triangular factor
 144       */
 145      public function getR()
 146      {
 147          $R = [];
 148          for ($i = 0; $i < $this->n; ++$i) {
 149              for ($j = 0; $j < $this->n; ++$j) {
 150                  if ($i < $j) {
 151                      $R[$i][$j] = $this->QR[$i][$j];
 152                  } elseif ($i == $j) {
 153                      $R[$i][$j] = $this->Rdiag[$i];
 154                  } else {
 155                      $R[$i][$j] = 0.0;
 156                  }
 157              }
 158          }
 159  
 160          return new Matrix($R);
 161      }
 162  
 163      //    function getR()
 164  
 165      /**
 166       * Generate and return the (economy-sized) orthogonal factor.
 167       *
 168       * @return Matrix orthogonal factor
 169       */
 170      public function getQ()
 171      {
 172          $Q = [];
 173          for ($k = $this->n - 1; $k >= 0; --$k) {
 174              for ($i = 0; $i < $this->m; ++$i) {
 175                  $Q[$i][$k] = 0.0;
 176              }
 177              $Q[$k][$k] = 1.0;
 178              for ($j = $k; $j < $this->n; ++$j) {
 179                  if ($this->QR[$k][$k] != 0) {
 180                      $s = 0.0;
 181                      for ($i = $k; $i < $this->m; ++$i) {
 182                          $s += $this->QR[$i][$k] * $Q[$i][$j];
 183                      }
 184                      $s = -$s / $this->QR[$k][$k];
 185                      for ($i = $k; $i < $this->m; ++$i) {
 186                          $Q[$i][$j] += $s * $this->QR[$i][$k];
 187                      }
 188                  }
 189              }
 190          }
 191  
 192          return new Matrix($Q);
 193      }
 194  
 195      //    function getQ()
 196  
 197      /**
 198       * Least squares solution of A*X = B.
 199       *
 200       * @param Matrix $B a Matrix with as many rows as A and any number of columns
 201       *
 202       * @return Matrix matrix that minimizes the two norm of Q*R*X-B
 203       */
 204      public function solve(Matrix $B)
 205      {
 206          if ($B->getRowDimension() == $this->m) {
 207              if ($this->isFullRank()) {
 208                  // Copy right hand side
 209                  $nx = $B->getColumnDimension();
 210                  $X = $B->getArray();
 211                  // Compute Y = transpose(Q)*B
 212                  for ($k = 0; $k < $this->n; ++$k) {
 213                      for ($j = 0; $j < $nx; ++$j) {
 214                          $s = 0.0;
 215                          for ($i = $k; $i < $this->m; ++$i) {
 216                              $s += $this->QR[$i][$k] * $X[$i][$j];
 217                          }
 218                          $s = -$s / $this->QR[$k][$k];
 219                          for ($i = $k; $i < $this->m; ++$i) {
 220                              $X[$i][$j] += $s * $this->QR[$i][$k];
 221                          }
 222                      }
 223                  }
 224                  // Solve R*X = Y;
 225                  for ($k = $this->n - 1; $k >= 0; --$k) {
 226                      for ($j = 0; $j < $nx; ++$j) {
 227                          $X[$k][$j] /= $this->Rdiag[$k];
 228                      }
 229                      for ($i = 0; $i < $k; ++$i) {
 230                          for ($j = 0; $j < $nx; ++$j) {
 231                              $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
 232                          }
 233                      }
 234                  }
 235                  $X = new Matrix($X);
 236  
 237                  return $X->getMatrix(0, $this->n - 1, 0, $nx);
 238              }
 239  
 240              throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
 241          }
 242  
 243          throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
 244      }
 245  }