vendor/phpoffice/phpexcel/Classes/PHPExcel/Shared/JAMA/QRDecomposition.php

   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