From 32f8bff63327b0d7ebf12c3fb8e192bf06077c14 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Tue, 7 Nov 2023 08:29:29 -0600
Subject: [PATCH 1/2] Test Eigenvectors can handle numerical imprecision

To mimic the error found in the QR algorithm, we have to test with matrices that have duplicate eigenvalues and introduce numerical precision errors.

To do this, a list of perturbed eigenvalues is passed to the eigenvectors method. The perturbation is achieved by adding a random +/- offset on an order of magnitude smaller than the default matrix error. This should allow the math to work out fine while causing the floating point comparison to fail.
---
 tests/LinearAlgebra/Eigen/EigenvectorTest.php | 57 +++++++++++++++++++
 1 file changed, 57 insertions(+)

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
      */

From 163c0f93eabcb3f6444220fbcdca5af8ea508cbd 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 2/2] 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 | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

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λ);