<?php /** * @package JAMA * * Cholesky decomposition class * * For a symmetric, positive definite matrix A, the Cholesky decomposition * is an lower triangular matrix L so that A = L*L'. * * If the matrix is not symmetric or positive definite, the constructor * returns a partial decomposition and sets an internal flag that may * be queried by the isSPD() method. * * @author Paul Meagher * @author Michael Bommarito * @version 1.2 */ class CholeskyDecomposition { /** * Decomposition storage * @var array * @access private */ private $L = array(); /** * Matrix row and column dimension * @var int * @access private */ private $m; /** * Symmetric positive definite flag * @var boolean * @access private */ private $isspd = true; /** * CholeskyDecomposition * * Class constructor - decomposes symmetric positive definite matrix * @param mixed Matrix square symmetric positive definite matrix */ public function __construct($A = null) { if ($A instanceof Matrix) { $this->L = $A->getArray(); $this->m = $A->getRowDimension(); for($i = 0; $i < $this->m; ++$i) { for($j = $i; $j < $this->m; ++$j) { for($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) { $sum -= $this->L[$i][$k] * $this->L[$j][$k]; } if ($i == $j) { if ($sum >= 0) { $this->L[$i][$i] = sqrt($sum); } else { $this->isspd = false; } } else { if ($this->L[$i][$i] != 0) { $this->L[$j][$i] = $sum / $this->L[$i][$i]; } } } for ($k = $i+1; $k < $this->m; ++$k) { $this->L[$i][$k] = 0.0; } } } else { throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException)); } } // function __construct() /** * Is the matrix symmetric and positive definite? * * @return boolean */ public function isSPD() { return $this->isspd; } // function isSPD() /** * getL * * Return triangular factor. * @return Matrix Lower triangular matrix */ public function getL() { return new Matrix($this->L); } // function getL() /** * Solve A*X = B * * @param $B Row-equal matrix * @return Matrix L * L' * X = B */ public function solve($B = null) { if ($B instanceof Matrix) { if ($B->getRowDimension() == $this->m) { if ($this->isspd) { $X = $B->getArrayCopy(); $nx = $B->getColumnDimension(); for ($k = 0; $k < $this->m; ++$k) { for ($i = $k + 1; $i < $this->m; ++$i) { for ($j = 0; $j < $nx; ++$j) { $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k]; } } for ($j = 0; $j < $nx; ++$j) { $X[$k][$j] /= $this->L[$k][$k]; } } for ($k = $this->m - 1; $k >= 0; --$k) { for ($j = 0; $j < $nx; ++$j) { $X[$k][$j] /= $this->L[$k][$k]; } for ($i = 0; $i < $k; ++$i) { for ($j = 0; $j < $nx; ++$j) { $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i]; } } } return new Matrix($X, $this->m, $nx); } else { throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException)); } } else { throw new PHPExcel_Calculation_Exception(JAMAError(MatrixDimensionException)); } } else { throw new PHPExcel_Calculation_Exception(JAMAError(ArgumentTypeException)); } } // function solve() } // class CholeskyDecomposition