Add Moore-Penrose inverse #469

Aweptimum · 2023-05-01T17:57:08Z

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 · 2023-05-01T18:00:24Z

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

Aweptimum · 2023-05-01T18:01:10Z

Oh, right, union types ;-;

markrogoyski · 2023-05-01T18:03:26Z

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 · 2023-05-16T14:09:31Z

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

Beakerboy · 2023-05-16T16:00:06Z

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 · 2023-05-16T16:38:10Z

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 · 2023-05-16T18:27:02Z

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 · 2023-05-16T20:50:38Z

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 · 2023-05-17T17:24:06Z

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

markrogoyski · 2023-05-23T13:55:40Z

@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 · 2023-05-23T14:42:58Z

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 · 2023-07-11T03:00:29Z

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

Mark

Aweptimum · 2023-10-26T13:58:07Z

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 · 2023-10-27T02:48:01Z

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 · 2023-10-27T13:27:55Z

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 · 2023-10-27T13:55:30Z

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.

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.

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`

Calculates the Moore-Penrose inverse of the matrix using SVD

Just an alias for inverse

If it isn't, create a permutation matrix, P, such that SP is diagonal. Multiply V by its inverse to balance the eqn.

After diagonalizing S and ensuring its values are positive, check that its values are sorted. If not, permute it

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.

Not available in php 7.0

to getStandardBasisIndex

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

Since we're checking for non-zero components using Arithmetic::almostEqual, it seems appropriate to pass in the error of the matrix

Adds the case from the pseudoinverse test

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.

Previously was just using the default error of `Support::isZero`

Aweptimum · 2023-12-04T21:27:44Z

@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 · 2023-12-05T07:01:31Z

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.

Aweptimum force-pushed the moore-penrose branch from a610dfa to 0845a30 Compare October 27, 2023 13:13

Aweptimum added 13 commits November 7, 2023 09:36

MatrixCatalog: Add pseudoInverse

67c9097

Add NumericMatrix->pseudoInverse()

50a3147

Calculates the Moore-Penrose inverse of the matrix using SVD

Add NumericDiagonalMatrix->pseudoInverse

8331a97

Just an alias for inverse

Test pseudoInverse in MatrixOperationsTest

48b5b6d

Add problem matrix to SVDTest

93eddf9

SVD: Diagonalize S

efdfcf6

If it isn't, create a permutation matrix, P, such that SP is diagonal. Multiply V by its inverse to balance the eqn.

SVD: Sort S's diagonal

e622e7f

After diagonalizing S and ensuring its values are positive, check that its values are sorted. If not, permute it

SVDTest: add a case for testing tabular sort

aa34e9a

SVD: remove union type from helper

fd4c22c

Not available in php 7.0

Rename SVD->isStandardBasisVector

b52e36a

to getStandardBasisIndex

Aweptimum added 6 commits December 4, 2023 13:26

SVD: allow zero rows in diagonalize

746303b

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

SVD: add error to getStandardBasisIndex

8399246

Since we're checking for non-zero components using Arithmetic::almostEqual, it seems appropriate to pass in the error of the matrix

Add failing case to pseudoInverse test

3ff591e

Add two failures to SVDTestCase

9d3c271

Adds the case from the pseudoinverse test

Numeric: use error in isRectangularDiagonal

77de387

Previously was just using the default error of `Support::isZero`

Aweptimum force-pushed the moore-penrose branch from fc653b3 to 77de387 Compare December 4, 2023 21:21

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add Moore-Penrose inverse #469

Add Moore-Penrose inverse #469

Aweptimum commented May 1, 2023

coveralls commented May 1, 2023 •

edited

Aweptimum commented May 1, 2023

markrogoyski commented May 1, 2023

Aweptimum commented May 16, 2023 •

edited

Beakerboy commented May 16, 2023

markrogoyski commented May 16, 2023

Aweptimum commented May 16, 2023 •

edited

Beakerboy commented May 16, 2023

Aweptimum commented May 17, 2023

markrogoyski commented May 23, 2023

Aweptimum commented May 23, 2023

markrogoyski commented Jul 11, 2023

Aweptimum commented Oct 26, 2023 •

edited

markrogoyski commented Oct 27, 2023

Aweptimum commented Oct 27, 2023 •

edited

Beakerboy commented Oct 27, 2023

Aweptimum commented Dec 4, 2023 •

edited

markrogoyski commented Dec 5, 2023

Add Moore-Penrose inverse #469

Are you sure you want to change the base?

Add Moore-Penrose inverse #469

Conversation

Aweptimum commented May 1, 2023

coveralls commented May 1, 2023 • edited

Aweptimum commented May 1, 2023

markrogoyski commented May 1, 2023

Aweptimum commented May 16, 2023 • edited

Beakerboy commented May 16, 2023

markrogoyski commented May 16, 2023

Aweptimum commented May 16, 2023 • edited

Beakerboy commented May 16, 2023

Aweptimum commented May 17, 2023

markrogoyski commented May 23, 2023

Aweptimum commented May 23, 2023

markrogoyski commented Jul 11, 2023

Aweptimum commented Oct 26, 2023 • edited

markrogoyski commented Oct 27, 2023

Aweptimum commented Oct 27, 2023 • edited

Beakerboy commented Oct 27, 2023

Aweptimum commented Dec 4, 2023 • edited

markrogoyski commented Dec 5, 2023

coveralls commented May 1, 2023 •

edited

Aweptimum commented May 16, 2023 •

edited

Aweptimum commented May 16, 2023 •

edited

Aweptimum commented Oct 26, 2023 •

edited

Aweptimum commented Oct 27, 2023 •

edited

Aweptimum commented Dec 4, 2023 •

edited