QRTest.php 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534
  1. <?php
  2. namespace MathPHP\Tests\LinearAlgebra\Decomposition;
  3. use MathPHP\LinearAlgebra\MatrixFactory;
  4. use MathPHP\Exception;
  5. use MathPHP\LinearAlgebra\Vector;
  6. use MathPHP\Tests\LinearAlgebra\Fixture\MatrixDataProvider;
  7. class QRTest extends \PHPUnit\Framework\TestCase
  8. {
  9. use MatrixDataProvider;
  10. /**
  11. * @test qrDecomposition property A = QR
  12. * @dataProvider dataProviderForQrDecompositionSquareMatricesWithSpecificResults
  13. * @dataProvider dataProviderForQrDecompositionNonSquareMatricesWithSpecificResults
  14. * @dataProvider dataProviderForSingularMatrix
  15. * @dataProvider dataProviderForNonsingularMatrix
  16. * @param array $A
  17. * @throws \Exception
  18. */
  19. public function testQrDecompositionPropertyAEqualsQR(array $A)
  20. {
  21. // Given
  22. $A = MatrixFactory::create($A);
  23. // And
  24. $qrDecomposition = $A->qrDecomposition();
  25. // When
  26. $QR = $qrDecomposition->Q->multiply($qrDecomposition->R);
  27. // Then A = QR
  28. $this->assertEqualsWithDelta($A->getMatrix(), $QR->getMatrix(), 0.00001);
  29. }
  30. /**
  31. * @test qrDecomposition properties Q is orthogonal and R is upper triangular
  32. * @dataProvider dataProviderForQrDecompositionSquareMatricesWithSpecificResults
  33. * @dataProvider dataProviderForSingularMatrix
  34. * @dataProvider dataProviderForNonsingularMatrix
  35. * @param array $A
  36. * @throws \Exception
  37. */
  38. public function testQrDecompositionPropertiesOfQR(array $A)
  39. {
  40. // Given
  41. $A = MatrixFactory::create($A);
  42. // When
  43. $qr = $A->qrDecomposition();
  44. // Then Q is orthogonal and R is upper triangular
  45. $this->assertTrue($qr->Q->isOrthogonal());
  46. $this->assertTrue($qr->R->isUpperTriangular());
  47. }
  48. /**
  49. * @test qrDecomposition returns the expected array of Q and R factorized matrices
  50. * This test could be removed. It is testing a specific implementation of QR decomposition.
  51. * @dataProvider dataProviderForQrDecompositionSquareMatricesWithSpecificResults
  52. * @dataProvider dataProviderForQrDecompositionNonSquareMatricesWithSpecificResults
  53. * @param array $A
  54. * @param array $Q
  55. * @param array $R
  56. * @throws \Exception
  57. */
  58. public function testQrDecompositionResultMatrices(array $A, array $Q, array $R)
  59. {
  60. // Given
  61. $A = MatrixFactory::create($A);
  62. $Q = MatrixFactory::create($Q);
  63. $R = MatrixFactory::create($R);
  64. // When
  65. $qr = $A->qrDecomposition();
  66. $qrQ = $qr->Q;
  67. $qrR = $qr->R;
  68. // And Q and R are expected solution to QR decomposition
  69. $this->assertEqualsWithDelta($R->getMatrix(), $qrR->getMatrix(), 0.00001);
  70. $this->assertEqualsWithDelta($Q->getMatrix(), $qrQ->getMatrix(), 0.00001);
  71. }
  72. /**
  73. * @test Orthonormal matrix Q has the property QᵀQ = I
  74. * @dataProvider dataProviderForQrDecompositionSquareMatricesWithSpecificResults
  75. * @dataProvider dataProviderForQrDecompositionNonSquareMatricesWithSpecificResults
  76. * @dataProvider dataProviderForSingularMatrix
  77. * @dataProvider dataProviderForNonsingularMatrix
  78. * @param array $A
  79. * @throws \Exception
  80. */
  81. public function testQrDecompositionOrthonormalMatrixQPropertyQTransposeQIsIdentity(array $A)
  82. {
  83. // Given
  84. $A = MatrixFactory::create($A);
  85. $I = MatrixFactory::identity(min($A->getM(), $A->getN()));
  86. // And
  87. $qr = $A->qrDecomposition();
  88. // When
  89. $QᵀQ = $qr->Q->transpose()->multiply($qr->Q);
  90. // Then QᵀQ = I
  91. $this->assertEqualsWithDelta($I->getMatrix(), $QᵀQ->getMatrix(), 0.000001);
  92. }
  93. /**
  94. * @test qrDecomposition property R = QᵀA
  95. * @dataProvider dataProviderForQrDecompositionSquareMatricesWithSpecificResults
  96. * @dataProvider dataProviderForQrDecompositionNonSquareMatricesWithSpecificResults
  97. * @dataProvider dataProviderForSingularMatrix
  98. * @dataProvider dataProviderForNonsingularMatrix
  99. * @param array $A
  100. * @throws \Exception
  101. */
  102. public function testQrDecompositionPropertyREqualsQTransposeA(array $A)
  103. {
  104. // Given
  105. $A = MatrixFactory::create($A);
  106. // And
  107. $qrDecomposition = $A->qrDecomposition();
  108. // When
  109. $QᵀA = $qrDecomposition->Q->transpose()->multiply($A);
  110. // Then R = QᵀA
  111. $this->assertEqualsWithDelta($qrDecomposition->R->getMatrix(), $QᵀA->getMatrix(), 0.00001);
  112. }
  113. /**
  114. * @test qrDecomposition property Qᵀ = Q⁻¹
  115. * @dataProvider dataProviderForQrDecompositionSquareMatricesWithSpecificResults
  116. * @dataProvider dataProviderForSingularMatrix
  117. * @dataProvider dataProviderForNonsingularMatrix
  118. * @param array $A
  119. * @throws \Exception
  120. */
  121. public function testQrDecompositionPropertyQTransposeEqualsQInverse(array $A)
  122. {
  123. // Given
  124. $A = MatrixFactory::create($A);
  125. // And
  126. $Q = $A->qrDecomposition()->Q;
  127. // When
  128. $Qᵀ = $Q->transpose();
  129. $Q⁻¹ = $Q->inverse();
  130. // Then Qᵀ = Q⁻¹
  131. $this->assertEqualsWithDelta($Qᵀ->getMatrix(), $Q⁻¹->getMatrix(), 0.00001);
  132. }
  133. /**
  134. * @test Solve
  135. * @dataProvider dataProviderForSolve
  136. * @param array $A
  137. * @param array $b
  138. * @param array $expected
  139. * @throws \Exception
  140. */
  141. public function testSolve(array $A, array $b, array $expected)
  142. {
  143. // Given
  144. $A = MatrixFactory::create($A);
  145. $QR = $A->qrDecomposition();
  146. // And
  147. $expected = new Vector($expected);
  148. // When
  149. $x = $QR->solve($b);
  150. // Then
  151. $this->assertEqualsWithDelta($expected, $x, 0.00001);
  152. }
  153. /**
  154. * @return array (A, Q, R)
  155. */
  156. public function dataProviderForQrDecompositionSquareMatricesWithSpecificResults(): array
  157. {
  158. return [
  159. [
  160. 'A' => [
  161. [0],
  162. ],
  163. 'Q' => [
  164. [1],
  165. ],
  166. 'R' => [
  167. [0],
  168. ],
  169. ],
  170. [
  171. 'A' => [
  172. [1],
  173. ],
  174. 'Q' => [
  175. [1],
  176. ],
  177. 'R' => [
  178. [1],
  179. ],
  180. ],
  181. [
  182. 'A' => [
  183. [0, 0],
  184. [0, 0],
  185. ],
  186. 'Q' => [
  187. [1, 0],
  188. [0, 1]
  189. ],
  190. 'R' => [
  191. [0, 0],
  192. [0, 0],
  193. ],
  194. ],
  195. [
  196. 'A' => [
  197. [1, 1],
  198. [1, 1],
  199. ],
  200. 'Q' => [
  201. [-0.7071068, -0.7071068],
  202. [-0.7071068, 0.7071068]
  203. ],
  204. 'R' => [
  205. [-1.414214, -1.414214],
  206. [0, 0],
  207. ],
  208. ],
  209. [
  210. 'A' => [
  211. [6, 3],
  212. [8, 4],
  213. ],
  214. 'Q' => [
  215. [-0.6, -0.8],
  216. [-0.8, 0.6]
  217. ],
  218. 'R' => [
  219. [-10, -5],
  220. [0, 0],
  221. ],
  222. ],
  223. [
  224. 'A' => [
  225. [2, -2, 18],
  226. [2, 1, 0],
  227. [1, 2, 0],
  228. ],
  229. 'Q' => [
  230. [-2 / 3, 2 / 3, 1 / 3],
  231. [-2 / 3, -1 / 3, -2 / 3],
  232. [-1 / 3, -2 / 3, 2 / 3],
  233. ],
  234. 'R' => [
  235. [-3, 0, -12],
  236. [ 0, -3, 12],
  237. [ 0, 0, 6],
  238. ],
  239. ],
  240. [
  241. 'A' => [
  242. [12, -51, 4],
  243. [ 6, 167, -68],
  244. [-4, 24, -41],
  245. ],
  246. 'Q' => [
  247. [ -0.85714286, 0.39428571, 0.33142857],
  248. [ -0.42857143, -0.90285714, -0.03428571],
  249. [0.28571429, -0.17142857, 0.94285714],
  250. ],
  251. 'R' => [
  252. [-14, -21, 14],
  253. [ 0, -175, 70],
  254. [ 0, 0, -35],
  255. ],
  256. ],
  257. [
  258. 'A' => [
  259. [4, 3, 7],
  260. [1, 3, 6],
  261. [8, 5, 7],
  262. ],
  263. 'Q' => [
  264. [-0.4444444, -0.1194133, -0.8878117],
  265. [-0.1111111, -0.9760737, 0.1869077],
  266. [-0.8888889, 0.1817158, 0.4205424],
  267. ],
  268. 'R' => [
  269. [-9, -6.111111, -10.000000],
  270. [0, -2.377882, -5.420324],
  271. [0, 0.000000, -2.149439],
  272. ],
  273. ],
  274. [
  275. 'A' => [
  276. [1, 2, 3, 2],
  277. [4, 5, 6, 2],
  278. [7, 8, 9, 2],
  279. [4, 5, 5, 6],
  280. ],
  281. 'Q' => [
  282. [-0.1104315, 0.8589557, 0.2886751, -4.082483e-01],
  283. [-0.4417261, 0.2342606, 0.2886751, 8.164966e-01],
  284. [-0.7730207, -0.3904344, 0.2886751, -4.082483e-01],
  285. [-0.4417261, 0.2342606, -0.8660254, -9.159340e-16],
  286. ],
  287. 'R' => [
  288. [-9.055385, -10.8222896, -12.1474679, -5.300713e+00],
  289. [ 0.000000, 0.9370426, 1.6398245, 2.811128e+00],
  290. [ 0.000000, 0.0000000, 0.8660254, -3.464102e+00],
  291. [ 0.000000, 0.0000000, 0.0000000, -5.329071e-15],
  292. ],
  293. ],
  294. [
  295. 'A' => [
  296. [7, 8, 9, 2],
  297. [1, 2, 3, 2],
  298. [4, -3, 2, 12],
  299. [4, 1, -6, 6],
  300. ],
  301. 'Q' => [
  302. [-0.7730207, -0.5413835, -0.2274010, -0.24006603],
  303. [-0.1104315, -0.2016919, -0.1581711, 0.96026411],
  304. [-0.4417261, 0.7890753, -0.4245139, 0.04501238],
  305. [-0.4417261, 0.2087688, 0.8620085, 0.13503714],
  306. ],
  307. 'R' => [
  308. [-9.055385, -5.521576, -5.521576, -9.7179743],
  309. [ 0.000000, -6.892909, -5.151990, 9.2353658],
  310. [ 0.000000, 0.000000, -8.542201, -0.6932606],
  311. [ 0.000000, 0.000000, 0.000000, 2.7907676],
  312. ],
  313. ],
  314. [
  315. 'A' => [
  316. [3, 7, 6, 4, 5],
  317. [2, 3, 6, 5, 8],
  318. [2, 3, 4, 1, 0],
  319. [3, 7, 6, 7, 7],
  320. [1, 3, 4, 9, 4],
  321. ],
  322. 'Q' => [
  323. [-0.5773503, 0.3086067, -0.2229113, 0.7003493, 0.1767767],
  324. [-0.3849002, -0.5657789, 0.5387023, 0.2150195, -0.4419417],
  325. [-0.3849002, -0.5657789, -0.2414872, -0.3010273, 0.6187184],
  326. [-0.5773503, 0.3086067, -0.2229113, -0.5713376, -0.4419417],
  327. [-0.1924501, 0.4114756, 0.7430376, -0.2150195, 0.4419417],
  328. ],
  329. 'R' => [
  330. [-5.196152, -10.969655, -11.5470054, -10.392305, -10.7772050],
  331. [ 0.000000, 2.160247, -0.3086067, 3.703280, 0.8229512],
  332. [ 0.000000, 0.000000, 2.5634798, 6.687339, 4.6068332],
  333. [ 0.000000, 0.000000, 0.0000000, -2.359071, 0.3624615],
  334. [ 0.000000, 0.000000, 0.0000000, 0.000000, -3.9774756],
  335. ],
  336. ],
  337. [
  338. 'A' => [
  339. [1, 0, 0, 0, 0],
  340. [0, 0, 1, 0, 0],
  341. [1, -7, 0, 4, 2],
  342. [0, 4, 2, -7, 1],
  343. [0, 2, 0, 1, -7],
  344. ],
  345. 'Q' => [
  346. [-0.7071068, -0.5246722, 0.3333983, -0.01889766, -0.3364633],
  347. [ 0.0000000, 0.0000000, -0.5298652, 0.63622126, -0.5607722],
  348. [-0.7071068, 0.5246722, -0.3333983, 0.01889766, 0.3364633],
  349. [ 0.0000000, -0.5996254, -0.6787037, -0.31811063, 0.2803861],
  350. [ 0.0000000, -0.2998127, 0.1905133, 0.70236308, 0.6168494],
  351. ],
  352. 'R' => [
  353. [-1.414214, 4.949747, 0.000000, -2.828427, -1.414214],
  354. [ 0.000000, -6.670832, -1.199251, 5.996254, 2.548408],
  355. [ 0.000000, 0.000000, -1.887273, 3.607846, -2.679094],
  356. [ 0.000000, 0.000000, 0.000000, 3.004728, -5.196857],
  357. [ 0.000000, 0.000000, 0.000000, 0.000000, -3.364633],
  358. ],
  359. ],
  360. ];
  361. }
  362. /**
  363. * @return array (A, Q, R)
  364. */
  365. public function dataProviderForQrDecompositionNonSquareMatricesWithSpecificResults(): array
  366. {
  367. return [
  368. [
  369. 'A' => [
  370. [0],
  371. [0],
  372. ],
  373. 'Q' => [
  374. [1],
  375. [0],
  376. ],
  377. 'R' => [
  378. [0],
  379. ],
  380. ],
  381. [
  382. 'A' => [
  383. [1],
  384. [1],
  385. ],
  386. 'Q' => [
  387. [-0.7071068],
  388. [-0.7071068],
  389. ],
  390. 'R' => [
  391. [-1.414214],
  392. ],
  393. ],
  394. [
  395. 'A' => [
  396. [2, -2, -3],
  397. [0, -6, -1],
  398. [0, 0, 1],
  399. [0, 0, 4],
  400. ],
  401. 'Q' => [
  402. [-1.0, 0.0, 0.0],
  403. [0.0, -1.0, 0.0],
  404. [0.0, 0.0, -1 / \sqrt(17)],
  405. [0.0, 0.0, -4 / \sqrt(17)],
  406. ],
  407. 'R' => [
  408. [-2.0, 2.0, 3.0],
  409. [0.0, 6.0, 1.0],
  410. [0.0, 0.0, -1 * \sqrt(17)],
  411. ],
  412. ],
  413. [
  414. 'A' => [
  415. [1,0,0],
  416. [0,0,0],
  417. [0,0,0],
  418. [0,0,0],
  419. ],
  420. 'Q' => [
  421. [-1,0,0],
  422. [0,1,0],
  423. [0,0,1],
  424. [0,0,0],
  425. ],
  426. 'R' => [
  427. [-1,0,0],
  428. [0,0,0],
  429. [0,0,0],
  430. ],
  431. ],
  432. [
  433. 'A' => [
  434. [3, 7, 6, 4, 5, 8],
  435. [2, 3, 6, 5, 8, 9],
  436. [2, 3, 4, 1, 0, 9],
  437. [3, 7, 6, 7, 7, 3],
  438. [1, 3, 4, 9, 4, 8],
  439. ],
  440. 'Q' => [
  441. [-0.5773503, 0.3086067, -0.2229113, 0.7003493, 0.1767767],
  442. [-0.3849002, -0.5657789, 0.5387023, 0.2150195, -0.4419417],
  443. [-0.3849002, -0.5657789, -0.2414872, -0.3010273, 0.6187184],
  444. [-0.5773503, 0.3086067, -0.2229113, -0.5713376, -0.4419417],
  445. [-0.1924501, 0.4114756, 0.7430376, -0.2150195, 0.4419417],
  446. ],
  447. 'R' => [
  448. [-5.196152, -10.969655, -11.5470054, -10.392305, -10.7772050, -14.818657],
  449. [ 0.000000, 2.160247, -0.3086067, 3.703280, 0.8229512, -3.497543],
  450. [ 0.000000, 0.000000, 2.5634798, 6.687339, 4.6068332, 6.167212],
  451. [ 0.000000, 0.000000, 0.0000000, -2.359071, 0.3624615, 1.394555],
  452. [ 0.000000, 0.000000, 0.0000000, 0.000000, -3.9774756, 5.214913],
  453. ],
  454. ],
  455. ];
  456. }
  457. /**
  458. * @test QR Decomposition invalid property
  459. * @throws \Exception
  460. */
  461. public function testQRDecompositionInvalidProperty()
  462. {
  463. // Given
  464. $A = MatrixFactory::create([
  465. [4, 1, -1],
  466. [1, 2, 1],
  467. [-1, 1, 2],
  468. ]);
  469. $qr = $A->qrDecomposition();
  470. // Then
  471. $this->expectException(Exception\MathException::class);
  472. // When
  473. $doesNotExist = $qr->doesNotExist;
  474. }
  475. /**
  476. * @test QR Decomposition solve incorrect type exception
  477. * @throws \Exception
  478. */
  479. public function testQRDecompositionSolveIncorrectTypeException()
  480. {
  481. // Given
  482. $A = MatrixFactory::create([
  483. [4, 1, -1],
  484. [1, 2, 1],
  485. [-1, 1, 2],
  486. ]);
  487. $qr = $A->qrDecomposition();
  488. // And
  489. $b = new \stdClass();
  490. // Then
  491. $this->expectException(Exception\IncorrectTypeException::class);
  492. // When
  493. $qr->solve($b);
  494. }
  495. }