For God so loved the world, that He gave His only begotten Son, that all who believe in Him should not perish but have everlasting life
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

148 lines
4.4 KiB

  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. * Matrix row and column dimension
  28. * @var int
  29. * @access private
  30. */
  31. private $m;
  32. /**
  33. * Symmetric positive definite flag
  34. * @var boolean
  35. * @access private
  36. */
  37. private $isspd = true;
  38. /**
  39. * CholeskyDecomposition
  40. *
  41. * Class constructor - decomposes symmetric positive definite matrix
  42. * @param mixed Matrix square symmetric positive definite matrix
  43. */
  44. public function __construct($A = null)
  45. {
  46. if ($A instanceof Matrix) {
  47. $this->L = $A->getArray();
  48. $this->m = $A->getRowDimension();
  49. for ($i = 0; $i < $this->m; ++$i) {
  50. for ($j = $i; $j < $this->m; ++$j) {
  51. for ($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
  52. $sum -= $this->L[$i][$k] * $this->L[$j][$k];
  53. }
  54. if ($i == $j) {
  55. if ($sum >= 0) {
  56. $this->L[$i][$i] = sqrt($sum);
  57. } else {
  58. $this->isspd = false;
  59. }
  60. } else {
  61. if ($this->L[$i][$i] != 0) {
  62. $this->L[$j][$i] = $sum / $this->L[$i][$i];
  63. }
  64. }
  65. }
  66. for ($k = $i+1; $k < $this->m; ++$k) {
  67. $this->L[$i][$k] = 0.0;
  68. }
  69. }
  70. } else {
  71. throw new PHPExcel_Calculation_Exception(JAMAError(ARGUMENT_TYPE_EXCEPTION));
  72. }
  73. } // function __construct()
  74. /**
  75. * Is the matrix symmetric and positive definite?
  76. *
  77. * @return boolean
  78. */
  79. public function isSPD()
  80. {
  81. return $this->isspd;
  82. } // function isSPD()
  83. /**
  84. * getL
  85. *
  86. * Return triangular factor.
  87. * @return Matrix Lower triangular matrix
  88. */
  89. public function getL()
  90. {
  91. return new Matrix($this->L);
  92. } // function getL()
  93. /**
  94. * Solve A*X = B
  95. *
  96. * @param $B Row-equal matrix
  97. * @return Matrix L * L' * X = B
  98. */
  99. public function solve($B = null)
  100. {
  101. if ($B instanceof Matrix) {
  102. if ($B->getRowDimension() == $this->m) {
  103. if ($this->isspd) {
  104. $X = $B->getArrayCopy();
  105. $nx = $B->getColumnDimension();
  106. for ($k = 0; $k < $this->m; ++$k) {
  107. for ($i = $k + 1; $i < $this->m; ++$i) {
  108. for ($j = 0; $j < $nx; ++$j) {
  109. $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
  110. }
  111. }
  112. for ($j = 0; $j < $nx; ++$j) {
  113. $X[$k][$j] /= $this->L[$k][$k];
  114. }
  115. }
  116. for ($k = $this->m - 1; $k >= 0; --$k) {
  117. for ($j = 0; $j < $nx; ++$j) {
  118. $X[$k][$j] /= $this->L[$k][$k];
  119. }
  120. for ($i = 0; $i < $k; ++$i) {
  121. for ($j = 0; $j < $nx; ++$j) {
  122. $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
  123. }
  124. }
  125. }
  126. return new Matrix($X, $this->m, $nx);
  127. } else {
  128. throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));
  129. }
  130. } else {
  131. throw new PHPExcel_Calculation_Exception(JAMAError(MATRIX_DIMENSION_EXCEPTION));
  132. }
  133. } else {
  134. throw new PHPExcel_Calculation_Exception(JAMAError(ARGUMENT_TYPE_EXCEPTION));
  135. }
  136. } // function solve()
  137. }