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.

[m93a]

Thanks for investigating the issue, Glen!
Upon reading this issue for the first time, I was quite frankly clueless. Now that you mention the the cause, I recall thinking the eigs API was quite cumbersome and counter-intuitive. That was the reason why I proposed the m.columns() method in #2177, which I originally intended to implement in PR #2155, but sadly never really finished because of a Babel issue.

However, my 2¢ are that changing the eigs API should at least be considered, as it confused even the project owner and the guy who implemented a part of it 😅

EDIT: Oh, now I recall even thinking about how to improve the API, and failing to. The problem is that in Math.js a two-dimensional array of numbers is automatically interpreted as a row matrix. So it is literally impossible to return an array of vectors and not have users confuse it with a row matrix. If there ever was Math.js 2, I would strongly argue against that decision, having now seen too many confusing edge-cases this leads to. But right now, it is baked deep inside Math.js, so there's no coming back.

EDIT 2: How about renaming the prop to vectorsAsColumns? That sounds about as specific as it can get.

[gwhitney]

So it is literally impossible to return an array of vectors and not have users confuse it with a row matrix

Right, I was trying to say that in my analysis, but may not have been clear enough. That's the central problem here.

I do think some changes are in order, and a PR will be needed for #2879 anyway, so there is a good opportunity here.
If we don't want this to be a breaking change, the easiest thing would be to add new clearer property names/values leaving the old ones alone for now but with them labeled as deprecated, and in the next breaking change retire the confusing properties.

Note that because of defective matrices, the eigenvector matrix will in generally not be nxn, and we do need to associate specific eigenvalues with specific eigenvectors. So one totally unambiguous data structure to pass back would be an array of objects, in order from largest eigenvalue (in abs value) to smallest, with each object having just two properties: value and vector (a single vector for each). The value property could repeat, but all of the vector properties would be orthogonal vectors. I think this would be the best thing to return, better than anything I suggested already in #2879, so I will link from there to here also. The only ingredient missing is what would the name of this property be? Maybe labeledVectors? Would be very glad for a better suggestion!

[josdejong]

Thanks for your clear explanation Glen, I now see that the results are indeed technically correct 😅, though the API and the actual values cause confusion.

So, we're mixing rows/columns when extracting the vector.

As a solution, I propose the following:

1. Let us add explanation in the documentation about this right now, that is just a few minutes of work. I'll do that today.
2. Making eigs output vectors that are easier on the eye requires a large rewrite if I understand it correctly. Or maybe we can add a minimalistic normalization step at the end or so. Anyway, it's not the biggest priority to me right now since eigs returns technically correct results. A PR improving it would be welcome of course.
3. We could spend time on an intermediate, non-breaking solution to improve the API of eigs, like adding an extra property vectorsAsColumns alongside the confusing vectors, but I prefer to skip this temporary solution and directly implement a better API in one go.
4. We indeed need to come up with a better API. I love the idea of Glen to change the API to an array with objects
   [{value, vector}, ...] for each pair, and then we have freedom to output the results for defective matrices in a meaningful way I think (see explanation by Glen). Any concrete proposal on an API for defective matrices? It will be a breaking change, but I'm convinced it is for the better.

Does that make sense?

[josdejong]

@gwhitney so I like your idea of an API [{value, vector}, ...] , but I'm not sure yet what cases there are regarding defective matrices. We need to cater for that too. Would that be like one eigenvalue having multiple vectors: [{value: ..., vectors: [...] }]?

[bartekleon]

@gwhitney so I like your idea of an API [{value, vector}, ...] , but I'm not sure yet what cases there are regarding defective matrices. We need to cater for that too. Would that be like one eigenvalue having multiple vectors: [{value: ..., vectors: [...] }]?

I think there are only two cases really:

One being the ideal one:
having n independent eigenvalues

Other being defective matrix:
having same eigenvalue (defective matrix) -- the eigenvectors can be different or the same.

I think there might be a case that you can't perform eigenvalue decomposition at all. (I do believe it might be a case of non-diagonalizable matrices. But teeechnically you still can do eigenvalue decomposition)

E.x
[0, 1]
[0, 0]

you have "two" eigenvalues being 0
you have only one eigenvector (1, 0) [for both cases]

But generally having eigenvalue = 0 means that matrix is not diagonalizable

[gwhitney]

I feel strongly that the data structure we pass back should not depend on whether the matrix is defective or not, so we should design a data structure that works in either case; the current one does not. Jos has already said he prefers a single breaking change to a non-breaking one followed by a breaking one. So I propose we pass back an array of values as we do now, with their algebraic multiplicities (so the vector will always have n entries for an nxn matrix), and then under the key 'eigenvectors' (just to change the key, so that anyone switching will immediately know there is a difference as there is no 'vectors' key), we pass back an array of plain objects, one for each independent eigenvector. Each of these plain objects has two keys: value and vector. The value property reiterates the eigenvalue for that eigenvector, but is necessary because you can't predict the geometric multiplicities from the arithmetic ones. The vector property gives a single eigenvector for that value.

An alternative might be a map from eigenvalues to eigenvectors, but then it seems as though the values would be arrays of vectors, which have the confusability property that started this whole thread.

I am not wedded to this proposal by any means, but I haven't thought of a better one.

[gwhitney]

Sorry I don't think I addressed the comment of @josdejong about how this API proposal handles defective matrices. For non-defective matrices, there is one item in the eigenvectors array of objects for each entry in the values array, e.g., if an eigenvalue occurs three times in the values list, there will be three entries in the eigenvectors array that have their value property equal to that eigenvalue. For defective matrices, there will be fewer entries in the eigenvectors array than in the values array; at least one of the eigenvalues will appear fewer times in the eigenvectors array than in the values array, since the number of occurrences in the values array is the algebraic multiplicity, but the number of occurrences in the eigenvectors array is the geometric multiplicity. In general, geometric multiplicity <= algebraic multiplicity. Hopefully I've made my proposal completely clear now. If not, feel free to ask questions, or propose alternative APIs. The key features are that the client of the function may want/need to know (A) what the eigenvalues are; (B) what the algebraic multiplicity and the geometric multiplicity of each is; and (C) for each eigenvalue, as many independent eigenvectors as the geometric multiplicity of that eigenvalue. Any API that can convey all of (A), (B), (C) in an unambiguous and easily-interpretable way is fine with me. The current proposal is just the best I've managed to come up with, partly because arrays of vectors just don't work too nicely in mathjs when what you want is a list of vectors, rather than a matrix...

[josdejong]

I feel strongly that the data structure we pass back should not depend on whether the matrix is defective or not

Yes that is a very good point, let's make sure of that.

The current API looks like this:

math.eigs([[2, 0, 0], [0, 1, 0], [0, 0, 5]])
// {
//   values: [1, 2, 5],
//   vectors: [  // Remember to take the columns, not the rows!
//     [ 0, 1, 0 ],
//     [ 1, 0, 0 ],
//     [ 0, 0, 1 ]
//   ]
// }

math.eigs([[2, 1], [0, 2]])
// Currently: Error
// Result should be: a double eigenvalue 2, with only one distinct eigenvector `[1; 0]`

One question: I'm not sure whether the index of the eigenvector currently correlates with the index of the correspondingeigenvalue, is that the case? I guess so?

@gwhitney If we output fewer eigenvectors than eigenvalues, how can we know when eigenvector belongs to which eigenvalue? I would like an API where that is clear. What do you have in mind in that regard?

Here some ideas with different output for math.eigs([[2, 1], [0, 2]]):

Idea 1

(I think too complex)

{
  values: [2, 2],
  vectors: [
    {
      vector: [1, 0],
      valueIndices: [0, 1]
    }
  ]
}

Idea 2

Just repeat double eigenvalues and eigenvectors, and inform about that via a flag

[
  {
    value: 2,
    vector: [1, 0],
    multiplicity: 2
  },
  {
    value: 2,
    vector: [1, 0],
    multiplicity: 2
  }
]

Idea 3

Do not repeat double eigenvalues, and inform about the multiplicity

[
  {
    value: 2,
    vector: [1, 0],
    multiplicity: 2
  }
]

Any thoughts on this? I'm not sure if it is handiest in practice to always have a number of eigenvalues equal to the size of your matrix, or not (that would be the only reason for Idea 2).

[bartekleon]

@josdejong you have have same eigenvalues having different eigenvectors so option 3 isn't really perfect.

Maybe something like:

Array[{
  value: vector
  vectors: vector[]
  multiplicity: number
}]

So:

[
  {
    value: 2,
    vectors: [[1, 0], [1, 0]],
    multiplicity: 2
  },
 {
    value: 3,
    vectors: [[1, 0], [0, 1]],
    multiplicity: 2
  },
  ​{​
    ​value​: ​1​,​
    ​vector​s: ​[[​1​,​ ​0​]​],​
    ​multiplicity​: ​1
  ​}​
]

Also general rule is to order eigenvalues by their real part (I believe)

[gwhitney]

Sorry if I have not been writing clearly. Altihough I understand that you are OK with breaking changes, I think that we should try to keep the interface as close as we can. In particular, I think we should continue to return as the values key a vector always of length n of the eigenvalues. That conveys two of the four pieces of information the client may want: (A) what the eigenvalues are and (B) the algebraic multiplicity of each. Also a significant amount of time all the client wants is the values, so they should be kept easy to access.

Then we need to replace the current vectors key, as it is ambiguous in defective cases, and prone to misinterpretation. The remaining info we need to convey is (C) the geometric multiplicity of each eigenvalue and (D) as many independent eigenvectors as exist. Since eigenvalues may repeat and we don't want to put eigenvectors in lists, a map from eigenvalues to eigenvectors is out. Instead we want a more general association of values and vectors, which we can pass back in a new key to make transition safer. I am proposing the key 'eigenvectors'. The value of this key would be a list of plain objects, each with a 'value' property giving the eigenvalue and a 'vector' property giving the eigenvector. An alternative would be just a list of pairs [value, vector] but that's more opaque. Anyhow, this gives you D from the vector part of each entry and C by counting how many times a given value appears.

To make this concrete, suppose M has eigenvalue 4 with algebraic multiplicity 3 but geometric multiplicity 2 and independent eigenvectors [1,0,0,0] and [0,1,0,0], and also eigenvalue -2 with algebraic multiplicity 1 and eigenvector [0,0,0,1]. Then I propose eigs would return:

{values: [4, 4, 4, -2],
 eigenvectors: [
    {value: 4, vector: [1,0,0,0]},
    {value: 4, vector: [0,1,0,0]},
    {value: -2, vector: [0,0,0,1]}
  ]
}

(but as I mention we could opt for the value of eigenvectors to be

[[4, [1,0,0,0]], [4, [0,1,0,0]], [-2, [0,0,0,1]]]

to be more compact but also more opaque.)

As to your ideas 1,2,3:
In idea 2, there is no benefit to repeating eigenvectors, so I would recommend against that. In idea 3, yes all of the info is there but we should rename 'multiplicity' to 'algebraicMultiplicity' to be clear. Then it has drawbacks that the algebraic multiplicity is repeated for each eigenvector with a given eigenvalue -- a bit redundant. Also I don't like the fact that it gets rid of the simple vector of n eigenvalues, a part of the current interface that is not broken. And idea 1 has some internal confusion; an eigenvector does not associate to two eigenvalues, not even when those are the same value.

As for bartekleon's suggestion, it has arrays of vectors which I think we are trying to avoid because of the structural confusion of array of vectors and matrix in mathjs, exacerbated by the fact that then the vectors look like rows when they are actually columns.

Hopefully I have finally been clearer. I am fine in the end with any return value that conveys all of A,B,C,D. Within that we should strive to maximize clarity and minimize redundancy. Thanks for considering.

[josdejong]

Thanks Glen and Bartosz.

From what I read here it looks like at least Matlab doesn't sort the returned eigenvalues. I'm not sure if there is value in sorting them. At least you can always sort them yourselves if needed. I'm OK with either sorting the output or not though if you guys think it is helpful.

I like Glen's proposal a lot:

{
  values: [4, 4, 4, -2],
  eigenvectors: [
    {value: 4, vector: [1,0,0,0]},
    {value: 4, vector: [0,1,0,0]},
    {value: -2, vector: [0,0,0,1]}
  ]
}

It solves the ambiguity that vectors currently has, it explicitly couples the value and vector that belong together. Providing an array with just the values makes sense, it's a common case to just need the eigenvalues I think. Also, I was wondering whether we should provide information about the multiplicity, but this is probably not necessary with the above interface: you can deduct the algebraic multiplicity by grouping values, and the geometric multiplicity by grouping eigenvectors by value. So all information is there.

One last thought: using names values and eigenvectors feels inconsistent, I would expect either values and vectors or eigenvalues and eigenvectors. People upgrading to this new interface have broken code in any case: either removing vectors (and adding a new eigenvalues), or changing vectors to an array with objects. Only in the first case we could throw an informative error when people access the old vectors property on the returned object, that would be very useful. Still, for the smoothest upgrade I think we should go with your proposal, and let the old property vectors throw an error, like:

{
  values: Array | Matrix
  eigenvectors: Array<{ value: number | BigNumber, vector: Array | Matrix}>
  vectors: void // throws an error explaining to use eigenvectors instead
}

What do you think?

[bartekleon]

Looks good to me :)

[gwhitney]

I do agree that the top level properties being values and eigenvectors is a bit unfortunate. But I don't see any reason to change values whereas I do feel like there is a need to change vectors to avoid obscuring the problems for any current clients when they change over if they don't pay attention to the CHANGES list in the release. But we could respond by changing both names to another pair that match better, even though it is more inconvenient for the (probably majority) of folks that use just the values. In that case, eigenvalues and eigenvectors makes the most sense to me, even though they are a bit verbose. Do you prefer that pair, or in the end do you think values/eigenvectors is best for the new interface?

And to your final point -- are you saying you want to, in the new version, return an object that has a getter/setter for the vectors property that throw errors pointing clients to the new eigenvectors property? That's a reasonable idea, I can make that so in the PR. (But then would we schedule that extra, non-useful stuff in the returned object to be removed say two rounds of breaking changes later, once "everyone" has made the transition?)

OK, I think we are super close --- I think Jos can just make his pronouncements on these final tiny points and then I can get started on a PR. Thanks to all!

[gwhitney]

I'm not sure if there is value in sorting them. At least you can always sort them yourselves if needed. I'm OK with either sorting the output or not though if you guys think it is helpful.

Ah, missed this point. I would not alter the order of the values property from whatever it is now. Even if you decide we should rename it, my goal for the PR is that there would be no change in what array of eigenvalues is produced. So if it's sorted now, it would stay sorted, and if it's not, it would stay that way, too. In my mind, the fewer collateral changes from bugfixes, the better.

[josdejong]

👍 sounds good. Ok let's do the following:

1. return an object with properties values and eigenvectors as discussed. It's indeed not ideal but I think the best compromise between improving the API and don't breaking everyone's existing code (when using values).
2. let's return a third property vectors, which is a getter that throws an error explaining to use eigenvectors instead.
3. Let's not add logic to sort the results.

[josdejong]

Fixed in v12.0.0 now via #3037