CholeskyTest.php 6.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214
  1. <?php
  2. namespace MathPHP\Tests\LinearAlgebra\Decomposition;
  3. use MathPHP\LinearAlgebra\MatrixFactory;
  4. use MathPHP\Exception;
  5. class CholeskyTest extends \PHPUnit\Framework\TestCase
  6. {
  7. /**
  8. * @testCase choleskyDecomposition computes the expected lower triangular matrix
  9. * @dataProvider dataProviderForCholeskyDecomposition
  10. * @param array $A
  11. * @param array $expected_L
  12. * @throws \Exception
  13. */
  14. public function testCholeskyDecomposition(array $A, array $expected_L)
  15. {
  16. // Given
  17. $A = MatrixFactory::create($A);
  18. $expected_L = MatrixFactory::create($expected_L);
  19. $expected_Lᵀ = $expected_L->transpose();
  20. // When
  21. $cholesky = $A->choleskyDecomposition();
  22. $L = $cholesky->L;
  23. $Lᵀ = $cholesky->LT;
  24. // Then
  25. $this->assertEqualsWithDelta($expected_L->getMatrix(), $L->getMatrix(), 0.00001);
  26. $this->assertEqualsWithDelta($expected_Lᵀ->getMatrix(), $Lᵀ->getMatrix(), 0.00001);
  27. // And LLᵀ = A
  28. $LLᵀ = $L->multiply($Lᵀ);
  29. $this->assertEqualsWithDelta($A->getMatrix(), $LLᵀ->getMatrix(), 0.00001);
  30. }
  31. /**
  32. * @return array
  33. */
  34. public function dataProviderForCholeskyDecomposition(): array
  35. {
  36. return [
  37. // Example data from Wikipedia
  38. [
  39. [
  40. [4, 12, -16],
  41. [12, 37, -43],
  42. [-16, -43, 98],
  43. ],
  44. [
  45. [2, 0, 0],
  46. [6, 1, 0],
  47. [-8, 5, 3],
  48. ],
  49. ],
  50. // Example data from rosettacode.org
  51. [
  52. [
  53. [25, 15, -5],
  54. [15, 18, 0],
  55. [-5, 0, 11],
  56. ],
  57. [
  58. [5, 0, 0],
  59. [3, 3, 0],
  60. [-1, 1, 3],
  61. ],
  62. ],
  63. [
  64. [
  65. [18, 22, 54, 42],
  66. [22, 70, 86, 62],
  67. [54, 86, 174, 134],
  68. [42, 62, 134, 106],
  69. ],
  70. [
  71. [4.24264, 0.00000, 0.00000, 0.00000],
  72. [5.18545, 6.56591, 0.00000, 0.00000],
  73. [12.72792, 3.04604, 1.64974, 0.00000],
  74. [9.89949, 1.62455, 1.84971, 1.39262],
  75. ],
  76. ],
  77. // Example data from https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/04LinearAlgebra/cholesky/
  78. [
  79. [
  80. [5, 1.2, 0.3, -0.6],
  81. [1.2, 6, -0.4, 0.9],
  82. [0.3, -0.4, 8, 1.7],
  83. [-0.6, 0.9, 1.7, 10],
  84. ],
  85. [
  86. [2.236067977499790, 0, 0, 0],
  87. [0.536656314599949, 2.389979079406345, 0, 0],
  88. [0.134164078649987, -0.197491268466351, 2.818332343581848, 0],
  89. [-0.268328157299975, 0.436823907370487, 0.646577012719190, 3.052723872310221],
  90. ],
  91. ],
  92. [
  93. [
  94. [9.0000, 0.6000, -0.3000, 1.5000],
  95. [0.6000, 16.0400, 1.1800, -1.5000],
  96. [-0.3000, 1.1800, 4.1000, -0.5700],
  97. [1.5000, -1.5000, -0.5700, 25.4500],
  98. ],
  99. [
  100. [3, 0, 0, 0],
  101. [0.2, 4, 0, 0],
  102. [-0.1, 0.3, 2, 0],
  103. [0.5, -0.4, -0.2, 5],
  104. ],
  105. ],
  106. // Example data created with http://calculator.vhex.net/post/calculator-result/cholesky-decomposition
  107. [
  108. [
  109. [2, -1],
  110. [-1, 2],
  111. ],
  112. [
  113. [1.414214, 0],
  114. [-0.707107, 1.224745],
  115. ],
  116. ],
  117. [
  118. [
  119. [1, -1],
  120. [-1, 4],
  121. ],
  122. [
  123. [1, 0],
  124. [-1, 1.732051],
  125. ],
  126. ],
  127. [
  128. [
  129. [6, 4],
  130. [4, 5],
  131. ],
  132. [
  133. [2.44949, 0],
  134. [1.632993, 1.527525],
  135. ],
  136. ],
  137. [
  138. [
  139. [4, 1, -1],
  140. [1, 2, 1],
  141. [-1, 1, 2],
  142. ],
  143. [
  144. [2, 0, 0],
  145. [0.5, 1.322876, 0],
  146. [-0.5, 0.944911, 0.92582],
  147. ],
  148. ],
  149. [
  150. [
  151. [9, -3, 3, 9],
  152. [-3, 17, -1, -7],
  153. [3, -1, 17, 15],
  154. [9, -7, 15, 44],
  155. ],
  156. [
  157. [3, 0, 0, 0],
  158. [-1, 4, 0, 0],
  159. [1, 0, 4, 0],
  160. [3, -1, 3, 5],
  161. ],
  162. ],
  163. ];
  164. }
  165. /**
  166. * @test choleskyDecomposition throws a MatrixException if the matrix is not positive definite
  167. * @throws \Exception
  168. */
  169. public function testCholeskyDecompositionException()
  170. {
  171. // Given
  172. $A = [
  173. [1, 2, 3],
  174. [2, 3, 4],
  175. [3, 4, 5],
  176. ];
  177. $A = MatrixFactory::create($A);
  178. // Then
  179. $this->expectException(Exception\MatrixException::class);
  180. // When
  181. $L = $A->choleskyDecomposition();
  182. }
  183. /**
  184. * @test Cholesky Decomposition invalid property
  185. * @throws \Exception
  186. */
  187. public function testCholeskyDecompositionInvalidProperty()
  188. {
  189. // Given
  190. $A = MatrixFactory::create([
  191. [4, 1, -1],
  192. [1, 2, 1],
  193. [-1, 1, 2],
  194. ]);
  195. $cholesky = $A->choleskyDecomposition();
  196. // Then
  197. $this->expectException(Exception\MathException::class);
  198. // When
  199. $doesNotExist = $cholesky->doesNotExist;
  200. }
  201. }