CholeskyDecomposition.php 4.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137
  1. <?php
  2. /**
  3. * 重庆赤晓店信息科技有限公司
  4. * https://www.chixiaodian.com
  5. * Copyright (c) 2023 赤店商城 All rights reserved.
  6. */
  7. class CholeskyDecomposition
  8. {
  9. /**
  10. * Decomposition storage
  11. * @var array
  12. * @access private
  13. */
  14. private $L = array();
  15. /**
  16. * Matrix row and column dimension
  17. * @var int
  18. * @access private
  19. */
  20. private $m;
  21. /**
  22. * Symmetric positive definite flag
  23. * @var boolean
  24. * @access private
  25. */
  26. private $isspd = true;
  27. /**
  28. * CholeskyDecomposition
  29. *
  30. * Class constructor - decomposes symmetric positive definite matrix
  31. * @param mixed Matrix square symmetric positive definite matrix
  32. */
  33. public function __construct($A = null)
  34. {
  35. if ($A instanceof Matrix) {
  36. $this->L = $A->getArray();
  37. $this->m = $A->getRowDimension();
  38. for ($i = 0; $i < $this->m; ++$i) {
  39. for ($j = $i; $j < $this->m; ++$j) {
  40. for ($sum = $this->L[$i][$j], $k = $i - 1; $k >= 0; --$k) {
  41. $sum -= $this->L[$i][$k] * $this->L[$j][$k];
  42. }
  43. if ($i == $j) {
  44. if ($sum >= 0) {
  45. $this->L[$i][$i] = sqrt($sum);
  46. } else {
  47. $this->isspd = false;
  48. }
  49. } else {
  50. if ($this->L[$i][$i] != 0) {
  51. $this->L[$j][$i] = $sum / $this->L[$i][$i];
  52. }
  53. }
  54. }
  55. for ($k = $i+1; $k < $this->m; ++$k) {
  56. $this->L[$i][$k] = 0.0;
  57. }
  58. }
  59. } else {
  60. throw new PHPExcel_Calculation_Exception(JAMAError(ARGUMENT_TYPE_EXCEPTION));
  61. }
  62. } // function __construct()
  63. /**
  64. * Is the matrix symmetric and positive definite?
  65. *
  66. * @return boolean
  67. */
  68. public function isSPD()
  69. {
  70. return $this->isspd;
  71. } // function isSPD()
  72. /**
  73. * getL
  74. *
  75. * Return triangular factor.
  76. * @return Matrix Lower triangular matrix
  77. */
  78. public function getL()
  79. {
  80. return new Matrix($this->L);
  81. } // function getL()
  82. /**
  83. * Solve A*X = B
  84. *
  85. * @param $B Row-equal matrix
  86. * @return Matrix L * L' * X = B
  87. */
  88. public function solve($B = null)
  89. {
  90. if ($B instanceof Matrix) {
  91. if ($B->getRowDimension() == $this->m) {
  92. if ($this->isspd) {
  93. $X = $B->getArrayCopy();
  94. $nx = $B->getColumnDimension();
  95. for ($k = 0; $k < $this->m; ++$k) {
  96. for ($i = $k + 1; $i < $this->m; ++$i) {
  97. for ($j = 0; $j < $nx; ++$j) {
  98. $X[$i][$j] -= $X[$k][$j] * $this->L[$i][$k];
  99. }
  100. }
  101. for ($j = 0; $j < $nx; ++$j) {
  102. $X[$k][$j] /= $this->L[$k][$k];
  103. }
  104. }
  105. for ($k = $this->m - 1; $k >= 0; --$k) {
  106. for ($j = 0; $j < $nx; ++$j) {
  107. $X[$k][$j] /= $this->L[$k][$k];
  108. }
  109. for ($i = 0; $i < $k; ++$i) {
  110. for ($j = 0; $j < $nx; ++$j) {
  111. $X[$i][$j] -= $X[$k][$j] * $this->L[$k][$i];
  112. }
  113. }
  114. }
  115. return new Matrix($X, $this->m, $nx);
  116. } else {
  117. throw new PHPExcel_Calculation_Exception(JAMAError(MatrixSPDException));
  118. }
  119. } else {
  120. throw new PHPExcel_Calculation_Exception(JAMAError(MATRIX_DIMENSION_EXCEPTION));
  121. }
  122. } else {
  123. throw new PHPExcel_Calculation_Exception(JAMAError(ARGUMENT_TYPE_EXCEPTION));
  124. }
  125. } // function solve()
  126. }