[Summary view] [Print] [Text view]

   1  <?php
   2  /**
   3   *    @package JAMA
   4   *
   5   *    Cholesky decomposition class
   6   *
   7   *    For a symmetric, positive definite matrix A, the Cholesky decomposition
   8   *    is an lower triangular matrix L so that A = L*L'.
   9   *
  10   *    If the matrix is not symmetric or positive definite, the constructor
  11   *    returns a partial decomposition and sets an internal flag that may
  12   *    be queried by the isSPD() method.
  13   *
  14   *    @author Paul Meagher
  15   *    @author Michael Bommarito
  16   *    @version 1.2
  17   */
  18  class CholeskyDecomposition
  19  {
  20      /**
  21       *    Decomposition storage
  22       *    @var array
  23       *    @access private
  24       */
  25      private $L = array();
  26  
  27      /**
  28       *    Matrix row and column dimension
  29       *    @var int
  30       *    @access private
  31       */
  32      private $m;
  33  
  34      /**
  35       *    Symmetric positive definite flag
  36       *    @var boolean
  37       *    @access private
  38       */
  39      private $isspd = true;
  40  
  41      /**
  42       *    CholeskyDecomposition
  43       *
  44       *    Class constructor - decomposes symmetric positive definite matrix
  45       *    @param mixed Matrix square symmetric positive definite matrix
  46       */
  47      public function __construct($A = null)
  48      {
  49          if ($A instanceof Matrix) {
  50              $this->L = $A->getArray();
  51              $this->m = $A->getRowDimension();
  52  
  53              for ($i = 0; $i < $this->m; ++$i) {
  54                  for ($j = $i; $j < $this->m; ++$j) {
  55                      for ($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
  56                          $sum -= $this->L[$i][$k] * $this->L[$j][$k];
  57                      }
  58                      if ($i == $j) {
  59                          if ($sum >= 0) {
  60                              $this->L[$i][$i] = sqrt($sum);
  61                          } else {
  62                              $this->isspd = false;
  63                          }
  64                      } else {
  65                          if ($this->L[$i][$i] != 0) {
  66                              $this->L[$j][$i] = $sum / $this->L[$i][$i];
  67                          }
  68                      }
  69                  }
  70  
  71                  for ($k = $i+1; $k < $this->m; ++$k) {
  72                      $this->L[$i][$k] = 0.0;
  73                  }
  74              }
  75          } else {
  76              throw new PHPExcel_Calculation_Exception(JAMAError(ARGUMENT_TYPE_EXCEPTION));
  77          }
  78      }    //    function __construct()
  79  
  80      /**
  81       *    Is the matrix symmetric and positive definite?
  82       *
  83       *    @return boolean
  84       */
  85      public function isSPD()
  86      {
  87          return $this->isspd;
  88      }    //    function isSPD()
  89  
  90      /**
  91       *    getL
  92       *
  93       *    Return triangular factor.
  94       *    @return Matrix Lower triangular matrix
  95       */
  96      public function getL()
  97      {
  98          return new Matrix($this->L);
  99      }    //    function getL()
 100  
 101      /**
 102       *    Solve A*X = B
 103       *
 104       *    @param $B Row-equal matrix
 105       *    @return Matrix L * L' * X = B
 106       */
 107      public function solve($B = null)
 108      {
 109          if ($B instanceof Matrix) {
 110              if ($B->getRowDimension() == $this->m) {
 111                  if ($this->isspd) {
 112                      $X  = $B->getArrayCopy();
 113                      $nx = $B->getColumnDimension();
 114  
 115                      for ($k = 0; $k < $this->m; ++$k) {
 116                          for ($i = $k + 1; $i < $this->m; ++$i) {
 117                              for ($j = 0; $j < $nx; ++$j) {
 118                                  $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
 119                              }
 120                          }
 121                          for ($j = 0; $j < $nx; ++$j) {
 122                              $X[$k][$j] /= $this->L[$k][$k];
 123                          }
 124                      }
 125  
 126                      for ($k = $this->m - 1; $k >= 0; --$k) {
 127                          for ($j = 0; $j < $nx; ++$j) {
 128                              $X[$k][$j] /= $this->L[$k][$k];
 129                          }
 130                          for ($i = 0; $i < $k; ++$i) {
 131                              for ($j = 0; $j < $nx; ++$j) {
 132                                  $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
 133                              }
 134                          }
 135                      }
 136  
 137                      return new Matrix($X, $this->m, $nx);
 138                  } else {
 139                      throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));
 140                  }
 141              } else {
 142                  throw new PHPExcel_Calculation_Exception(JAMAError(MATRIX_DIMENSION_EXCEPTION));
 143              }
 144          } else {
 145              throw new PHPExcel_Calculation_Exception(JAMAError(ARGUMENT_TYPE_EXCEPTION));
 146          }
 147      }    //    function solve()
 148  }