`eigs` throws on a basic case (apparently any defective matrix)

[nguyenvukhang]

Describe the bug

Trying to find the eigenvalues and eigenvectors of the matrix [[2, 1], [0, 2]] throws an error. Expected result: eigenvector of [1, 0] with eigenvalue 2.

To Reproduce Steps to reproduce the behavior.

import { eigs } from "mathjs";

const A = [ [5, 2.3], [2.3, 1] ];
const B = [ [2.0, 1.0], [0.0, 2.0] ];

console.log(eigs(A));
// console.log(eigs(B));

Uncommenting the only commented line should throw the Error. (Matrix A is there to make sure that it at least runs ok on some input)

[josdejong]

Thanks for reporting, I would indeed expect the second case to work too. When I execute eigs(B) it throws

Error: Dimension mismatch (1 != 2)

Anyone able to do debug what's going on here? Help would be welcome.

[SinanAkkoyun]

They are bot symmetric, did you try bigger, asymmetric matrices?

[brunoSnoww]

After testing some triangular matrices of different dimensions I was receiving the same error. If in [[2, 1], [0, 2]] we substitute the 0 for 0.0000001 it works, so I guess that doComplexEigs is having a hard time for triangular matrices. However for such matrices the eigenvalues are simply the entries on the diagonal and for getting the eigenvectors we could use another algorithm like SVD. Maybe we could check inside the eigs function if the matrix is triangular, and if so, handle it like that ?

[josdejong]

Thanks for debugging this issue Bruno. So your proposal is to add a special case for triangular matrices? Is this issue is really only occurring for triangular matrices you think? In that case it makes sense to implement a solution for that special case.

[bartekleon]

looks like math for symetric real case is broken. The values are just wrong. a revised version of #3016 should probably fix most of issues, since it is using QR algorithm which loves triangular matrices, and works for symetric as well

[gwhitney]

The observation that the difficulty seems to be with triangular matrices is an interesting one given that the algorithm currently in use begins by performing transformations to make its argument triangular. My guess is that's the part that is stumbling (maybe it's hard to transform something already triangular into triangular form because of something like a division by zero, that sort of thing) and perhaps there's a fix to be had by just testing off the bat if the matrix is triangular, and then proceeding with the rest of the algorithm. I will investigate a bit.

[gwhitney]

Well, I was wrong. The simplest case to look at is [1,1;0,1] from #2984 and I think the characterization there is closer to the mark. That is to say, the existing code does notice that the matrix is triangular, and proceeds immediately to finding the eigenvectors. It also immediately finds the eigenvector (1,0), and crashes looking for another eigenvector -- which does not in fact exist. So the problem lies in that fruitless search. Note [2,1; 0,2] also has just one eigenvector, hence that seems to be the common thread. I will investigate more.

[gwhitney]

In case it is of future interest, a fairly complete reference seems to be http://www.sml.ece.upatras.gr/images/UploadedFiles/efarmosmenes-ypologistikes-methodoi/MatrixEigenvalueProblems-DAHLQUISTch09.pdf

[gwhitney]

OK, the initial error is simply a shape mismatch at in orthogonalComplement in function/matrix/eigs/complexEigs.js. It has the expression:

subtract(v, multiply(divideScalar(dot(w, v), dot(w, w)), w))

but v is a plain vector (1-dim array/matrix of length n) while w is a column vector (i.e., n-by-1 2-dim array/matrix), so the subtraction doesn't really make sense. (It actually produces a value via broadcasting or something like that, but it is not a helpful value and then the calculation goes haywire.)

Reshaping at that point allows the algorithm to proceed further. However, it then quickly runs up against another difficulty: the current algorithm appears to belabor under the misapprehension that all nxn matrices have n independent eigenvectors. While that's almost always true (in the sense that a small random perturbation to any matrix will almost surely have a basis of eigenvectors), it is not literally true for all matrices. The ones that do not have that many eigenvectors are called "defective," and it seems the current algorithm has no chance of succeeding with such matrices. Indeed, with the reshaping fix in place, it crashes again trying to solve a singular system of linear equations; the system is singular precisely because [1,1; 0,1] is defective.

This raises the question of what do we ideally want eigs([1,1,0,0; 0,1,0,0; 0,0,0.5,0;0,0,0,0.5]) for example to return? This is a matrix with two distinct eigenvalues, 1 and 0.5. Right now, eigs calculates the eigenvalues [1, 1, 0.5, 0.5] reflecting the fact that the algebraic multiplicity of each eigenvalue is 2 (the algebraic multiplicities always add up to n when the matrix is nxn). But the geometric multiplicity of eigenvalue 1 is only one, in that there is only one independent eigenvector with eigenvalue 1, e.g. [1,0,0,0]. On the other hand, the geometric multiplicity of eigenvalue 0.5 is two, since both [0,0,1,0] and [0,0,0,1] are eigenvectors. So what should eigs return? if it returns

P1 = {values: [1,1,0.5,0.5], vectors:[1,0,0; 0,0,0; 0,1,0; 0,0,1]}

then you will know that there are just three eigenvectors, but you will not be given the information as to which vector has which eigenvalue: it's ambiguous as to whether the middle column has eigenvalue 1 or 0.5. Conversely, if it returns

P2 = {values: [1,0.5,0.5], vectors:[1,0,0; 0,0,0; 0,1,0; 0,0,1]}

then you will know the eigenvalue for each eigenvector, but you will be left in the dark as to whether 1 has algebraic multiplicity 2 or actually 0.5 has algebraic multiplicity 3 (both are possible for these eigenvalues/eigenvectors with different input matrices). So it seems that we need to extend the return value of eigs to provide more information than it does now (to make the change non-breaking, we want to return both properties (values and vectors) with the exact values they have now in the cases that work, but I don't see it as a problem to add a third property).

So please anyone interested in this issue suggest what values and vectors should be in defective cases and what additional property should be returned (in all cases). To my mind, the vectors property should always be all of the eigenvectors that exist. So to keep the convenient positional correspondence between values and vectors, we should list each eigenvalue with its geometric multiplicity, not algebraic. In other words, use P2 above as the base. Then the additional information needed is the algebraic multiplicity of each eigenvalue. In what format should we return that? It could be an object whose keys are the eigenvalues, e.g. in this case algebraicMultiplicities: {1: 2, 0.5: 2}. But is an object with number keys too weird? We could instead use a row vector positionally corresponding to values at the cost of some redundancy, e.g. algebraicMultiplicities: [2, 2, 2]. Or maybe there is a better format than either? Please weigh in.

When the question of what to actually return is dealt with, I believe I can reasonably easily whip up a PR to handle defective cases and produce the desired return value... Thanks for your thoughts.

[gwhitney]

Oh a decent reference about defective matrices seems to be https://web.mit.edu/18.06/www/Spring17/jordan-vectors.pdf

[gwhitney]

@m93a and I cooked up a likely better suggestion in https://github.com/josdejong/mathjs/issues/3014#issuecomment-1705737957. If we went that way, we could leave the values property exactly as it is with algebraic multiplicities, leave the vectors property alone for now but mark it as deprecated, and add the new property described there (which still needs a name, maybe labeledVectors). Then in the next breaking change we could take out the confusing vectors property. Please let me know what you think, or make another proposal, and I can move forward. Thanks!

[bartekleon]

As a side-note, if possible, we should add some tests for QR decomposition as well. The complex case (at least) uses it. And I have some feelings that QR algorithm might be failing for some of the matrices we feed it. While that might not be a case, testing it doesn't hurt.

[gwhitney]

As a side-note, if possible, we should add some tests for QR decomposition as well.

I think a PR adding tests for QR decomposition would be welcome.

[bartekleon]

@gwhitney also what are your thoughts of using SVD to detect defective matrices / using SVD to calculate eigenvalues / vectors? (I am not very proficient in this part of mathematics so I would love to have second opinion from someone more knowledgable)

[gwhitney]

I am not experienced in the numerical aspects of eigenvalue/eigenvector computation. The current algorithm appears to be operable, so I am reluctant to change it given my lack of expertise. Therefore, my plan is to add detection of defective matrices to the current algorithm, and improve the interface so that it is clear how to obtain the eigenvectors and the which eigenvalue each eigenvector is associated with. Anyone who is experienced in such numerical matters is welcome to submit a PR providing another method of producing eigenvectors and eigenvalues, along with tests that demonstrate that it behaves better/faster/more accurately/etc. than the current algorithm. But that's beyond the scope of what I am able to do.

If you have comments on the matter here, as to what exactly the eigs function should return in the case of a defective matrix, please share so that I can continue to proceed with resolving this issue. Thanks.

[bartekleon]

It would be nice to have some tests, especially edge cases. From what I read, but correct me if I am wrong, these cases are including:

Defective Matrices: Some matrices are called "defective" because they cannot be diagonalized, even if they are square. This occurs when there are not enough linearly independent eigenvectors to form a complete basis for the vector space.

We should be able to detect defective matrices using SVD (tho preferably should find some less intense test)

Zero Eigenvalues: A matrix may have zero as one or more of its eigenvalues. This indicates that the matrix is singular (not invertible). The corresponding eigenvectors are then the null space (kernel) of the matrix.

Should we warn in this case?? I don't know what we should do.

Diagonalizability: A square matrix is diagonalizable if it has a complete set of linearly independent eigenvectors. Not all matrices are diagonalizable. For example, a matrix with repeated eigenvalues may not be diagonalizable.

Same as above.

Also this case, although I don't fully understand what it means : https://en.wikipedia.org/wiki/QR_algorithm#Finding_eigenvalues_versus_finding_eigenvectors

I will be reading / learning about decompositions / eigenvalues/vectors this week and hopefully will find some more information / fixes to current code. As of today I have QR algorithm working which can be ported into the project in case current one will have issues (we need some testing, especially strange / edge cases for current QR implementation).

[josdejong]

Thanks a lot for debugging this Glen. So if I understand it correctly, the solution will be to handle defective matrices in eigs. It makes sense to combine this with an improved API for eigs discussed in #3014. Thanks for your offer @gwhitney to come up with a PR to handle defective cases!

[bartekleon]

TLDR: Having SVD, computing EVD should be trivial

@gwhitney @josdejong It seems that algorithms like QR, LU, EVD (eigenvalue decomposition) are prone to having some cases which are not computable (as we already know). From what I have read, SVD (singular value decomposition) will always work resulting in 3 matrices - U, Σ (sigma), and V* (complex conjugate transpose of V). From there we can just look at the results in sigma, and if any of the values is zero or close to zero (precision) then the eigenvalues and eigenvectors of given matrix cannot be computed.

The small problem as of now is understanding how to write SVD algorithm. There seem to be different algorithms, one looking quite tedious and annoying to write (basically getting eigenvalues and eigenvectors of A*A and AA* - which, from what I understand, is always possible and constructing SVD from it using simple matrix multiplication), the other one making something like repeated power iteration to get SVD.

Upside of it is that we will have another feature of math.js being SVD decomposition. Downside is that we need to code it and it will hurdle the performance.

edit: It seems in order to have SVD we still need QR/householder decomposition + finding eigenvalues and eigenvectors (yes, this is kinda cyclic)

I think everything would look like this:

_EVD - to calculate eigenvalues and vectors (this is for where we know the results exist) SVD - calculates singular value decomposition EVD - calculates eigenvalues and vectors / throws if i can't / returns that matrix is singular / defective

EVD algorithm:

• run SVD and check if any of singular values is ~ 0
• if zero - throw / return
• else run _EVD

[bartekleon]

It sure took me a while, but I finally got eigenvalues and vectors working.... Now only thing left is to test it... A lot... Also found two small places that can be optimised (in other parts of code)

[bartekleon]

During coding found another (related) bug in the code. Schur decomposition is not working properly for some cases: [ [0, 1], [-1, 0] ] just keeps spinning around, while the result should be: Q = {{i,-i},{1,1}} T = {{-i,0},{0,i}} see: https://www.dcode.fr/matrix-schur-decomposition

and validation

https://www.wolframalpha.com/input?i2d=true&i=%7B%7Bi%2C-i%7D%2C%7B1%2C1%7D%7D%7B%7B-i%2C0%7D%2C%7B0%2Ci%7D%7D%7B%7BDivide%5B1%2C2i%5D%2CDivide%5B1%2C2%5D%7D%2C%7B-Divide%5B1%2C2i%5D%2CDivide%5B1%2C2%5D%7D%7D

(I haven't manage to make wolframalpha calculate inverse of matrix so I calculated it beforehand) To validate it's true: https://www.wolframalpha.com/input?i2d=true&i=inverse+%7B%7Bi%2C-i%7D%2C%7B1%2C1%7D%7D

Edit: It seems that for some reason both of the results are not proper schur decomposition... According to wikipedia (well, hard to tell how much this is trustworthy now...) Q should be unitary (which in our case is, in case of dcode it is not) T (or U) should be upper triangular (in our case it is not, in case of dcode it is)

[dcode decomposition looks better though, since it has the eigenvalues on it's diagonal and works with inverse, while our/wolframlpha only works if we use complex conjugate (no eigenvalues on diagonal tho)]

edit 2: it does seem like both answers are wrong. Can be tested in python:

edit 3: found quite nice solution for Schur decomposition https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/Eigenvalues/ComplexSchur.h?ref_type=heads

This should handle some additional cases in testing (still failing a few tests)

[gwhitney]

During coding found another (related) bug in the code. Schur decomposition is not working properly for some cases: [ [0, 1], [-1, 0] ] just keeps spinning around,

Thanks. Before this issue is closed, either that should be solved or be filed as its own issue if not.

[bartekleon]

During coding found another (related) bug in the code. Schur decomposition is not working properly for some cases: [ [0, 1], [-1, 0] ] just keeps spinning around,

Thanks. Before this issue is closed, either that should be solved or be filed as its own issue if not.

Well, you most likely use Schur decomposition to calculate eigenvalues, so I feel it is very much related. Or to be precise, if i had working Schur decomposition, Eigenvalue/vector decomposition should just work

Unfortunately from what I researched, we need Schur decomposition using QR decomposition with shifting, otherwise it just won't work. Normal QR decomposition using householder or even Givens rotations won't work for these rare cases.

@gwhitney You can actually read about the issue I am currently having here: https://faculty.ucmerced.edu/mhyang/course/eecs275/lectures/lecture17.pdf page 23

I feel like it should be:

• Hessenberg decomposition (to improve convergence and stability)
• Schur decomposition using QR decomposition with shifting (see: https://gitlab.com/libeigen/eigen/-/blob/master/Eigen/src/Eigenvalues/ComplexSchur.h?ref_type=heads)
• From Schur decomposition eigenvalues are trivial (diagonal of uppertriangular matrix from Schur)
• Then just power inversion method or something similar to get eigenvectors

I think the Eigen CPP project is the best to put reference on. The results from there are correct and code is in readable language

(I think I am almost finished rewriting this schur decomposition, though it is quite tough)

[josdejong]

So it looks like we have two possible approaches here:

1. handle defective matrices in the current implementation
2. replace the implementation with an SVD algorithm, and from what I understand from Bart osz this comes with it's own challenges 😅.

It sounds promising to work on SVD as a long term solution, and we can work on that in the sideline, but I prefer to focus on trying to fix the issue around defective matrices in the current implementation first. I expect that will require less effort than a complete new implementation, and fixing on short term would be nice.

[bartekleon]

I have almost finished fixing the current implementation:

• fixed Schur decomposition
• get eigenvalues correctly for ALL cases

I am still struggling with eigenvectors thought. There are so many cases...

• technically the Q matrix of Schur decomposition (A=QUQ^T) should have eigenvectors as columns, but it is not always the case it seems...
• then there is inverse power method / power method - they only give one pair of eigenvalue / vector
• inverse power method with shifting - fails when determinant is 0 (which is almost always the case....)

If anyone has some iterative methods / resources to get eigenvectors when you already have eigenvalues I would be grateful

Also @josdejong it appeared that I misunderstood what SVD does. So it only helps in classifying matrices, but does not really help in finding eigenvalues/vectors (you get eigenvalues of AA and AA which is completely something else)

[bartekleon]

Small update @josdejong : It seems while my implementation of Schur does work on 2x2, it still fails for some bigger matrices - that being said, we really need to rewrite / check all tests... It appears they were smacked to "just pass" while the results are incorrect... While the first test in Schur does hold the A = QUQ^T equation, it misses the fact that Q should be unitary matrix... So I will take a look into that a little more.

In fact I will focus on making Schur implementation work first, before working on eigenvectors (since Schur is a base algorithm for eigen decomposition)

(I will check all decompositions, especially for bigger matrices if they are behaving similar to the python scipy implementation [which uses Lapack AFAIK and is proven to give correct results])

[gwhitney]

Sorry I was just waiting on a decision on what the correct API for eigs should be before repairing the current algorithm. Has that been decided? Whatever algorithm is used, the API must be changed to handle all cases -- it's really inadequate for defective matrices. I think we should fix the API and the current implementation before embarking on any major internal rewrite. Jos, did you get enough feedback to select a specific API for eigs going forward? When we have a clear design, I am happy to make it so.

[bartekleon]

Sorry I was just waiting on a decision on what the correct API for eigs should be before repairing the current algorithm. Has that been decided? Whatever algorithm is used, the API must be changed to handle all cases -- it's really inadequate for defective matrices. I think we should fix the API and the current implementation before embarking on any major internal rewrite. Jos, did you get enough feedback to select a specific API for eigs going forward? When we have a clear design, I am happy to make it so.

As for API, whenever you decide on it I can just plug it in, this shouldn't be an issue really so don't worry. I am just working on trying to fix / workout actually working implementation. AFAIK defective matrices still have eigenvalues, so it should not fail. Eigenvectors on the other hand are not always required to exist

example of how scipy handles some defective matrix case.

[josdejong]

Thanks for the update @bartekleon .

@gwhitney shall we discuss the details of the API further in #3014?

[bartekleon]

Writing everything from scratch made me finally understand what is going on and what we are "sitting on".

So we do support all non defective cases. We support some defective cases to some degree.

First thing that would need to be done (to make stuff a little less confusing) is normalise the vectors - basically divide each vector by it's last element so it is 1 (seems to be done universally, and makes sense since we getting something like ax = y, instead of ax = by)

As for inverseIterate throwing error - it is happening when we find one vector via usolve and it is passed further (because it looks like [[a], [b]] instead [a, b]). It is very easily fixable, but it doesn't resolve any issue sadly

As for defective cases: apparently if we calculate (A - eI) and get rank of it we can determine if we get linearly independent eigenvectors or not. But as for calculation of these, still don't know universal algorithm

Side note: it seems that our eigs algorithm does Schur decomposition, so my suggestion is to move most of the code to schur.js Although I would need to test it a lot since with so many different cases it might have been a fluke...

[gwhitney]

First thing that would need to be done (to make stuff a little less confusing) is normalise the vectors - basically divide each vector by it's last element so it is 1 (seems to be done universally, and makes sense since we getting something like ax = y, instead of ax = by)

I beg to differ, that is not the first thing that should be done in response to this bug; it is a thing we might want to do at some point, there's a very reasonable argument for it, but it is a breaking change to working code. Definitely the PR that addresses this bug should not change any working test cases for eigs, only add new test cases that failed before but now pass. Then we can open a new issue about eigs normalization. For that, my advice would be a two-step process, first allowing normalization by a flag so it is non-breaking, and later removing the unnormalized version in a breaking change, but Jos is generally in favor of a more direct breaking way of moving forward, so we would probably do it that way. But that's a separate discussion to have at a separate time.

As for inverseIterate throwing error - it is happening when we find one vector via usolve and it is passed further (because it looks like [[a], [b]] instead [a, b]). It is very easily fixable, but it doesn't resolve any issue sadly

As for defective cases: apparently if we calculate (A - eI) and get rank of it we can determine if we get linearly independent eigenvectors or not.

Agreed on both counts. And in fact we are already doing part of the rank computation (for the submatrices of the almost-triangular form of A that correspond to each eigenvalue in turn). We certainly find out when those submatrices are not full rank (well, by crashing at the moment, but still). So we know when any eigenvector has lower geometric multiplicity than algebraic. So that's why I say, once we can finally settle on a desired eigs interface, it will not be too much work to fix the current implementation.

[josdejong]

Fixed in v12.0.0 now via #3037