New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
eigs
throws on a basic case (apparently any defective matrix)
#2879
Comments
Thanks for reporting, I would indeed expect the second case to work too. When I execute
Anyone able to do debug what's going on here? Help would be welcome. |
They are bot symmetric, did you try bigger, asymmetric matrices? |
After testing some triangular matrices of different dimensions I was receiving the same error. If in |
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. |
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 |
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. |
Well, I was wrong. The simplest case to look at is |
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 |
OK, the initial error is simply a shape mismatch at in orthogonalComplement in function/matrix/eigs/complexEigs.js. It has the expression:
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
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
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 ( So please anyone interested in this issue suggest what 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. |
eigs
throws on a basic case.eigs
throws on a basic case (apparently any defective matrix)
Oh a decent reference about defective matrices seems to be https://web.mit.edu/18.06/www/Spring17/jordan-vectors.pdf |
@m93a and I cooked up a likely better suggestion in #3014 (comment). If we went that way, we could leave the |
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. |
I think a PR adding tests for QR decomposition would be welcome. |
@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) |
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. |
It would be nice to have some tests, especially edge cases. 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). |
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. edit: I think everything would look like this: _EVD - to calculate eigenvalues and vectors (this is for where we know the results exist) EVD algorithm:
|
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) |
During coding found another (related) bug in the code. Schur decomposition is not working properly for some cases: and validation (I haven't manage to make wolframalpha calculate inverse of matrix so I calculated it beforehand) 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...) [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 This should handle some additional cases in testing (still failing a few tests) |
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:
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) |
So it looks like we have two possible approaches here:
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. |
I have almost finished fixing the current implementation:
I am still struggling with eigenvectors thought. There are so many cases...
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 A*A which is completely something else) |
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... 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]) |
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 |
Thanks for the update @bartekleon . @gwhitney shall we discuss the details of the API further in #3014? |
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: Side note: |
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.
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. |
Previously, attempting to take the `eigs` of any defective matrix was doomed to fail in an attempt to solve a singular linear system. This PR detects the situation (as best as it can given the inherent numerical instability of the current methods used) and handles it. Note that in such cases, it's not possible to return a square matrix whose columns are the eigenvectors corresponding to the returned eigenvalues. In light of that fact and issue josdejong#3014, this PR also changes the return value of `eigs` so that the eigenvectors are passed back in a property `eigenvectors` which is an array of plain objects `{value: e, vector: v}`. Note that this PR makes the ancillary changes of correcting the spelling of the filename which was "realSymetric.js," and replacing the now-unnecessary auxiliary function "createArray" therein with `Array(size).fill(element)`. The rationale for performing these changes not strictly related to the issues at hand is that this file is rarely touched and with the level of maintenance hours we have at hand, it's more efficient to do these small refactorings in parallel with the actual bugfixes, which are orthogonal and so will not be obfuscated by this refactor. Note `git diff` does properly track the file name change. However, it also makes a potentially more pervasive change: in order for the numerically-sensitive algorithm to work, it changes the condition on when two very close (double) numbers are "nearlyEqual" from differing by less than DBL_EPSILON to differing by less than or equal to DBL_EPSILON. Although this may change other behaviors than the ones primarily being addressed, I believe it is an acceptable change because (a) It preserves all tests. (b) DBL_EPSILON is well below the standard config.epsilon anyway (c) I believe there are extant issues noting the odd/inconsistent behavior of nearlyEqual near 0 anyway, so I believe this will be overhauled in the future in any case. If so, the eigenvector computation will make a good test that a future nearlyEqual algorithm is working well. To be clear, the direct motivation for the change is that there are multiple cases in the eigenvector computation in which a coefficient that is "supposed" to be zero comes out to precisely DBL_EPSILON, which is fairly unsurprising given that these coefficients are produced by subtracting an eigenvalue from a diagonal entry of a matrix, which is likely to be essentially equal to that eigenvalue. As many tests of defective matrices as I could readily find by web searching have been added as unit tests (and one more in the typescript type testing). An additional case I found still fails, but in the _eigenvalue_ computation rather than the _eigenvector_ search, so that was deemed beyond the scope of this PR and has been filed as issue josdejong#3036. Resolves josdejong#2879. Resolves josdejong#2927. Resolves josdejong#3014.
Fixed in |
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 eigenvalue2
.To Reproduce
Steps to reproduce the behavior.
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)
The text was updated successfully, but these errors were encountered: