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 01/19] 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 02/19] 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λ);
 

From 67c9097ada84e4f75ac6a3de90a0b69b1ff3be0b Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Thu, 23 Mar 2023 15:27:43 -0400
Subject: [PATCH 03/19] MatrixCatalog: Add pseudoInverse

---
 src/LinearAlgebra/MatrixCatalog.php | 30 +++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/src/LinearAlgebra/MatrixCatalog.php b/src/LinearAlgebra/MatrixCatalog.php
index d015003b5..9e5dfa45c 100644
--- a/src/LinearAlgebra/MatrixCatalog.php
+++ b/src/LinearAlgebra/MatrixCatalog.php
@@ -15,6 +15,9 @@ class MatrixCatalog
     /** @var Matrix<T> inverse */
     private $A⁻¹;
 
+    /** @var NumericMatrix pseudoInverse */
+    private $A⁺;
+
     /** @var Reduction\RowEchelonForm */
     private $REF;
 
@@ -43,6 +46,7 @@ class MatrixCatalog
      * DERIVED MATRICES
      *  - transpose
      *  - inverse
+     *  - pseudo-inverse
      **************************************************************************/
 
     // TRANSPOSE
@@ -99,6 +103,32 @@ public function getInverse(): Matrix
         return $this->A⁻¹;
     }
 
+    // PSEUDO-INVERSE
+
+    /**
+     * @param Matrix $A⁺
+     */
+    public function addPseudoInverse(Matrix $A⁺): void
+    {
+        $this->A⁺ = $A⁺;
+    }
+
+    /**
+     * @return bool
+     */
+    public function hasPseudoInverse(): bool
+    {
+        return isset($this->A⁺);
+    }
+
+    /**
+     * @return Matrix
+     */
+    public function getPseudoInverse(): Matrix
+    {
+        return $this->A⁺;
+    }
+
     /**************************************************************************
      * MATRIX REDUCTIONS
      *  - ref (row echelon form)

From 50a3147290e7068c740ca2bc1ccc6acc5793ef10 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Thu, 23 Mar 2023 15:28:37 -0400
Subject: [PATCH 04/19] Add NumericMatrix->pseudoInverse()

Calculates the Moore-Penrose inverse of the matrix using SVD
---
 src/LinearAlgebra/NumericMatrix.php | 59 +++++++++++++++++++++++++++++
 1 file changed, 59 insertions(+)

diff --git a/src/LinearAlgebra/NumericMatrix.php b/src/LinearAlgebra/NumericMatrix.php
index 52970dced..4984678f8 100644
--- a/src/LinearAlgebra/NumericMatrix.php
+++ b/src/LinearAlgebra/NumericMatrix.php
@@ -1599,6 +1599,65 @@ public function inverse(): NumericMatrix
         return $A⁻¹;
     }
 
+    /**
+     * Moore-Penrose inverse, A⁺
+     * Used for non-square or singular matrices
+     * 
+     * Uses the SVD method of construction
+     * Given SVD of A = USVᵀ
+     * A⁺ = VS⁻¹Uᵀ
+     * 
+     * @return NumericMatrix A⁺
+     */
+    public function pseudoInverse(): NumericMatrix
+    {
+        if ($this->catalog->hasPseudoInverse()) {
+            return $this->catalog->getPseudoInverse();
+        }
+        
+        $SVD = $this->svd();
+
+        $U = $SVD->U;
+        $S = $SVD->S;
+        $V = $SVD->V; // already transposed
+        $D = $SVD->D;
+
+        // Manually construct the inverse of S (in case it's singular)
+        $D⁻¹ = [];
+
+        foreach ($D->getVector() as $element)
+        {
+            $D⁻¹[] = abs($element) < 0.0001 ? 0 : 1 / $element;
+        }
+
+        $m = $S->getM();
+        $n = $S->getN();
+
+        $s = MatrixFactory::zero($m, $n)->transpose()->getMatrix();
+
+        for ($i = 0; $i < $n; $i++)
+        {
+            for ($j = 0; $j < $m; $j++)
+            {
+                if ($i === $j) {
+                    $s[$i][$j] = array_shift($D⁻¹);
+                } else {
+                    $s[$i][$j] = 0;
+                }
+            }
+        }
+
+        $S⁻¹ = MatrixFactory::createNumeric($s);
+
+        $Uᵀ  = $U->transpose();
+
+        $A⁺ = $V->multiply($S⁻¹)->multiply($Uᵀ);
+
+        $this->catalog->addPseudoInverse($A⁺);
+
+        return $A⁺;
+    }
+
     /**
      * Cofactor matrix
      * A matrix where each element is a cofactor.

From 8331a97351f9fc69d7f3e00edffb8958a6caf71d Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Fri, 24 Mar 2023 12:30:12 -0400
Subject: [PATCH 05/19] Add NumericDiagonalMatrix->pseudoInverse

Just an alias for inverse
---
 src/LinearAlgebra/NumericDiagonalMatrix.php | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/src/LinearAlgebra/NumericDiagonalMatrix.php b/src/LinearAlgebra/NumericDiagonalMatrix.php
index dc1e424f7..1d16cba85 100644
--- a/src/LinearAlgebra/NumericDiagonalMatrix.php
+++ b/src/LinearAlgebra/NumericDiagonalMatrix.php
@@ -80,4 +80,14 @@ public function inverse(): NumericMatrix
     {
         return MatrixFactory::diagonal(Single::reciprocal($this->getDiagonalElements()));
     }
+
+    /**
+     * pseudoInverse identical to inverse
+     *
+     * @return NumericMatrix
+     */
+    public function pseudoInverse(): NumericMatrix
+    {
+        return $this->inverse();
+    }
 }

From 48b5b6d1e5e02d0194b71cb4bd36eca34b29eee8 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Fri, 24 Mar 2023 12:31:31 -0400
Subject: [PATCH 06/19] Test pseudoInverse in MatrixOperationsTest

---
 .../Matrix/Numeric/MatrixOperationsTest.php   | 78 +++++++++++++++++++
 1 file changed, 78 insertions(+)

diff --git a/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php b/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php
index 5024d942d..b52100631 100644
--- a/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php
+++ b/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php
@@ -367,6 +367,84 @@ public function dataProviderForInverseExceptionDetIsZero(): array
         ];
     }
 
+    /**
+     * Should return the inverse when the matrix is square, so reuse the inverse test dataset 
+     * 
+     * @test         pseudoInverse
+     * @dataProvider dataProviderForInverse
+     * @dataProvider dataProviderForPseudoInverse
+     * @param        array $A
+     * @param        array $A⁺
+     * @throws       \Exception
+     */
+    public function testPseudoInverse(array $A, array $A⁺)
+    {
+        // Given
+        $A  = MatrixFactory::createNumeric($A);
+        $A⁺ = MatrixFactory::createNumeric($A⁺);
+
+        // When
+        $pseudo      = $A->pseudoInverse();
+        $pseudoAgain = $A->pseudoInverse();
+
+        // Then
+        $this->assertEqualsWithDelta($A⁺, $pseudo, 0.001); // Test calculation
+        $this->assertEqualsWithDelta($A⁺, $pseudoAgain, 0.001); // Test class attribute
+    }
+    
+    /**
+     * @return array
+     */
+    public function dataProviderForPseudoInverse(): array
+    {
+        return [
+            [
+                [
+                    [0, 0, 0, 1,0, 0],
+                    [0, 1,0, 0, 1,0],
+                    [0, 0, 1,0, 0, 1],
+                    [0, 0, 0, 0, 0, 0],
+                    [0, 0, -40, 0, 0, -160],
+                    [0, 40, 0, 0, 160, 0]
+                ],
+                [
+                    [0, 0, 0, 0, 0, 0],
+                    [0, 4/3, 0, 0, 0, -1/120],
+                    [0, 0, 4/3, 0, 1/120, 0],
+                    [1, 0, 0, 0, 0, 0],
+                    [0, -1/3, 0, 0, 0, 1/120],
+                    [0, 0, -1/3, 0, -1/120, 0]
+                ]
+            ],
+            [
+                [
+                    [0, 1, 0, 0],
+                    [1, 0, 0, 0],
+                    [0, 0, 1, 0],
+                ],
+                [
+                    [0, 1, 0],
+                    [1, 0, 0],
+                    [0, 0, 1],
+                    [0, 0, 0],
+                ],
+            ],
+            [
+                [
+                    [0, 1, 0, 0],
+                    [1, 1, 0, 0],
+                    [0, 0, 1, 0],
+                ],
+                [
+                    [-1, 1, 0],
+                    [1, 0, 0],
+                    [0, 0, 1],
+                    [0, 0, 0],
+                ],
+            ],
+        ];
+    }
+
     /**
      * @test         cofactor
      * @dataProvider dataProviderForCofactor

From 93eddf9cc1af7c67e7d8c19fd863a6f42b62e860 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Fri, 28 Apr 2023 12:57:14 -0400
Subject: [PATCH 07/19] Add problem matrix to SVDTest

---
 tests/LinearAlgebra/Decomposition/SVDTest.php | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/tests/LinearAlgebra/Decomposition/SVDTest.php b/tests/LinearAlgebra/Decomposition/SVDTest.php
index da9037d90..70bfd04d8 100644
--- a/tests/LinearAlgebra/Decomposition/SVDTest.php
+++ b/tests/LinearAlgebra/Decomposition/SVDTest.php
@@ -211,6 +211,20 @@ public function dataProviderForSVD(): array
                     ],
                 ],
             ],
+            [
+                [
+                    [0,1,0,0],
+                    [1,0,0,0],
+                    [0,0,1,0]
+                ],
+                [
+                    'S' => [
+                        [1,0,0,0],
+                        [0,1,0,0],
+                        [0,0,1,0]
+                    ]
+                ]
+            ],
             // Idempotent
             [
                 [

From efdfcf648a35c7927c5db68b0676bc60cc10f0fd Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 1 May 2023 10:45:21 -0400
Subject: [PATCH 08/19] SVD: Diagonalize S

If it isn't, create a permutation matrix, P, such that SP is diagonal.
Multiply V by its inverse to balance the eqn.
---
 src/LinearAlgebra/Decomposition/SVD.php | 129 ++++++++++++++++++++++++
 1 file changed, 129 insertions(+)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index 3bb807833..d914a5462 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -2,7 +2,9 @@
 
 namespace MathPHP\LinearAlgebra\Decomposition;
 
+use MathPHP\Arithmetic;
 use MathPHP\Exception;
+use MathPHP\Exception\MatrixException;
 use MathPHP\LinearAlgebra\NumericMatrix;
 use MathPHP\LinearAlgebra\MatrixFactory;
 use MathPHP\LinearAlgebra\Vector;
@@ -107,6 +109,14 @@ public static function decompose(NumericMatrix $M): SVD
         // A rectangular diagonal matrix
         $S = $U->transpose()->multiply($M)->multiply($V);
 
+        // If S is non-diagonal, try permuting S to be diagonal
+        if (!$S->isRectangularDiagonal()) {
+            $P = self::diagonalizeColumnar($S);
+
+            $S = $S->multiply($P);            // Permute columns of S
+            $V = $P->inverse()->multiply($V); // Permute corresponding rows of V
+        }
+
         $diag = $S->getDiagonalElements();
 
         // If there is a negative singular value, we need to adjust the signs of columns in U
@@ -123,6 +133,125 @@ public static function decompose(NumericMatrix $M): SVD
         return new SVD($U, $S, $V);
     }
 
+    /**
+     * Returns a permutation matrix, P, such that the product of SP is diagonal
+     * 
+     * @param NumericMatrix $S the matrix to diagonalize
+     * 
+     * @return NumericMatrix a matrix, P, that will diagonalize S (by shuffling cols)
+     */
+    private static function diagonalizeColumnar(NumericMatrix $S): NumericMatrix
+    {
+        if ($S->isRectangularDiagonal()) {
+            return MatrixFactory::identity($S->getN());
+        }
+
+        // Create an identity matrix to store permutations in
+        $P = MatrixFactory::identity($S->getN())->asVectors();
+
+        // Push all zero-columns to the right
+        $cols = $S->asVectors();
+
+        $zeroCols = [];
+
+        foreach ($cols as $i => $colVector)
+        {
+            // Each column should contain 1 non-zero element
+            $isZero = Arithmetic::almostEqual((float) $colVector->l2Norm(), 0);
+
+            $zeroCols[$i] = $isZero ? 0 : 1;
+        }
+
+        arsort($zeroCols, SORT_NUMERIC);
+
+        $zeroMap = array_keys($zeroCols);
+
+        uksort($P, function ($left, $right) use ($zeroMap) {
+            $leftPos = $zeroMap[$left];
+            $rightPos = $zeroMap[$right];
+
+            return $leftPos >= $rightPos;
+        });
+
+        // Only check the columns that contain diagonal entries
+        $rowBound = $S->getM() - 1;
+        $colBound = count($S->getDiagonalElements()) - 1;
+        $cols = $S->submatrix(0,0, $rowBound, $colBound)->asVectors();
+
+        $nonDiagonalValues = [];
+
+        foreach ($cols as $i => $colVector)
+        {
+            // Each column should contain 1 non-zero element
+            $j = self::isStandardBasisVector($colVector);
+
+            if ($j === false) {
+                throw new MatrixException("S Matrix in SVD is not orthogonal:\n" . (string) $S);
+            }
+
+            if ($i === $j) {
+                continue;
+            } else {
+                $nonDiagonalValues[$i] = ['value' => $colVector[$j], 'row' => $j];
+            }
+        }
+
+        // Now create a sort order
+        $order = range(0, $S->getN() - 1);
+
+        foreach ($nonDiagonalValues as $col => $elem)
+        {
+            $row = $elem['row'];
+            $order[$row] = $col;
+        }
+
+        $map = array_flip($order);
+
+        // Need to make column ($i of $nonDiagonalValues) = row ($j)
+        // order = [1=>2, 2=>1, 3=>3]
+        uksort($P, function ($left, $right) use ($map) {
+            $leftPos = $map[$left];
+            $rightPos = $map[$right];
+
+            return $leftPos >= $rightPos;
+        });
+
+        // fromVectors treats the array as column vectors, so the matrix needs to be transposed
+        return MatrixFactory::createFromVectors($P);
+    }
+
+    /**
+     * Checks that a vector has a single non-zero entry
+     * 
+     * @param Vector $v
+     * 
+     * @return int|false The index of the non-zero entry or false if either:
+     *      1. There are multiple non-zero entries
+     *      2. The vector is a zero vector 
+     */
+    private static function isStandardBasisVector(Vector $v): int|false
+    {
+        if ($v->l2Norm() === 0) {
+            return false;
+        }
+
+        // Vectors don't have negative indices
+        $index = -1;
+    
+        foreach ($v->getVector() as $i => $component)
+        {
+            if (!Arithmetic::almostEqual($component, 0)) {
+                if ($index === -1) {
+                    $index = $i;
+                } else { // If we already found a non-zero component, then return false
+                    return false;
+                }
+            }
+        }
+
+        return $index;
+    }
+
     /**
      * Get U, S, or V matrix, or D vector
      *

From e622e7f8fb98d27dbe2d23448ee3d79792d3152b Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 1 May 2023 10:46:41 -0400
Subject: [PATCH 09/19] SVD: Sort S's diagonal

After diagonalizing S and ensuring its values are positive,
check that its values are sorted. If not, permute it
---
 src/LinearAlgebra/Decomposition/SVD.php | 60 +++++++++++++++++++++++++
 1 file changed, 60 insertions(+)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index d914a5462..d7ea17938 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -130,6 +130,16 @@ public static function decompose(NumericMatrix $M): SVD
             $S = $signature->multiply($S);
         }
 
+        // Check the elements are in descending order
+        if (!self::isDiagonalDescending($S)) {
+            $P = self::sortDiagonal($S);
+            /** @var NumericMatrix */
+            $Pᵀ = $P->transpose();
+
+            $S = $Pᵀ->multiply($S)->multiply($P);
+            $U = $P->multiply($U)->multiply($Pᵀ);
+        }
+
         return new SVD($U, $S, $V);
     }
 
@@ -252,6 +262,56 @@ private static function isStandardBasisVector(Vector $v): int|false
         return $index;
     }
 
+    /**
+     * Returns a permutation matrix that sorts its diagonal values in descending order
+     * 
+     * @param NumericMatrix $S singular matrix
+     * 
+     * @return NumericMatrix a permutation matrix such that PᵀSP is diagonal
+     */
+    private static function sortDiagonal(NumericMatrix $S): NumericMatrix
+    {
+        // Get diagonal, pad it by columns, and sort it
+        $diagonal = $S->getDiagonalElements();
+        
+        // Pad
+        $padLength = $S->getN() - count($diagonal);
+
+        $diagonal = array_merge($diagonal, array_fill(0, $padLength, 0)); // Pick 0 because the numbers should all be positive
+
+        // arsort preserves the indices
+        arsort($diagonal, SORT_NUMERIC);
+
+        // ... so we can create a position map from the keys
+        $map = array_keys($diagonal);
+
+        $P = MatrixFactory::identity($S->getM())->asVectors();
+
+        uksort($P, function ($left, $right) use ($map) {
+            $leftPos = $map[$left];
+            $rightPos = $map[$right];
+
+            return $leftPos >= $rightPos;
+        });
+
+        return MatrixFactory::createFromVectors($P);
+    }
+
+    /**
+     * Checks if the elements of a diagonal matrix are in descending order
+     * 
+     * @param NumericMatrix $S the matrix to check
+     * 
+     * @return bool
+     */
+    private static function isDiagonalDescending(NumericMatrix $S): bool
+    {
+        $diagonal = $S->getDiagonalElements();
+        $sorted = array_values($diagonal); rsort($sorted, SORT_NUMERIC);
+
+        return $diagonal === $sorted;
+    }
+
     /**
      * Get U, S, or V matrix, or D vector
      *

From 379bdb44604491df98aa73b2a015075382a7b18a Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 1 May 2023 12:27:19 -0400
Subject: [PATCH 10/19] SVD: Refactor diagonalization for 'tall' matrices

The method of diagonalization permuted columns. However, if the row count exceeds the column count, there's a chance the column permutation won't diagonalize S. To address this, the diagonalization has been refactored to sort on the bigger dimension.
---
 src/LinearAlgebra/Decomposition/SVD.php | 79 +++++++++++++++++--------
 1 file changed, 54 insertions(+), 25 deletions(-)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index d7ea17938..207d497dd 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -111,10 +111,15 @@ public static function decompose(NumericMatrix $M): SVD
 
         // If S is non-diagonal, try permuting S to be diagonal
         if (!$S->isRectangularDiagonal()) {
-            $P = self::diagonalizeColumnar($S);
-
-            $S = $S->multiply($P);            // Permute columns of S
-            $V = $P->inverse()->multiply($V); // Permute corresponding rows of V
+            ['sort'=>$sort, 'P'=>$P] = self::diagonalize($S);
+            // Depending on the value of $sort, we either permute the rows or columns of $S 
+            if ($sort === 'm') {
+                $S = $P->multiply($S);            // Permute rows of S
+                $U = $U->multiply($P->inverse()); // Permute corresponding columns of U
+            } elseif ($sort === 'n') {
+                $S = $S->multiply($P);            // Permute columns of S
+                $V = $P->inverse()->multiply($V); // Permute corresponding rows of V
+            }
         }
 
         $diag = $S->getDiagonalElements();
@@ -148,26 +153,45 @@ public static function decompose(NumericMatrix $M): SVD
      * 
      * @param NumericMatrix $S the matrix to diagonalize
      * 
-     * @return NumericMatrix a matrix, P, that will diagonalize S (by shuffling cols)
+     * @return array{'sort': string, 'P': NumericMatrix} a matrix, P, that will diagonalize S. Multiplication order defined by sort
+     * If 'm', then pre-multiply
+     * If 'n', then post-multiply
      */
-    private static function diagonalizeColumnar(NumericMatrix $S): NumericMatrix
+    private static function diagonalize(NumericMatrix $S): array
     {
         if ($S->isRectangularDiagonal()) {
             return MatrixFactory::identity($S->getN());
         }
 
+        $sort = '';
+        $vecMethod = '';
+        $max = 0;
+        $min = 0;
+
+        if ($S->getM() >= $S->getN()) {
+            $sort = 'm'; // rows
+            $vecMethod = 'asRowVectors';
+            $max = $S->getM();
+            $min = $S->getN();
+        } else {
+            $sort = 'n'; // columns
+            $vecMethod = 'asVectors';
+            $max = $S->getN();
+            $min = $S->getM();
+        }
+
         // Create an identity matrix to store permutations in
-        $P = MatrixFactory::identity($S->getN())->asVectors();
+        $P = MatrixFactory::identity($max)->{$vecMethod}();
 
         // Push all zero-columns to the right
-        $cols = $S->asVectors();
+        $vectors = $S->{$vecMethod}();
 
         $zeroCols = [];
 
-        foreach ($cols as $i => $colVector)
+        foreach ($vectors as $i => $vector)
         {
             // Each column should contain 1 non-zero element
-            $isZero = Arithmetic::almostEqual((float) $colVector->l2Norm(), 0);
+            $isZero = Arithmetic::almostEqual((float) $vector->l2Norm(), 0);
 
             $zeroCols[$i] = $isZero ? 0 : 1;
         }
@@ -184,16 +208,15 @@ private static function diagonalizeColumnar(NumericMatrix $S): NumericMatrix
         });
 
         // Only check the columns that contain diagonal entries
-        $rowBound = $S->getM() - 1;
-        $colBound = count($S->getDiagonalElements()) - 1;
-        $cols = $S->submatrix(0,0, $rowBound, $colBound)->asVectors();
+        $vectors = $S->submatrix(0,0, $min-1, $min-1)->{$vecMethod}();
 
         $nonDiagonalValues = [];
 
-        foreach ($cols as $i => $colVector)
+        /** @var Vector */
+        foreach ($vectors as $i => $vector)
         {
             // Each column should contain 1 non-zero element
-            $j = self::isStandardBasisVector($colVector);
+            $j = self::isStandardBasisVector($vector);
 
             if ($j === false) {
                 throw new MatrixException("S Matrix in SVD is not orthogonal:\n" . (string) $S);
@@ -202,17 +225,17 @@ private static function diagonalizeColumnar(NumericMatrix $S): NumericMatrix
             if ($i === $j) {
                 continue;
             } else {
-                $nonDiagonalValues[$i] = ['value' => $colVector[$j], 'row' => $j];
+                $nonDiagonalValues[$i] = ['value' => $vector[$j], 'j' => $j];
             }
         }
 
         // Now create a sort order
-        $order = range(0, $S->getN() - 1);
+        $order = range(0, $min - 1);
 
-        foreach ($nonDiagonalValues as $col => $elem)
+        foreach ($nonDiagonalValues as $i => $elem)
         {
-            $row = $elem['row'];
-            $order[$row] = $col;
+            $entry = $elem['j'];
+            $order[$entry] = $i;
         }
 
         $map = array_flip($order);
@@ -220,14 +243,20 @@ private static function diagonalizeColumnar(NumericMatrix $S): NumericMatrix
         // Need to make column ($i of $nonDiagonalValues) = row ($j)
         // order = [1=>2, 2=>1, 3=>3]
         uksort($P, function ($left, $right) use ($map) {
-            $leftPos = $map[$left];
-            $rightPos = $map[$right];
+            $leftPos = isset($map[$left]) ? $map[$left] : INF; // sorts in ascending order, so just use inf
+            $rightPos = isset($map[$right]) ? $map[$right] : INF;
 
-            return $leftPos >= $rightPos;
+            return $leftPos <=> $rightPos;
         });
 
-        // fromVectors treats the array as column vectors, so the matrix needs to be transposed
-        return MatrixFactory::createFromVectors($P);
+        $P = MatrixFactory::createFromVectors($P);
+        
+        // fromVectors treats the array as column vectors, so the matrix might need to be transposed
+        if ($sort === 'm') {
+            $P = $P->transpose();
+        }
+
+        return ['sort'=>$sort, 'P' => $P];
     }
 
     /**

From aa34e9a6ebcdc90ef1e64311db156b8b2dc03e34 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 1 May 2023 12:27:43 -0400
Subject: [PATCH 11/19] SVDTest: add a case for testing tabular sort

---
 tests/LinearAlgebra/Decomposition/SVDTest.php | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/tests/LinearAlgebra/Decomposition/SVDTest.php b/tests/LinearAlgebra/Decomposition/SVDTest.php
index 70bfd04d8..d13884185 100644
--- a/tests/LinearAlgebra/Decomposition/SVDTest.php
+++ b/tests/LinearAlgebra/Decomposition/SVDTest.php
@@ -211,6 +211,7 @@ public function dataProviderForSVD(): array
                     ],
                 ],
             ],
+            // Test SVD columnar sort
             [
                 [
                     [0,1,0,0],
@@ -225,6 +226,23 @@ public function dataProviderForSVD(): array
                     ]
                 ]
             ],
+            // Test SVD tabular sort
+            [
+                [
+                    [0, 1, 0],
+                    [1, 0, 0],
+                    [0, 0, 1],
+                    [0, 0, 0],
+                ],
+                [
+                    'S' => [
+                        [1, 0, 0],
+                        [0, 1, 0],
+                        [0, 0, 1],
+                        [0, 0, 0]
+                    ]
+                ]
+            ],
             // Idempotent
             [
                 [

From fd4c22c6040dba5f1238eb99aa779d4748be5a7f Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 1 May 2023 14:13:25 -0400
Subject: [PATCH 12/19] SVD: remove union type from helper

Not available in php 7.0
---
 src/LinearAlgebra/Decomposition/SVD.php | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index 207d497dd..03a78538f 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -218,7 +218,7 @@ private static function diagonalize(NumericMatrix $S): array
             // Each column should contain 1 non-zero element
             $j = self::isStandardBasisVector($vector);
 
-            if ($j === false) {
+            if ($j === -1) {
                 throw new MatrixException("S Matrix in SVD is not orthogonal:\n" . (string) $S);
             }
 
@@ -264,11 +264,11 @@ private static function diagonalize(NumericMatrix $S): array
      * 
      * @param Vector $v
      * 
-     * @return int|false The index of the non-zero entry or false if either:
+     * @return int The index of the non-zero entry or -1 if either:
      *      1. There are multiple non-zero entries
      *      2. The vector is a zero vector 
      */
-    private static function isStandardBasisVector(Vector $v): int|false
+    private static function isStandardBasisVector(Vector $v): int
     {
         if ($v->l2Norm() === 0) {
             return false;
@@ -283,7 +283,7 @@ private static function isStandardBasisVector(Vector $v): int|false
                 if ($index === -1) {
                     $index = $i;
                 } else { // If we already found a non-zero component, then return false
-                    return false;
+                    return -1;
                 }
             }
         }

From b52e36a61838fdd106bc8475518cd62555d82393 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Wed, 17 May 2023 08:20:55 -0400
Subject: [PATCH 13/19] Rename SVD->isStandardBasisVector

to getStandardBasisIndex
---
 src/LinearAlgebra/Decomposition/SVD.php | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index 03a78538f..82d1eb9c4 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -216,7 +216,7 @@ private static function diagonalize(NumericMatrix $S): array
         foreach ($vectors as $i => $vector)
         {
             // Each column should contain 1 non-zero element
-            $j = self::isStandardBasisVector($vector);
+            $j = self::getStandardBasisIndex($vector);
 
             if ($j === -1) {
                 throw new MatrixException("S Matrix in SVD is not orthogonal:\n" . (string) $S);
@@ -260,7 +260,7 @@ private static function diagonalize(NumericMatrix $S): array
     }
 
     /**
-     * Checks that a vector has a single non-zero entry
+     * Checks that a vector has a single non-zero entry and returns its index
      * 
      * @param Vector $v
      * 
@@ -268,7 +268,7 @@ private static function diagonalize(NumericMatrix $S): array
      *      1. There are multiple non-zero entries
      *      2. The vector is a zero vector 
      */
-    private static function isStandardBasisVector(Vector $v): int
+    private static function getStandardBasisIndex(Vector $v): int
     {
         if ($v->l2Norm() === 0) {
             return false;
@@ -282,7 +282,7 @@ private static function isStandardBasisVector(Vector $v): int
             if (!Arithmetic::almostEqual($component, 0)) {
                 if ($index === -1) {
                     $index = $i;
-                } else { // If we already found a non-zero component, then return false
+                } else { // If we already found a non-zero component, then return -1
                     return -1;
                 }
             }

From 746303bcb150a8d38137813362e3bd7efbe17595 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 4 Dec 2023 13:23:35 -0600
Subject: [PATCH 14/19] SVD: allow zero rows in diagonalize

Because 0 can be a singular value, the S matrix can have zero-entries along the diagonal.
To allow this, we first check if the vector magnitude is 0 before calling getStandardBasisIndex
---
 src/LinearAlgebra/Decomposition/SVD.php | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index 82d1eb9c4..4fe7c9540 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -215,8 +215,14 @@ private static function diagonalize(NumericMatrix $S): array
         /** @var Vector */
         foreach ($vectors as $i => $vector)
         {
-            // Each column should contain 1 non-zero element
-            $j = self::getStandardBasisIndex($vector);
+            $ε = $S->getError();
+
+            // Each column should contain up to 1 non-zero element
+            if (Arithmetic::almostEqual($vector->l2Norm(), 0, $ε)) {
+                $j = $i;
+            } else {
+                $j = self::getStandardBasisIndex($vector);
+            }
 
             if ($j === -1) {
                 throw new MatrixException("S Matrix in SVD is not orthogonal:\n" . (string) $S);

From 8399246fd67107a99b421d12d103e538d863758c Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 4 Dec 2023 13:25:06 -0600
Subject: [PATCH 15/19] SVD: add error to getStandardBasisIndex

Since we're checking for non-zero components using Arithmetic::almostEqual, it seems appropriate to pass in the error of the matrix
---
 src/LinearAlgebra/Decomposition/SVD.php | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index 4fe7c9540..17c7d4019 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -221,7 +221,7 @@ private static function diagonalize(NumericMatrix $S): array
             if (Arithmetic::almostEqual($vector->l2Norm(), 0, $ε)) {
                 $j = $i;
             } else {
-                $j = self::getStandardBasisIndex($vector);
+                $j = self::getStandardBasisIndex($vector, $ε);
             }
 
             if ($j === -1) {
@@ -274,7 +274,7 @@ private static function diagonalize(NumericMatrix $S): array
      *      1. There are multiple non-zero entries
      *      2. The vector is a zero vector 
      */
-    private static function getStandardBasisIndex(Vector $v): int
+    private static function getStandardBasisIndex(Vector $v, float $ε): int
     {
         if ($v->l2Norm() === 0) {
             return false;
@@ -285,7 +285,7 @@ private static function getStandardBasisIndex(Vector $v): int
     
         foreach ($v->getVector() as $i => $component)
         {
-            if (!Arithmetic::almostEqual($component, 0)) {
+            if (!Arithmetic::almostEqual($component, 0, $ε)) {
                 if ($index === -1) {
                     $index = $i;
                 } else { // If we already found a non-zero component, then return -1

From 3ff591edf21dce90ea4f0ff7fa821220d4991a40 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Fri, 27 Oct 2023 08:15:24 -0500
Subject: [PATCH 16/19] Add failing case to pseudoInverse test

---
 .../Matrix/Numeric/MatrixOperationsTest.php   | 20 +++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php b/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php
index b52100631..4140ac778 100644
--- a/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php
+++ b/tests/LinearAlgebra/Matrix/Numeric/MatrixOperationsTest.php
@@ -442,6 +442,26 @@ public function dataProviderForPseudoInverse(): array
                     [0, 0, 0],
                 ],
             ],
+
+            // Test Real Failing Matrix
+            [
+                [
+                    [1,0,0,1,0,0],
+                    [0,1,0,0,1,0],
+                    [0,0,1,0,0,1],
+                    [0,0,0,0,0,0],
+                    [0,0,-48.5,0,0,-1681.2],
+                    [0,48.5,0,0,1681.2,0]
+                ],
+                [
+                    [0.5,0,0,0,0,0.],
+                    [0,1.02971,0,0,0,-0.000612482],
+                    [0,0,1.02971,0,0.000612482,0.],
+                    [0.5,0,0,0,0,0.],
+                    [0,-0.0297054,0,0,0,0.000612482],
+                    [0,0,-0.0297054,0,-0.000612482,0.]
+                ]
+            ]
         ];
     }
 

From 9d3c27178fc7939de1df37ebbed2af30030aac36 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Fri, 27 Oct 2023 08:16:07 -0500
Subject: [PATCH 17/19] Add two failures to SVDTestCase

Adds the case from the pseudoinverse test
---
 tests/LinearAlgebra/Decomposition/SVDTest.php | 26 +++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/tests/LinearAlgebra/Decomposition/SVDTest.php b/tests/LinearAlgebra/Decomposition/SVDTest.php
index d13884185..b1ff2da66 100644
--- a/tests/LinearAlgebra/Decomposition/SVDTest.php
+++ b/tests/LinearAlgebra/Decomposition/SVDTest.php
@@ -5,6 +5,7 @@
 use MathPHP\Functions\Support;
 use MathPHP\LinearAlgebra\MatrixFactory;
 use MathPHP\Exception;
+use MathPHP\LinearAlgebra\Eigenvalue;
 use MathPHP\LinearAlgebra\NumericMatrix;
 use MathPHP\LinearAlgebra\Vector;
 use MathPHP\Tests\LinearAlgebra\Fixture\MatrixDataProvider;
@@ -243,6 +244,27 @@ public function dataProviderForSVD(): array
                     ]
                 ]
             ],
+            // Test Real Failing Matrix
+            [
+                [
+                    [1,0,0,1,0,0],
+                    [0,1,0,0,1,0],
+                    [0,0,1,0,0,1],
+                    [0,0,0,0,0,0],
+                    [0,0,-48.5,0,0,-1681.2],
+                    [0,48.5,0,0,1681.2,0]
+                ],
+                [
+                    'S' => [
+                        [1681.89974,0,0,0,0,0],
+                        [0,1681.89974,0,0,0,0],
+                        [0,0,1.41421,0,0,0],
+                        [0,0,0,0.970748,0,0],
+                        [0,0,0,0,0.970748,0],
+                        [0,0,0,0,0,0]
+                    ]
+                ]
+            ],
             // Idempotent
             [
                 [
@@ -334,6 +356,10 @@ public function testSVDProperties(array $A)
         $V = $svd->V;
         $D = $svd->D;
 
+        if (!$svd->getV()->isOrthogonal()) {
+            $p = true;
+        }
+
         // Then U and V are orthogonal
         $this->assertTrue($svd->getU()->isOrthogonal());
         $this->assertTrue($svd->getV()->isOrthogonal());

From 9764083febcba028ef137065b5b1585306126d41 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 4 Dec 2023 10:12:36 -0600
Subject: [PATCH 18/19] Check for floating point in SVD diagonal

If the diagonal has duplicate eigenvalues that are floats, the sort operation can shuffle them around. The strict equivalence then fails when it technically should not.

Instead, we can compare the diagonal entries to the sorted entries using `AlmostEqual` with the matrix error for tolerance.
---
 src/LinearAlgebra/Decomposition/SVD.php | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/LinearAlgebra/Decomposition/SVD.php b/src/LinearAlgebra/Decomposition/SVD.php
index 17c7d4019..01d5ee7b4 100644
--- a/src/LinearAlgebra/Decomposition/SVD.php
+++ b/src/LinearAlgebra/Decomposition/SVD.php
@@ -344,7 +344,15 @@ private static function isDiagonalDescending(NumericMatrix $S): bool
         $diagonal = $S->getDiagonalElements();
         $sorted = array_values($diagonal); rsort($sorted, SORT_NUMERIC);
 
-        return $diagonal === $sorted;
+        // Compare sorted using matrix error (in case duplicate, floating-point eigenvalues) 
+        $n = count($diagonal);
+        for ($i = 0; $i < $n; $i++) {
+            if (!Arithmetic::almostEqual($diagonal[$i], $sorted[$i], $S->getError())) {
+                return false;
+            }
+        }
+
+        return true;
     }
 
     /**

From 77de387a964ad5dc06679b4686481ee50cf19a86 Mon Sep 17 00:00:00 2001
From: Sam <40273116+Aweptimum@users.noreply.github.com>
Date: Mon, 4 Dec 2023 13:07:08 -0600
Subject: [PATCH 19/19] Numeric: use error in isRectangularDiagonal

Previously was just using the default error of `Support::isZero`
---
 src/LinearAlgebra/NumericMatrix.php | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/LinearAlgebra/NumericMatrix.php b/src/LinearAlgebra/NumericMatrix.php
index 4984678f8..e8982761b 100644
--- a/src/LinearAlgebra/NumericMatrix.php
+++ b/src/LinearAlgebra/NumericMatrix.php
@@ -542,7 +542,7 @@ public function isRectangularDiagonal(): bool
         $n = $this->n;
         for ($i = 0; $i < $m; $i++) {
             for ($j = 0; $j < $n; $j++) {
-                if ($i !== $j && !Support::isZero($this->A[$i][$j])) {
+                if ($i !== $j && !Support::isZero($this->A[$i][$j], $this->getError())) {
                     return false;
                 }
             }