`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.