[Summary view]

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

   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($A)
  60      {
  61          if ($A instanceof Matrix) {
  62              // Initialize.
  63              $this->QR = $A->getArray();
  64              $this->m = $A->getRowDimension();
  65              $this->n = $A->getColumnDimension();
  66              // Main loop.
  67              for ($k = 0; $k < $this->n; ++$k) {
  68                  // Compute 2-norm of k-th column without under/overflow.
  69                  $nrm = 0.0;
  70                  for ($i = $k; $i < $this->m; ++$i) {
  71                      $nrm = hypo($nrm, $this->QR[$i][$k]);
  72                  }
  73                  if ($nrm != 0.0) {
  74                      // Form k-th Householder vector.
  75                      if ($this->QR[$k][$k] < 0) {
  76                          $nrm = -$nrm;
  77                      }
  78                      for ($i = $k; $i < $this->m; ++$i) {
  79                          $this->QR[$i][$k] /= $nrm;
  80                      }
  81                      $this->QR[$k][$k] += 1.0;
  82                      // Apply transformation to remaining columns.
  83                      for ($j = $k + 1; $j < $this->n; ++$j) {
  84                          $s = 0.0;
  85                          for ($i = $k; $i < $this->m; ++$i) {
  86                              $s += $this->QR[$i][$k] * $this->QR[$i][$j];
  87                          }
  88                          $s = -$s / $this->QR[$k][$k];
  89                          for ($i = $k; $i < $this->m; ++$i) {
  90                              $this->QR[$i][$j] += $s * $this->QR[$i][$k];
  91                          }
  92                      }
  93                  }
  94                  $this->Rdiag[$k] = -$nrm;
  95              }
  96          } else {
  97              throw new CalculationException(Matrix::ARGUMENT_TYPE_EXCEPTION);
  98          }
  99      }
 100  
 101      //    function __construct()
 102  
 103      /**
 104       *    Is the matrix full rank?
 105       *
 106       * @return bool true if R, and hence A, has full rank, else false
 107       */
 108      public function isFullRank()
 109      {
 110          for ($j = 0; $j < $this->n; ++$j) {
 111              if ($this->Rdiag[$j] == 0) {
 112                  return false;
 113              }
 114          }
 115  
 116          return true;
 117      }
 118  
 119      //    function isFullRank()
 120  
 121      /**
 122       * Return the Householder vectors.
 123       *
 124       * @return Matrix Lower trapezoidal matrix whose columns define the reflections
 125       */
 126      public function getH()
 127      {
 128          $H = [];
 129          for ($i = 0; $i < $this->m; ++$i) {
 130              for ($j = 0; $j < $this->n; ++$j) {
 131                  if ($i >= $j) {
 132                      $H[$i][$j] = $this->QR[$i][$j];
 133                  } else {
 134                      $H[$i][$j] = 0.0;
 135                  }
 136              }
 137          }
 138  
 139          return new Matrix($H);
 140      }
 141  
 142      //    function getH()
 143  
 144      /**
 145       * Return the upper triangular factor.
 146       *
 147       * @return Matrix upper triangular factor
 148       */
 149      public function getR()
 150      {
 151          $R = [];
 152          for ($i = 0; $i < $this->n; ++$i) {
 153              for ($j = 0; $j < $this->n; ++$j) {
 154                  if ($i < $j) {
 155                      $R[$i][$j] = $this->QR[$i][$j];
 156                  } elseif ($i == $j) {
 157                      $R[$i][$j] = $this->Rdiag[$i];
 158                  } else {
 159                      $R[$i][$j] = 0.0;
 160                  }
 161              }
 162          }
 163  
 164          return new Matrix($R);
 165      }
 166  
 167      //    function getR()
 168  
 169      /**
 170       * Generate and return the (economy-sized) orthogonal factor.
 171       *
 172       * @return Matrix orthogonal factor
 173       */
 174      public function getQ()
 175      {
 176          $Q = [];
 177          for ($k = $this->n - 1; $k >= 0; --$k) {
 178              for ($i = 0; $i < $this->m; ++$i) {
 179                  $Q[$i][$k] = 0.0;
 180              }
 181              $Q[$k][$k] = 1.0;
 182              for ($j = $k; $j < $this->n; ++$j) {
 183                  if ($this->QR[$k][$k] != 0) {
 184                      $s = 0.0;
 185                      for ($i = $k; $i < $this->m; ++$i) {
 186                          $s += $this->QR[$i][$k] * $Q[$i][$j];
 187                      }
 188                      $s = -$s / $this->QR[$k][$k];
 189                      for ($i = $k; $i < $this->m; ++$i) {
 190                          $Q[$i][$j] += $s * $this->QR[$i][$k];
 191                      }
 192                  }
 193              }
 194          }
 195  
 196          return new Matrix($Q);
 197      }
 198  
 199      //    function getQ()
 200  
 201      /**
 202       * Least squares solution of A*X = B.
 203       *
 204       * @param Matrix $B a Matrix with as many rows as A and any number of columns
 205       *
 206       * @return Matrix matrix that minimizes the two norm of Q*R*X-B
 207       */
 208      public function solve($B)
 209      {
 210          if ($B->getRowDimension() == $this->m) {
 211              if ($this->isFullRank()) {
 212                  // Copy right hand side
 213                  $nx = $B->getColumnDimension();
 214                  $X = $B->getArrayCopy();
 215                  // Compute Y = transpose(Q)*B
 216                  for ($k = 0; $k < $this->n; ++$k) {
 217                      for ($j = 0; $j < $nx; ++$j) {
 218                          $s = 0.0;
 219                          for ($i = $k; $i < $this->m; ++$i) {
 220                              $s += $this->QR[$i][$k] * $X[$i][$j];
 221                          }
 222                          $s = -$s / $this->QR[$k][$k];
 223                          for ($i = $k; $i < $this->m; ++$i) {
 224                              $X[$i][$j] += $s * $this->QR[$i][$k];
 225                          }
 226                      }
 227                  }
 228                  // Solve R*X = Y;
 229                  for ($k = $this->n - 1; $k >= 0; --$k) {
 230                      for ($j = 0; $j < $nx; ++$j) {
 231                          $X[$k][$j] /= $this->Rdiag[$k];
 232                      }
 233                      for ($i = 0; $i < $k; ++$i) {
 234                          for ($j = 0; $j < $nx; ++$j) {
 235                              $X[$i][$j] -= $X[$k][$j] * $this->QR[$i][$k];
 236                          }
 237                      }
 238                  }
 239                  $X = new Matrix($X);
 240  
 241                  return $X->getMatrix(0, $this->n - 1, 0, $nx);
 242              }
 243  
 244              throw new CalculationException(self::MATRIX_RANK_EXCEPTION);
 245          }
 246  
 247          throw new CalculationException(Matrix::MATRIX_DIMENSION_EXCEPTION);
 248      }
 249  }