eigs function b̶r̶o̶k̶e̶n̶ has a confusing interface

[been-jamming]

math.eigs produces correct eigenvalues with incorrect eigenvectors in some cases. Note that the example in the documentation still seems to give the correct results.

To Reproduce
Execute the following:

T = [[1, 2], [4, 3]];
ans = math.eigs(T);
ans.vectors[0];

Result:

Array [ -0.9701425001453324, -2.4253562503633317 ]

Expected:

Array [-2, 1]

Edit: It looks like the correct eigenvector is actually [ -1, 1 ]

[josdejong]

Hm, that looks wrong indeed. Thanks for reporting.

Anyone able to dive into this bug with eigs? Help would be very welcome.

[bartekleon]

I can provide some code, but I don't know how general it is. It would need to be ported to math.js syntax. @josdejong do you want me to make PR or post it somewhere here so someone could take it on?

[josdejong]

Thanks! Both will help, though I guess a PR will be closest to a concrete fix.

[bartekleon]

Created a draft, so someone can get inspired :)

[bartekleon]

@been-jamming is it [-2, 1] for sure? wolfram alpha shows different results
https://www.wolframalpha.com/input?i=eigs%28%5B%5B1%2C+2%5D%2C+%5B4%2C+3%5D%5D%29

(I just want to be sure)

[been-jamming]

@bartekleon You're right, [-2, 1] is the eigenvector for the transpose. It looks like the correct eigenvector would be [-1, 1] instead.

[josdejong]

Thanks for the PRs (#3016, #3017) and your feedbacks @bartekleon and @been-jamming. I indeed have the feeling too that there is some more generic issue in play, and that we should figure that out rather than implementing a special case for a 2x2 matrix.

I did try out how math.eigs([[1, 2], [4, 3]]) evaluates in previous versions of mathjs. According to Wolfram Alpha, the correct result is an eigenvalue 5 with eigenvector [1,2], and an eigenvalue -1 with eigenvector [-1,1].

1. mathjs@6.6.0, PR Adding eigs function  #1705. The result is an error:

   Input matrix is not symmetric
   

2. mathjs@9.4.0, PR Extend eigs to a general matrix #1741. The result is an error:

   Failed to find eigenvectors for the following eigenvalues: -1
   

3. mathjs@9.4.2, PR Implemented iterative eigenvalue finder for eigs #2237. The result is (round-off errors removed):

   {
     "values": [-1, 5], 
     "vectors": [[0, 1], [1, 0]]
   }
   

4. mathjs@10.3.0, PR fix(eigs): Handle matrices with complex entries #2445. Should be non-related but I just wanted to verify. The result still is:

   {
     "values": [-1, 5], 
     "vectors": [[0, 1], [1, 0]]
   }
   

5. mathjs@10.4.1, PR fix(eigs): Track transforms in generating complex eigenvectors #2496. The result is:

   {
     "values": [-1, 5],
     "vectors": [[-0.9701425001453323, 2.4253562503633317], [0.9701425001453324, 4.850712500726662]]
   }
   

Maybe this gives a clue? @m93a or @gwhitney any idea?

[bartekleon]

I will try making more generalized version this week. If possible, test cases would be nice. (I don't know how big / many / edge cases I should make so contributions are welcome)

[josdejong]

Thanks Bartosz!

[gwhitney]

Sorry I didn't have a chance to look at this issue before now. I do not think there is any serious difficulty with the current behavior of eigs in mathjs, just some misunderstandings:

(1) eigs(A) by design returns the matrix E whose columns are eigenvectors. That's the canonical way they are presented, because then AE equals the matrix ΛE, where Λ = diag(λ) is the diagonal matrix of eigenvalues. However if in mathjs you let E = eigs(A).vectors, then E[0] is the first row of E (because that's how mathjs indexes matrices) It doesn't make sense to put eigenvectors in rows, because they are naturally multiplied by by A as columns, not rows. But there's also no practical way to alter the way mathjs indexes matrices (everything else would get more confusing). And it seems unwise to pass back an array of column vectors, because arrays of vectors are pretty much the same thing as matrices, and so it would be too easy to confuse that array of column vectors as a matrix, and hence you'd be interpreting columns as rows, which would be wrong.

So the quick upshot is that to get the eigenvector corresponding to the first eigenvalue, you want column(E,0) not E(0) = row(E(1))

Now I do agree this is nevertheless somewhat confusing, even though I hope I've justified that it's really the only reasonable way for the function to work. So a couple of things we might do:

(a) rename the property of eigs(A) from vectors (which may be inviting been-jamming's interpretation of that property) to something like evMatrix that at least makes it clear the result is a matrix, inviting the question of which way do the eigenvectors go? and/or
(b) put some prominent comments in the appropriate spots of the documentation that the way to get the nth eigenvector of A is column(eigs(A).vectors, n) not eigs(A).vectors[n].

(2) Then it was hard to recognize the transposition that was going on because of the "funny" numbers in the returned E. After all, we know that "the" eigenvectors of A are (1,2) and (-1, 1). But the fact is there is no such thing as "the" eigenvectors: if v is an eigenvector, then so is kv for any scalar k. Hence (1/3, 2/3) and (1,-1) are just as good a pair of eigenvectors, as are 2.42535(1,2) and -0.97014(-1,1) -- which just happens to be the pair of eigenvectors that this particular algorithm finds.

I do not think it's worth significant effort to replace mathjs's current eigs function, given its maturity (been through a lot of overhauls above!) unless someone in this conversation has a proposal that will be demonstrably faster/more accurate/better in some other metric. But I think a PR for either of (a) or (b) would be a good idea -- renaming the issue to reflect that.