Replace array_search with floating-point filter in `Eigenvector::eige…

…nvectors` (#473)

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

* 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

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

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