eig to return eigenvalues in order of descending absolute value

[mhauru]

Eigenvalue decomposition of a positive-definite Hermitian matrix is a special case of SVD. However, svd returns the singular values in descending order, where as eig returns them in ascending order. The order of singular values is fixed popular convention, so I think for consistency eig should adhere to the same order.

Further more, I would prefer the ordering for eig to be descending absolute value. This would also make it consistent with the default behavior of eigs and with the fact that singular values of Hermitian matrices are the absolute values of the eigenvalues.

I don't know how this would affect perfomance, i.e. can LAPACK do this for us naturally or would we need to sort (and permute the eigenvectors) afterwards.

[andreasnoack]

We follow LAPACK here and I think it is unlikely that we want to deviate from them. In Julia, it's quite easy for the user to sort the values in the order they like and there is not a huge cost concern here since the sorting time is fairly small compared to the compute time for the eigenvalue decomposition anyway.

[nalimilan]

As a data point, R always uses a descending order, but both MATLAB and NumPy directly return what LAPACK computes. IIUC, LAPACK does not actually guarantee any order for eigenvalues, though it does for singular values.

Given that (as @andreasnoack noted) the cost of sorting is low compared with the computation, how about sorting them automatically? Seeing the number of questions about this on the forums, it's certainly confusing that the values are sorted most of the time in practice, but not always. Sorting could prevent a few mistakes.

[andreasnoack]

The values from the symmetric eigenvalue problem are always sorted but for the general eigenvalue problem they are not sorted. In general, they are complex so sorting opens a discussion about how to sort them.

I'm not sure there is a single sorting that will make all users happy. It is true that the positive definite eigenproblem is equivalent to the singular value problem but what is actually solved is the general Hermitian/symmetric eigenvalue problem which is more general since the eigenvalues can be negative. As @mhauru points out, we could just sort according to largest absolute value but I'm not sure that it generally is what users want. Regarding eigs, then I think the sorting depends on which values you ask for e.g. :SM, :LM.

[mhauru]

For eigs the default, however, is :LM, and I think there's a good case for that. The largest magnitude eigenvalues are in many ways the most important ones, the ones that tell you the most about your matrix. For instance they are the ones you want to keep when you approximate your matrix by a smaller rank matrix. Maybe I'm too biased by my own use case (I often just switch from svd to eig for performance reasons when my matrix is Hermitian), but for me it seems clear that descending magnitude is the most natural way to sort the eigenvalues.

I do admit that conforming to how LAPACK works has value in itself. Whether that outweighs the benefits of what I consider to be conceptually the right thing to do -- sorting by descending magnitude -- is something that I leave up to the judgement of more established Julia contributors.

If the sorting is implemented as an additional step after the decomposition, meaning there's a subleading increase in running time, I think it should be made an optional/keyword argument that defaults to "yes, do sort". eig is the bottleneck in a lot of algorithms and I wouldn't want to impose even subleading time costs on anybody without the possibility to opt out.

[nalimilan]

I would say that almost any sorting method is better than no guaranteed order at all (though it's logical to be consistent with eigs). Note that Julia is very different from LAPACK since a single function eigfact wraps several Fortran routines. It could be very confusing that the order of returned eigenvalues varies depending on the type of the input matrix.

In any case, the docs need to mention this. Currently, they implicitly assume the user knows which order is used, without defining it:

If A is Symmetric, Hermitian or SymTridiagonal, it is possible to calculate
only a subset of the eigenvalues by specifying either a UnitRange irange
covering indices of the sorted eigenvalues...

[jiahao]

For eigs the default, however, is :LM, and I think there's a good case for that. The largest magnitude eigenvalues are in many ways the most important ones, the ones that tell you the most about your matrix.

eigs and eig are completely different computations underneath though. Iterative methods such as the implicitly restarted Arnoldi method used by eigs have a hard time computing anything other than :LM; therefore, :LM is the default. The dense methods provided by eig/eigfact have no such restriction.

As @andreasnoack says, sorting after the fact is cheap, and there is no "correct" order, so this really is a documentation issue as @nalimilan describes rather than a problem in what we currently return.

[juliohm]

I am experiencing similar issue with eigvals.

Consider the following example:

A = rand(4,4)
eigvals(A)

4-element Array{Complex{Float64},1}:
2.28902+0.0im
-0.0557381+0.0im
-0.321139+0.207219im
-0.321139-0.207219im

These are ordered in descending magnitude.

Now, continue with:

B = A'*A
eigvals(B)

4-element Array{Float64,1}:
0.0014481
0.119291
0.370974
5.41973

These are ordered in ascending magnitude.

Could you confirm this is expected and that the behavior won't be changed?

[andreasnoack]

These are ordered in descending magnitude.

The eigenvalues of the non-symmetric case are not sorted.

These are ordered in ascending magnitude.

True. The eigenvalues of the symmetric case are sorted in ascending order.

Both behaviors are inherited from LAPACK.

Could you confirm this is expected and that the behavior won't be changed?

It's expected but I don't know if it will change. Personally, I think you should sort your eigenvalues in the way that fits your problem.

[juliohm]

Thank you @andreasnoack, I am enforcing magnitude order in my user code now.

[stevengj]

From a pedagogical perspective (e.g. when teaching linear algebra), the "random" ordering is a continual annoyance. Almost any consistent ordering would be better than what we have now, in my opinion, e.g. sorting lexicographically in descending order by (abs(λ),real(λ),imag(λ)).

Since re-ordering is virtually free compared to the cost of eigs or eigvals, why not? If the user wants a different order, they can always re-order again.

[StefanKarpinski]

💯 agree

[stevengj]

Closed by JuliaLang/julia#21598