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