[Summary view]

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

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