@@ -494,6 +494,167 @@ public function dataProviderForLuDecomposition(): array
494494 [0 , 1 , 0 , 0 ],
495495 ],
496496 ],
497+ // Near-singular matrix with very small determinant
498+ [
499+ [
500+ [1.0 , 2.0 , 3.0 ],
501+ [4.0 , 5.0 , 6.0 ],
502+ [7.0 , 8.0 , 9.00001 ]
503+ ],
504+ [
505+ [1 , 0 , 0 ],
506+ [0.14285714 , 1 , 0 ],
507+ [0.57142857 , 0.5 , 1 ]
508+ ],
509+ [
510+ [7.0 , 8.0 , 9.00001 ],
511+ [0 , 0.857142857 , 1.71428429 ],
512+ [0 , 0 , -5.0e-6 ]
513+ ],
514+ [
515+ [0 , 0 , 1 ],
516+ [1 , 0 , 0 ],
517+ [0 , 1 , 0 ]
518+ ]
519+ ],
520+ // Near-singular 2x2 matrix
521+ [
522+ [
523+ [2.0 , 4.0 ],
524+ [1.0 , 2.00001 ]
525+ ],
526+ [
527+ [1 , 0 ],
528+ [0.5 , 1 ]
529+ ],
530+ [
531+ [2 , 4 ],
532+ [0 , 1.0E-5 ]
533+ ],
534+ [
535+ [1 , 0 ],
536+ [0 , 1 ]
537+ ]
538+ ],
539+ // Hilbert matrix (ill-conditioned)
540+ [
541+ [
542+ [1.0 , 0.5 , 0.3333333333 ],
543+ [0.5 , 0.3333333333 , 0.25 ],
544+ [0.3333333333 , 0.25 , 0.2 ]
545+ ],
546+ [
547+ [1.0 , 0 , 0 ],
548+ [0.5 , 1.0 , 0 ],
549+ [0.3333333333 , 1 , 1 ]
550+ ],
551+ [
552+ [1.0 , 0.5 , 0.3333333333 ],
553+ [0 , 0.08333333 , 0.08333333335 ],
554+ [0 , 0 , 0.00555556 ]
555+ ],
556+ [
557+ [1 , 0 , 0 ],
558+ [0 , 1 , 0 ],
559+ [0 , 0 , 1 ]
560+ ]
561+ ],
562+ // Diagonally dominant 3x3
563+ [
564+ [
565+ [10.0 , 1.0 , 2.0 ],
566+ [1.0 , 15.0 , 3.0 ],
567+ [2.0 , 1.0 , 20.0 ]
568+ ],
569+ [
570+ [1 , 0 , 0 ],
571+ [0.1 , 1 , 0 ],
572+ [0.2 , 0.05369128 , 1 ]
573+ ],
574+ [
575+ [10.0 , 1 , 2 ],
576+ [0 , 14.9 , 2.8 ],
577+ [0 , 0 , 19.44966443 ]
578+ ],
579+ [
580+ [1 , 0 , 0 ],
581+ [0 , 1 , 0 ],
582+ [0 , 0 , 1 ]
583+ ]
584+ ],
585+ // Diagonally dominant 4x4
586+ [
587+ [
588+ [20.0 , 1.0 , 2.0 , 3.0 ],
589+ [1.0 , 25.0 , 2.0 , 1.0 ],
590+ [2.0 , 1.0 , 30.0 , 2.0 ],
591+ [3.0 , 1.0 , 2.0 , 35.0 ]
592+ ],
593+ [
594+ [1 , 0 , 0 , 0 ],
595+ [0.05 , 1 , 0 , 0 ],
596+ [0.1 , 0.03607214 , 1 , 0 ],
597+ [0.15 , 0.03406814 , 0.05500135 , 1 ]
598+ ],
599+ [
600+ [20 , 1 , 2 , 3 ],
601+ [0 , 24.95 , 1.9 , 0.85 ],
602+ [0 , 0 , 29.73146293 , 1.66933868 ],
603+ [0 , 0 , 0 , 34.42922621 ]
604+ ],
605+ [
606+ [1 , 0 , 0 , 0 ],
607+ [0 , 1 , 0 , 0 ],
608+ [0 , 0 , 1 , 0 ],
609+ [0 , 0 , 0 , 1 ]
610+ ]
611+ ],
612+ // Numerical stability test - large and small numbers
613+ [
614+ [
615+ [1000000.0 , 0.000001 , 1.0 ],
616+ [0.000001 , 1000000.0 , 1.0 ],
617+ [1.0 , 1.0 , 1000000.0 ]
618+ ],
619+ [
620+ [1 , 0 , 0 ],
621+ [0 , 1 , 0 ],
622+ [0.000001 , 1.e-06 , 1 ]
623+ ],
624+ [
625+ [1000000 , 0.000001 , 1 ],
626+ [0 , 1000000 , 1 ],
627+ [0 , 0 , 1.e+06 ]
628+ ],
629+ [
630+ [1 , 0 , 0 ],
631+ [0 , 1 , 0 ],
632+ [0 , 0 , 1 ]
633+ ]
634+ ],
635+ // Numerical stability test - moderate scaling
636+ [
637+ [
638+ [100.0 , 0.01 , 1.0 ],
639+ [0.01 , 100.0 , 1.0 ],
640+ [1.0 , 1.0 , 100.0 ]
641+ ],
642+ [
643+ [1 , 0 , 0 ],
644+ [0.0001 , 1.0 , 0 ],
645+ [0.01 , 0.0099990001 , 1 ]
646+ ],
647+ [
648+ [100 , 0.01 , 1 ],
649+ [0 , 99.999999 , 0.9999 ],
650+ [0 , 0 , 99.980002 ]
651+ ],
652+ [
653+ [1 , 0 , 0 ],
654+ [0 , 1 , 0 ],
655+ [0 , 0 , 1 ]
656+ ]
657+ ],
497658 ];
498659 }
499660
@@ -562,4 +723,53 @@ public function testLUDecompositionSolveIncorrectTypeError()
562723 // When
563724 $ LU ->solve ($ b );
564725 }
726+
727+ /**
728+ * @test Permutation matrix verification - P should be orthogonal
729+ * @dataProvider dataProviderForLUDecomposition
730+ * @param array $A
731+ * @throws \Exception
732+ */
733+ public function testPermutationMatrixProperties (array $ A )
734+ {
735+ // Given
736+ $ A = MatrixFactory::create ($ A );
737+
738+ // When
739+ $ LU = $ A ->luDecomposition ();
740+
741+ // Then - P should be orthogonal (P * P^T = I)
742+ $ P = $ LU ->P ;
743+ $ PT = $ P ->transpose ();
744+ $ PxPT = $ P ->multiply ($ PT );
745+ $ I = MatrixFactory::identity ($ A ->getN ());
746+
747+ $ this ->assertEqualsWithDelta ($ I ->getMatrix (), $ PxPT ->getMatrix (), 0.0001 );
748+
749+ // And - determinant should be ±1
750+ $ det = abs ($ P ->det ());
751+ $ this ->assertEqualsWithDelta (1.0 , $ det , 0.0001 );
752+ }
753+
754+ /**
755+ * @test Round-trip accuracy test - solve Ax=b, then verify A*x ≈ b
756+ * @dataProvider dataProviderForSolve
757+ * @param array $A
758+ * @param array $b
759+ * @throws \Exception
760+ */
761+ public function testRoundTripAccuracy (array $ A , array $ b )
762+ {
763+ // Given
764+ $ A = MatrixFactory::create ($ A );
765+ $ b_vec = new Vector ($ b );
766+ $ LU = $ A ->luDecomposition ();
767+
768+ // When
769+ $ x = $ LU ->solve ($ b );
770+
771+ // Then - A*x should equal b
772+ $ A_times_x = $ A ->multiply ($ x );
773+ $ this ->assertEqualsWithDelta ($ b_vec ->getVector (), $ A_times_x ->getColumn (0 ), 0.0001 );
774+ }
565775}
0 commit comments