From cd45a97139641f7b7c9cbe286146e9a9f15a412a Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Mon, 6 Nov 2023 10:36:47 -0600 Subject: [PATCH 01/10] Replace array_search with floating-point filter array_search seems to fail in most cases when looking for a float in an array of floats. And even if it does find a match, if the key is 0, php evaluates `!0` to true. To find a match, we can instead loop through and compare the numbers with `Arithmetic::almostEqual` and then explicitly check if `$key === false` --- src/LinearAlgebra/Eigenvector.php | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/LinearAlgebra/Eigenvector.php b/src/LinearAlgebra/Eigenvector.php index abed66e87..d71e38eb0 100644 --- a/src/LinearAlgebra/Eigenvector.php +++ b/src/LinearAlgebra/Eigenvector.php @@ -2,6 +2,7 @@ namespace MathPHP\LinearAlgebra; +use MathPHP\Arithmetic; use MathPHP\Exception; use MathPHP\Exception\MatrixException; use MathPHP\Functions\Map\Single; @@ -66,8 +67,15 @@ public static function eigenvectors(NumericMatrix $A, array $eigenvalues = []): foreach ($eigenvalues as $eigenvalue) { // If this is a duplicate eigenvalue, and this is the second instance, the first // pass already found all the vectors. - $key = \array_search($eigenvalue, \array_column($solution_array, 'eigenvalue')); - if (!$key) { + $key = false; + foreach (\array_column($solution_array, 'eigenvalue') as $i => $v) + { + if (Arithmetic::almostEqual($v, $eigenvalue, $A->getError())) { + $key = $i; + break; + } + } + if ($key === false) { $Iλ = MatrixFactory::identity($number)->scalarMultiply($eigenvalue); $T = $A->subtract($Iλ); From 7c75688dcd9fa3cb1a6403f673213b1d143584a0 Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Mon, 6 Nov 2023 10:59:05 -0600 Subject: [PATCH 02/10] Make 8.2 linter happy --- src/LinearAlgebra/Eigenvector.php | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/LinearAlgebra/Eigenvector.php b/src/LinearAlgebra/Eigenvector.php index d71e38eb0..de9c8de63 100644 --- a/src/LinearAlgebra/Eigenvector.php +++ b/src/LinearAlgebra/Eigenvector.php @@ -68,8 +68,7 @@ public static function eigenvectors(NumericMatrix $A, array $eigenvalues = []): // If this is a duplicate eigenvalue, and this is the second instance, the first // pass already found all the vectors. $key = false; - foreach (\array_column($solution_array, 'eigenvalue') as $i => $v) - { + foreach (\array_column($solution_array, 'eigenvalue') as $i => $v) { if (Arithmetic::almostEqual($v, $eigenvalue, $A->getError())) { $key = $i; break; From 22ffbbff823a44ae1c86dccd6a02abc31b9d07fb Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Tue, 31 Oct 2023 14:02:23 -0500 Subject: [PATCH 03/10] Add NumericaMatrix->upperHessenberg Converts a matrix to hessenberg form using householder reflections --- src/LinearAlgebra/NumericMatrix.php | 60 +++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/src/LinearAlgebra/NumericMatrix.php b/src/LinearAlgebra/NumericMatrix.php index 52970dced..7d5a2d107 100644 --- a/src/LinearAlgebra/NumericMatrix.php +++ b/src/LinearAlgebra/NumericMatrix.php @@ -2,6 +2,7 @@ namespace MathPHP\LinearAlgebra; +use MathPHP\Arithmetic; use MathPHP\Functions\Map; use MathPHP\Functions\Support; use MathPHP\Exception; @@ -1478,6 +1479,7 @@ function (NumericMatrix $augmented_matrix, NumericMatrix $matrix) { * - covarianceMatrix * - adjugate * - householder + * - upperHessenberg **************************************************************************/ /** @@ -1886,6 +1888,64 @@ public function householder(): NumericMatrix return Householder::transform($this); } + /** + * Uses householder to convert the matrix to upper hessenberg form + * + * @return NumericMatrix + * + * @throws Exception\MathException + */ + public function upperHessenberg(): NumericMatrix + { + $n = $this->getM(); + + $hessenberg = $this; + $identity = MatrixFactory::identity($n); + + for ($i = 0; $i < $n-2; $i++) + { + $slice = $this->submatrix($i+1, $i, $n-1, $i); + $x = $slice->asVectors()[0]; + $sign = $x[0] >= 0 ? 1 : -1; + $x1 = $sign * $x->l2norm(); + if (Arithmetic::almostEqual($x1, 0, $this->getError())) { + continue; + } + $u = \array_fill(0, $n-$i-1, 0); + $u[0] = $x1; + $u = new Vector($u); + + $v = $u->subtract($x); + $v = MatrixFactory::createFromVectors([$v]); + $vt = $v->transpose(); + + $vvt = $v->multiply($vt); + $vtv = $vt->multiply($v)[0][0]; + $P = $vvt->scalarDivide($vtv); + + $H = MatrixFactory::identity($P->getM()); + $H = $H->subtract($P->scalarMultiply(2)); + + // augment + $offset = $n - $H->getM(); + $elems = $identity->getMatrix(); + for ($i = 0; $i < $H->getM(); $i++) + { + for ($j = 0; $j < $H->getN(); $j++) + { + $elems[$i + $offset][$j + $offset] = $H[$i][$j]; + } + } + + $H = MatrixFactory::create($elems); + + // Multiply + $hessenberg = $H->multiply($hessenberg->multiply($H)); + } + + return $hessenberg; + } + /************************************************************************** * MATRIX VECTOR OPERATIONS - Return a Vector * - vectorMultiply From 9fb768a98acc55c2925ebf9baab8e44c6107d6ae Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Tue, 31 Oct 2023 14:42:02 -0500 Subject: [PATCH 04/10] Test uppserHessenberg --- .../Matrix/Numeric/MatrixOperationsTest.php | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php b/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php index 5024d942d..b66eb0454 100644 --- a/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php +++ b/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php @@ -1485,4 +1485,42 @@ public function dataProviderForRank(): array ], ]; } + + /** + * @test upperHessenberg reduction + * @dataProvider dataProviderForUpperHessenberg + * @param array $A + * @throws \Exception + */ + public function testUpperHessenberg(array $A, array $H) + { + // Given + $A = MatrixFactory::create($A); + $H = MatrixFactory::create($H); + + // When + $actual = $A->upperHessenberg(); + + // Then + $this->assertTrue($actual->isUpperHessenberg()); + $this->assertEqualsWithDelta($H, $actual, 0.0000001); + } + + public static function dataProviderForUpperHessenberg(): \Iterator + { + yield [ + 'A' => [ + [3, 0, 0, 0], + [0, 1, 0, 1], + [0, 0, 2, 0], + [0, 1, 0, 1], + ], + 'H' => [ + [3, 0, 0, 0], + [0, 1, 1, 0], + [0, 1, 1, 0], + [0, 0, 0, 2], + ] + ]; + } } From 44637c37bbfbc118afedf05fcb2542ef1585bc93 Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Mon, 30 Oct 2023 15:18:50 -0500 Subject: [PATCH 05/10] Add qrAlgorithm --- src/LinearAlgebra/Eigenvalue.php | 40 ++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/src/LinearAlgebra/Eigenvalue.php b/src/LinearAlgebra/Eigenvalue.php index c4db0ef4c..692922041 100644 --- a/src/LinearAlgebra/Eigenvalue.php +++ b/src/LinearAlgebra/Eigenvalue.php @@ -2,20 +2,24 @@ namespace MathPHP\LinearAlgebra; +use MathPHP\Arithmetic; use MathPHP\Exception; use MathPHP\Expression\Polynomial; use MathPHP\Functions\Support; +use MathPHP\LinearAlgebra\Decomposition\QR; class Eigenvalue { public const CLOSED_FORM_POLYNOMIAL_ROOT_METHOD = 'closedFormPolynomialRootMethod'; public const POWER_ITERATION = 'powerIteration'; public const JACOBI_METHOD = 'jacobiMethod'; + public const QR_ALGORITHM = 'qrAlgorithm'; private const METHODS = [ self::CLOSED_FORM_POLYNOMIAL_ROOT_METHOD, self::POWER_ITERATION, self::JACOBI_METHOD, + self::QR_ALGORITHM, ]; /** @@ -261,4 +265,40 @@ public static function powerIteration(NumericMatrix $A, int $iterations = 1000): return [$max_ev]; } + + public static function qrAlgorithm(NumericMatrix $A, array &$values = null): array + { + self::checkMatrix($A); + + if ($A->getM() === 1) { + $values[] = $A[0][0]; + return $values; + } + + $values = $values ?? []; + $e = $A->getError(); + $n = $A->getM(); + $H = $A; + $ident = MatrixFactory::identity($n); + $iters = $A->getM() * $A->getN() * 2; + + while (!Arithmetic::almostEqual($H[$n-1][$n-2], 0, $e)) + { + $qr = QR::decompose($H); + $H = $qr->R->multiply($qr->Q); + + // Check if we can shift by the last diagonal element + #$tail = $H[$n-1][$n-2]; + #if (Arithmetic::almostEqual($tail, 0, $e)) { + # $eig = $H[$n-1][$n-1]; + # $shift = $H[$n-2][$n-2]; + # $values[] = $eig; + # $H = $H->add($ident->scalarMultiply($shift)); + # self::qrAlgorithm($H->submatrix(0, 0, $n-2, $n-2), $values); + # break; + #} + } + + return $H->getDiagonalElements(); + } } From 51ba07e60ecb0fa3ebccdc7091cb7f00a4cdd7c0 Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Tue, 31 Oct 2023 09:16:41 -0500 Subject: [PATCH 06/10] test QRAlgorithm eigenvalues --- tests/LinearAlgebra/Eigen/EigenvalueTest.php | 25 ++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/tests/LinearAlgebra/Eigen/EigenvalueTest.php b/tests/LinearAlgebra/Eigen/EigenvalueTest.php index c176283c9..0b1fc95a9 100644 --- a/tests/LinearAlgebra/Eigen/EigenvalueTest.php +++ b/tests/LinearAlgebra/Eigen/EigenvalueTest.php @@ -109,6 +109,31 @@ public function testPowerIteration(array $A, array $S, float $max_abs_eigenvalue $this->assertEqualsWithDelta([$max_abs_eigenvalue], $eigenvalues, 0.0001); } + /** + * @test qrAlgorithm returns the expected eigenvalues + * @dataProvider dataProviderForEigenvalues + * @dataProvider dataProviderForLargeMatrixEigenvalues + * @dataProvider dataProviderForSymmetricEigenvalues + * @param array $A + * @param array $S + * @param float $max_abs_eigenvalue maximum absolute eigenvalue + * @throws \Exception + */ + public function testQRAlgorithm(array $A, array $S) + { + // Given + $A = MatrixFactory::create($A); + + // When + $eigenvalues = Eigenvalue::qrAlgorithm($A); + + sort($S); + sort($eigenvalues); + + // Then + $this->assertEqualsWithDelta($S, $eigenvalues, 0.0001); + } + /** * @test Matrix eigenvalues using powerIterationMethod returns the expected eigenvalues * @dataProvider dataProviderForEigenvalues From edb1bb31fbca6b2e9a5de023774cbf42b58408fd Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Tue, 31 Oct 2023 10:35:05 -0500 Subject: [PATCH 07/10] Shift QR by Wilkinson, refactor a bit The Wilkinson shift accelerates the rate of convergence for the QR algorithm - it's in the static calcShift helper Additionally, moved the iteration to an inner helper function to hide the $values parameter from the public method --- src/LinearAlgebra/Eigenvalue.php | 75 ++++++++++++++++++++++---------- 1 file changed, 52 insertions(+), 23 deletions(-) diff --git a/src/LinearAlgebra/Eigenvalue.php b/src/LinearAlgebra/Eigenvalue.php index 692922041..eb374e92a 100644 --- a/src/LinearAlgebra/Eigenvalue.php +++ b/src/LinearAlgebra/Eigenvalue.php @@ -266,39 +266,68 @@ public static function powerIteration(NumericMatrix $A, int $iterations = 1000): return [$max_ev]; } - public static function qrAlgorithm(NumericMatrix $A, array &$values = null): array + public static function qrAlgorithm(NumericMatrix $A): array { self::checkMatrix($A); - if ($A->getM() === 1) { - $values[] = $A[0][0]; - return $values; + if ($A->isUpperHessenberg()) { + $H = $A; + } else { + $H = $A->upperHessenberg(); } + return self::qrIteration($H); + } + + private static function qrIteration(NumericMatrix $A, array &$values = null): array + { $values = $values ?? []; $e = $A->getError(); $n = $A->getM(); - $H = $A; + + if ($A->getM() === 1) { + $values[] = $A[0][0]; + return $values; + } + $ident = MatrixFactory::identity($n); - $iters = $A->getM() * $A->getN() * 2; - - while (!Arithmetic::almostEqual($H[$n-1][$n-2], 0, $e)) - { - $qr = QR::decompose($H); - $H = $qr->R->multiply($qr->Q); - - // Check if we can shift by the last diagonal element - #$tail = $H[$n-1][$n-2]; - #if (Arithmetic::almostEqual($tail, 0, $e)) { - # $eig = $H[$n-1][$n-1]; - # $shift = $H[$n-2][$n-2]; - # $values[] = $eig; - # $H = $H->add($ident->scalarMultiply($shift)); - # self::qrAlgorithm($H->submatrix(0, 0, $n-2, $n-2), $values); - # break; - #} + + do { + // Pick shift + $shift = self::calcShift($A); + $A = $A->subtract($ident->scalarMultiply($shift)); + + // decompose + $qr = QR::decompose($A); + $A = $qr->R->multiply($qr->Q); + + // shift back + $A = $A->add($ident->scalarMultiply($shift)); + } while (!Arithmetic::almostEqual($A[$n-1][$n-2], 0, $e)); // subdiagonal entry needs to be 0 + + // Check if we can deflate + $eig = $A[$n-1][$n-1]; + $values[] = $eig; + self::qrIteration($A->submatrix(0, 0, $n-2, $n-2), $values); + + return array_reverse($values); + } + + private static function calcShift(NumericMatrix $A): float + { + $n = $A->getM() - 1; + $sigma = .5 * ($A[$n-1][$n-1] - $A[$n][$n]); + $sign = $sigma >= 0 ? 1 : -1; + + $numerator = ($sign * ($A[$n][$n-1])**2); + $denominator = abs($sigma) + sqrt($sigma**2 + ($A[$n-1][$n-1])**2); + + try { + $fraction = $numerator / $denominator; + } catch (\DivisionByZeroError $error) { + $fraction = 0; } - return $H->getDiagonalElements(); + return $A[$n][$n] - $fraction; } } From 51a36e9b5d6e147f46eadc3d5f21dc69de7479c4 Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Wed, 1 Nov 2023 10:02:52 -0500 Subject: [PATCH 08/10] Add qr method to Eigenvector Just gets the eigenvalues from the qr algorithm and feeds them to the eigenvectors method --- src/LinearAlgebra/Eigenvector.php | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/LinearAlgebra/Eigenvector.php b/src/LinearAlgebra/Eigenvector.php index de9c8de63..af0b110f2 100644 --- a/src/LinearAlgebra/Eigenvector.php +++ b/src/LinearAlgebra/Eigenvector.php @@ -10,6 +10,13 @@ class Eigenvector { + public static function qrAlgorithm(NumericMatrix $A) + { + $eigenvalues = Eigenvalue::qrAlgorithm($A); + + return self::eigenvectors($A, $eigenvalues); + } + /** * Calculate the Eigenvectors for a matrix * From 6b445d91b1fc2b74cb9df99b6aaaab2d33e729be Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Wed, 1 Nov 2023 10:03:04 -0500 Subject: [PATCH 09/10] Test Eigenvector::qrAlgorithm --- tests/LinearAlgebra/Eigen/EigenvectorTest.php | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/LinearAlgebra/Eigen/EigenvectorTest.php b/tests/LinearAlgebra/Eigen/EigenvectorTest.php index 0a6e52817..e52f0a6e6 100644 --- a/tests/LinearAlgebra/Eigen/EigenvectorTest.php +++ b/tests/LinearAlgebra/Eigen/EigenvectorTest.php @@ -234,4 +234,25 @@ public function testMatrixEigenvectorInvalidMethodException() // When $A->eigenvectors($invalidMethod); } + + /** + * @test qrAlgorithm + * @dataProvider dataProviderForEigenvector + * @param array $A + * @param array $B + */ + public function testQRAlgorithm(array $A, array $S): void + { + // Given + $A = MatrixFactory::create($A); + $S = MatrixFactory::create($S); + + // When + $eigenvectors = Eigenvector::qrAlgorithm($A); + + // Then + $this->assertEqualsWithDelta($S, $eigenvectors, 0.0001, sprintf( + "Eigenvectors unequal:\nExpected:\n" . (string) $S . "\nActual:\n" . (string) $eigenvectors + )); + } } From 2f08e84d9e7a77e9a96ce2a30cfe12275b9447be Mon Sep 17 00:00:00 2001 From: Sam <40273116+Aweptimum@users.noreply.github.com> Date: Wed, 1 Nov 2023 10:10:24 -0500 Subject: [PATCH 10/10] Add failing eigenvector cases to eigenvalues They work, so the problem must be in the eigenvectors method --- tests/LinearAlgebra/Eigen/EigenvalueTest.php | 30 ++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tests/LinearAlgebra/Eigen/EigenvalueTest.php b/tests/LinearAlgebra/Eigen/EigenvalueTest.php index 0b1fc95a9..af47150dd 100644 --- a/tests/LinearAlgebra/Eigen/EigenvalueTest.php +++ b/tests/LinearAlgebra/Eigen/EigenvalueTest.php @@ -220,6 +220,24 @@ public function dataProviderForEigenvalues(): array [2, 2, 1], 2, ], + [ + [ + [2, 0, 1], + [2, 1, 2], + [3, 0, 4] + ], + [5, 1, 1], + 5, + ], + [ // Matrix has duplicate eigenvalues. no solution on the axis + [ + [2, 2, -3], + [2, 5, -6], + [3, 6, -8], + ], + [-3, 1, 1], + -3 + ], [ [ [1, 2, 1], @@ -276,6 +294,18 @@ public function dataProviderForLargeMatrixEigenvalues(): array ], [4, 3, 2, -2, 1, -1], 4, + ], + [ // Failing case + [ + [2,0,0,0,0,0], + [0,2,0,0,0,1729.7], + [0,0,2,0,-1729.7,0], + [0,0,0,0,0,0], + [0,0,-1729.7,0,2.82879*10**6,0], + [0,1729.7,0,0,0,2.82879*10**6] + ], + [2828791.05765, 0.94235235527, 0.94235235527, 2828791.05765, 2, 0], + 2828791.05765 ] ]; }