[Summary view] [Print] [Text view]

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