Add Moore-Penrose inverse

[Aweptimum]

Closes #468 and fixes the edge-case found in the SVD implementation. Helper methods were added to the SVD class to do two things:

1. Diagonalize S if it's not rectangular diagonal
2. Sort the entries of S if they aren't in descending order

You can remove the commit that sorts the values in descending order if you wish. It was done to fix a technicality that isn't represented in the test data, and I intuited the procedure (not sure it works in all cases). It worked on paper with a dummy S-matrix, but I could not figure out how to create a matrix that yielded an S-matrix with unordered values when decomposed.

[coveralls]

coverage: 99.638% (-0.3%) from 99.921%
when pulling 77de387 on Aweptimum:moore-penrose
into da47751 on markrogoyski:develop.

[Aweptimum]

Oh, right, union types ;-;

[markrogoyski]

FYI: I have to fix the static analysis in dev. It may still be failing on an unrelated PLS issue. Ignore any PHPStan failures for that for the time being.

[Aweptimum]

I've pushed the change, seems the CI is happy with the syntax now (all unit tests passed) and only the formatting failed.

[Beakerboy]

I see a push where isStandardBasisVector() returns an int. Is there precedent in other places of MathPHP to return more than just True/False from a isFoo() named method? I seem to remember us using the pattern:

if isFoo($bar) {
    $j = $bar.getFoo();
}

[markrogoyski]

I see a push where isStandardBasisVector() returns an int. Is there precedent in other places of MathPHP to return more than just True/False from a isFoo() named method? I seem to remember us using the pattern:

if isFoo($bar) {
    $j = $bar.getFoo();
}

Preferably the method name should indicate what it returns. So a query phrased as a question should return a boolean and only that. If you want to get a value, just ask for the value. If it is part of the public interface, getValueName() is fine.

[Aweptimum]

I could rename it to getStandardBasisIndex since that's what it's grabbing. I just wanted to avoid looping twice - once to check, once to find the index and return it.

[Beakerboy]

If it’s computationally intensive, the isFoo() method could set a class attribute, and the getter can check if it is set and only calculate it if it is not.

[Aweptimum]

I chose to rename it. Let me know if there's anything else that needs changing.

[markrogoyski]

@Aweptimum,

Thank you for the pull request implementing matrix pseudo inverse and fixing SVD behavior.

Give me some time to learn about the math here and I'll come up with some additional test cases to verify the changes. Thanks for your patience.

Mark

[Aweptimum]

Cool, thanks Mark! If you have any questions let me know. For the Pseudo-Inverse it should deal with any singular matrix, but I'm not sure what sort of matrix the SVD couldn't handle before (was it any permutation matrix or just that one?)

[markrogoyski]

Hey. I haven't forgotten about this. I will get to this soon. Sorry for the wait.

Mark

[Aweptimum]

I've been using my forked branch in a project and have got two failing test cases:

1. According to Wolfram, you can have a 0 in the diagonal as long as it's in the trailing side, so I need to modify the standard basis helper to be fine with zero-vectors. Not sure what to rename it to.
2. I have some other matrix that generates a very wonky inverse. The values are correct, but some of them are in the wrong column and have the wrong sign.
   I'm working on fixing both cases now

Edit:
So with case 2, it seems that the eigenvectors matrix, U (MM^T), is incorrect. This is MM^T):

[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, 2828785.69, 0]
[0, 1729.7, 0, 0, 0, 2828785.69]

These are its eigenvectors according to wolfram:

v1 = [0,  0.00061146305, 0, 0, 0, 1]
v2 = [0, 0, -0.00061146305, 0, 1, 0]
v3 = [1, 0, 0, 0, 0, 0]
v4 = [0, -1635.422, 0, 0, 0, 1]
v5 = [0, 0, 1635.422, 0, 1, 0]
v6 = [0, 0, 0, 1, 0, 0]

This is the value of $MMᵀ->eigenvectors() inside of SVD->decompose:

[-0, -0, 1, -0, -0, 0]
[-0.00043237025078941, -0.00043237025078941, 0, -0.70710664899714, -0.70710664899714, 0]
[-0.00043237025078941, 0.00043237025078941, 0, -0.70710664899714, 0.70710664899714, 0]
[-0, -0, 0, -0, -0, 1]
[0.70710664899714, -0.70710664899714, 0, -0.00043237025078941, 0.00043237025078941, 0]
[-0.70710664899714, -0.70710664899714, 0, 0.00043237025078941, 0.00043237025078941, 0]

I added the matrix to the EigenvalueTest's LargeMatrixProvider and I think the calculated values are correct - the power method spits out 2828791.05765 which is consistent with the eigenvalues in that wolfram link.

Should I commit the test cases and push them so you guys can examine them?

[markrogoyski]

Hi @Aweptimum,

Thank you for continuing to work on this. Yes, please push any test cases you have. Unit tests are useful for not only validation but also understanding how things should work.

Thanks again for your continued support in making MathPHP better.

Mark

[Aweptimum]

Alright, pushed. I'm really not sure why this particular matrix fails, most of the ones in my code are of the same form:

[
    [1, 0, 0, 1, 0, 0],
    [0, 1, 0, 0, 1, 0],
    [0, 0, 1, 0, 0, 1],
    [0, -z1, y1, 0, y2, -z2]
    [z1, 0, -x1, z2, 0, -x2]
    [-y1, x1, 0, -y2, x2, 0]
]

It represents a bearing couple on a shaft, with x, y, and z being the displacements of the bearings wrt the origin. I'm using geometry transforms so that this matrix is always evaluated as if it were aligned with the x-axis, so the z and y displacements are always 0. I've successfully solved 5 systems with the same sort of matrix, but this particular one is breaking it:

[
    [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]
]

@Beakerboy did you write the eigenvector solver? Would you be able to provide any insight? I read a bit on the power method and I don't understand how the singular values of a 6x6 matrix can be calculated if the one available eigenvalue algorithm spits out a single value.

[Beakerboy]

If I remember, the power method returns the largest eigenvalue. You can then use this value to reformulate the matrix in a way that removes it…then rerun the power method to get the new largest value.

[Aweptimum]

@markrogoyski I think I tracked the problem down. I experimented with removing zero rows from my A matrix and corresponding zero entries in the b vector in my project. The resulting inverse was still incorrect, and I tracked it down to the isDiagonalDescending. It was sorting the diagonal entries and then using a strict compare against the original diagonal elements. Because the case I introduced has multiple duplicate singular values that are floating points, the sorted values didn't exactly match the original and the method failed.

I replaced the strict array comparison with a loop that uses Arithmetic::almostEquals to check if the elements are within the matrix error, and all the SVD/Pseudoinverse tests pass now!

[markrogoyski]

Hi @Aweptimum,

Thank you for continuing to work to improve the SVD function.

I wrote a bunch of new unit tests on SVD and most of them are passing with your changes. However, I did get one to fail. Maybe you can take a look at it. Here is how to reproduce it. Just add this test case to the data provider.

            [
                [
                    [1, 0, 1, 0],
                    [0, 1, 0, 1],
                ],
                [
                    'S' => [
                        [\sqrt(2), 0, 0, 0],
                        [0, \sqrt(2), 0, 0],
                    ],
                ],
            ],

Exception that is thrown:

1) MathPHP\Tests\LinearAlgebra\Decomposition\SVDTest::testBug
MathPHP\Exception\MatrixException: S Matrix in SVD is not orthogonal:
[1, 1, 0, 0]
[1, -1, 0, 0]

/src/LinearAlgebra/Decomposition/SVD.php:228
/src/LinearAlgebra/Decomposition/SVD.php:114
/src/LinearAlgebra/NumericMatrix.php:2922
/tests/LinearAlgebra/Decomposition/SVDTest.php:488

This is where I got this test case from, which shows the steps and solution: https://atozmath.com/example/MatrixEv.aspx?q=svd&q1=E2

And for reference, here is the same matrix in R and Python.

R

> A <- rbind(c(1, 0, 1, 0), c(0, 1, 0, 1))
> svd(A)
$d
[1] 1.414214 1.414214

$u
     [,1] [,2]
[1,]    1    0
[2,]    0    1

$v
          [,1]      [,2]
[1,] 0.7071068 0.0000000
[2,] 0.0000000 0.7071068
[3,] 0.7071068 0.0000000
[4,] 0.0000000 0.7071068

Python

In [2]: A = [[1, 0, 1, 0], [0, 1, 0, 1]]

In [3]: svd = scipy.linalg.svd(A)

In [6]: svd
Out[6]:
(array([[1., 0.],
        [0., 1.]]),
 array([1.41421356, 1.41421356]),
 array([[ 0.70710678,  0.        ,  0.70710678,  0.        ],
        [ 0.        ,  0.70710678,  0.        ,  0.70710678],
        [-0.70710678,  0.        ,  0.70710678,  0.        ],
        [ 0.        , -0.70710678,  0.        ,  0.70710678]]))

By the way, here are some of the new tests I'm planning to add that do pass, if you want to add them to your PR.

           [
                [
                    [3, 2, 2],
                    [2, 3, -2],
                ],
                [
                    'S' => [
                        [5, 0, 0],
                        [0, 3, 0],
                    ],
                ],
            ],
            [
                [
                    [2, 4],
                    [1, 3],
                    [0, 0],
                    [0, 0],
                ],
                [
                    'S' => [
                        [5.4649857, 0],
                        [0, 0.3659662],
                        [0, 0],
                        [0, 0],
                    ],
                ],
            ],
            [
                [
                    [0, 0, 0],
                    [0, 0, 0],
                    [0, 0, 0],
                ],
                [
                    'S' => [
                        [0, 0, 0],
                        [0, 0, 0],
                        [0, 0, 0],
                    ],
                ],
            ],
            [
                [
                    [1, 0, 0],
                    [0, 0, 0],
                    [0, 0, 0],
                ],
                [
                    'S' => [
                        [1, 0, 0],
                        [0, 0, 0],
                        [0, 0, 0],
                    ],
                ],
            ],
            [
                [
                    [2, 0, 0],
                    [0, 0, 0],
                    [0, 0, 0],
                ],
                [
                    'S' => [
                        [2, 0, 0],
                        [0, 0, 0],
                        [0, 0, 0],
                    ],
                ],
            ],
            [
                [
                    [0, 1, 0],
                    [0, 0, 0],
                    [0, 0, 0],
                ],
                [
                    'S' => [
                        [1, 0, 0],
                        [0, 0, 0],
                        [0, 0, 0],
                    ],
                ],
            ],
            [
                [
                    [-4, -7],
                    [1, 4],
                ],
                [
                    'S' => [
                        [9, 0],
                        [0, 1],
                    ],
                ],
            ],
            [
                [
                    [3, 0],
                    [4, 5],
                ],
                [
                    'S' => [
                        [\sqrt(45), 0],
                        [0, \sqrt(5)],
                    ],
                ],
            ],
            [
                [
                    [0, 1, 0, 0],
                    [0, 0, 2, 0],
                    [0, 0, 0, 3],
                    [0, 0, 0, 0],
                ],
                [
                    'S' => [
                        [3, 0, 0, 0],
                        [0, 2, 0, 0],
                        [0, 0, 1, 0],
                        [0, 0, 0, 0],
                    ],
                ],
            ],
            [
                [
                    [0, 1, 0, 0],
                    [0, 0, 2, 0],
                    [0, 0, 0, 3],
                    [1 / 60_000, 0, 0, 0],
                ],
                [
                    'S' => [
                        [3, 0, 0, 0],
                        [0, 2, 0, 0],
                        [0, 0, 1, 0],
                        [0, 0, 0, 1.666667e-05],
                    ],
                ],
            ],
            [
                [
                    [4, 0],
                    [3, -5],
                ],
                [
                    'S' => [
                        [\sqrt(40), 0],
                        [0, \sqrt(10)],
                    ],
                ],
            ],

I'll keep trying other various unit test cases. And thanks again for continuing to move forward with this.