Add Moore-Penrose inverse

Aweptimum commented 1 year ago

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

Diagonalize S if it's not rectangular diagonal
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 commented 1 year ago

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

Aweptimum commented 1 year ago

Oh, right, union types ;-;

markrogoyski commented 1 year ago

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 commented 1 year ago

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

Beakerboy commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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

markrogoyski commented 1 year ago

@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 commented 1 year ago

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 commented 11 months ago

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

Mark

Aweptimum commented 8 months ago

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

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.
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 commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

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 commented 7 months ago

@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 commented 7 months ago

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.

markrogoyski / math-php

Add Moore-Penrose inverse #469