diff --git a/src/LinearAlgebra/Eigenvector.php b/src/LinearAlgebra/Eigenvector.php index abed66e87..de9c8de63 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,14 @@ 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λ); diff --git a/tests/LinearAlgebra/Eigen/EigenvectorTest.php b/tests/LinearAlgebra/Eigen/EigenvectorTest.php index 0a6e52817..b57cddd86 100644 --- a/tests/LinearAlgebra/Eigen/EigenvectorTest.php +++ b/tests/LinearAlgebra/Eigen/EigenvectorTest.php @@ -131,6 +131,63 @@ public function dataProviderForEigenvector(): array ]; } + /** + * @test eigenvector can handle numerical precision errors + * @dataProvider dataProviderForPerturbedEigenvalues + * @param array $A + * @param array $E + * @param array $S + */ + public function testEigenvectorsPerturbedEigenvalues(array $A, array $E, array $S) + { + // Perturb E + foreach ($E as $i => $component) { + $E[$i] = $component + (random_int(-1, 1) * 10**-12); + } + + // Given + $A = MatrixFactory::create($A); + $S = MatrixFactory::create($S); + + // When + $eigenvectors = Eigenvector::eigenvectors($A, $E); + + // Then + $this->assertEqualsWithDelta($S, $eigenvectors, 0.0001); + } + + public function dataProviderForPerturbedEigenvalues(): array + { + return [ + [ // Matrix has duplicate eigenvalues. One vector is on an axis. + [ + [2, 0, 1], + [2, 1, 2], + [3, 0, 4], + ], + [5, 1, 1], + [ + [1 / \sqrt(14), 0, \M_SQRT1_2], + [2 / \sqrt(14), 1, 0], + [3 / \sqrt(14), 0, -1 * \M_SQRT1_2], + ] + ], + [ // Matrix has duplicate eigenvalues. no solution on the axis + [ + [2, 2, -3], + [2, 5, -6], + [3, 6, -8], + ], + [-3, 1, 1], + [ + [1 / \sqrt(14), 1 / \M_SQRT3, 5 / \sqrt(42)], + [2 / \sqrt(14), 1 / \M_SQRT3, -4 / \sqrt(42)], + [3 / \sqrt(14), 1 / \M_SQRT3, -1 / \sqrt(42)], + ] + ], + ]; + } + /** * @test eigenvectors throws a BadDataException when the matrix is not square */