RFC: sort eigenvalues in a canonical order

[stevengj]

As discussed in JuliaLang/LinearAlgebra.jl#324, this PR sorts the eigenvalues (and corresponding eigenvectors) in a canonical order: in ascending order for real eigenvalues (consistent with what LAPACK does automatically in the Hermitian case, for which no sorting is required), and in ascending order lexicographically by (real(λ),imag(λ)) for complex eigenvalues (chosen so that it is consistent with the real case when the eigenvalues are real).

The basic principle here is that virtually any consistent, comprehensible order is preferable to random ordering, and the cost of the sorting is negligible for matrices of any significant size. (I measured a 15% overhead for a 10x10 complex matrix, and for bigger matrices it quickly becomes unobservable.)

Along the way, I noticed that eigvals(::Diagonal{Matrix}) seemed wrong: it returned an array of arrays of eigenvalues, rather than an array of eigenvalues (which are still well-defined scalars for block matrices, not vectors). I changed it, and sorted it to be consistent with other eigvals routines. Moved to JuliaLang/LinearAlgebra.jl#594

(I implemented an internal permutecols!! function to permute the columns of a matrix in-place with a given permutation vector. I tried the alternative of permuting the rows one by-one using permute!! on a view, but it was many times slower.)

See also JuliaLang/LinearAlgebra.jl#343, #9655, #8441.

Updates:

• As per the comments below, I added a sortby keyword to the various eig functions. If you pass sortby=nothing, it suppresses sorting. If you pass another by function, e.g. sortby=abs, that changes the sorting accordingly.

• I also changed the documentation to say that sortby may not be a keyword for special matrix types like SymTridiagonal that may implement their own sorting convention. SymTridiagonal already didn't accept the other eigen keywords like permute and scale, and allowing special matrices to decide on their preferred eigenvalue ordering seemed like the most sensible and conservative approach to me.

[andreasnoack]

But only after 0.6 has been branched.

[ararslan]

Would something like

sort!(mapreduce(eigvals, vcat, D.diag), by=eigsortby)

be any more efficient, as it avoids the splatting?

[stevengj]

That would be worse, I think, because it would call vcat pairwise over and over and generate lots of temporary arrays.

[ararslan]

Ah I forgot about the pairwise thing. Sorry for the noise.

[tkelman]

will there be a fallback for eigvecs that does the right thing?

[stevengj]

yes, the fallback calls eigfact(M)[:vectors], which gets the sorted eigenvectors.

[tkelman]

Is calling lapack functions directly (which won't work for all array types) the only way here that an expert user could avoid the sorting overhead if they had an application that involved many small matrices and didn't need an ordering?

[stevengj]

@tkelman, yes.

[tkelman]

That's an unfortunate loss of generality then. If wanting an ordered output is primarily for pedagogical purposes, then I think it would be worth having a more convenient way of opting out of the sorting overhead. Other than the dense lapack types, do our other implementations also return unsorted results right now? If not, the sorting should maybe only be put in the code path for those dense types, rather than in the generic paths.

[andreasnoack]

That's an unfortunate loss of generality then

I have been wondering the same thing but I've come to the conclusion that the generality isn't in demand here. As already mentioned, the symmetric solver in LAPACK sorts, see line 536 and I'm not aware of any complaints about that. Also, my guess is that the general eigenvalue problem is mainly for teaching and exploration. People who are really interested in performance and accuracy would/should use the Schur form anyway.

[stevengj]

@nanosoldier runbenchmarks(ALL, vs = ":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

[Sacha0]

@nanosoldier runbenchmarks("linalg" || "problem", vs = ":master")

[nanosoldier]

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

[stevengj]

Maybe we shouldn't sort eigenvalues for the diagonal, bidiagonal, and triangular special cases where they can be just read off the diagonal?

[tkelman]

That's making my concern about generic code even worse. If we sort, we should do it predictably as part of the documented contract of what the generic eig* functions should do for all types, including user defined ones to be a correct implementation. Maybe a kwarg as a way of opting out of sorting?

[andreasnoack]

I'd like to avoid adding a keyword argument if possible. We already have a generic function for sorting and a type for the eigendecomposition so would it make sense to overload sort!(::Eigen) for this? Retrieving the sorting permutation and applying it to the values and vectors would be a bit annoying but wrapping the call in a sort! is simple and transparent. We should then document that, in general, eigenvalues and vectors are not sorted but that it is easy and cheap to do with sort!. It doesn't give you sorted values by default but I think @stevengj's comment about triangular matrices suggests that adding sorting to the definition might not be the best thing even pedagogically. At least I think it is a useful point that any permutation of the values and vectors represents an eigendecomposition.

Another thing is that it might be a bit weird to error on sort(complex.(randn(3), randn(3))) while happily sorting complex eigenvalues by default.

[nalimilan]

Adding a keyword argument isn't a big deal, and it's easier to discover than sort!. It wouldn't take a lot a space in the docstring since we need to mention the default ordering strategy anyway.

[StefanKarpinski]

Sorting by default with a keyword to opt out in the rare case that it is too expensive to sort seems like a good balance to me.

[stevengj]

I updated this PR for master. As per the comments above, I added a sortby keyword to the various eig functions. If you pass sortby=nothing, it suppresses sorting. If you pass another by function, e.g. sortby=abs, that changes the sorting accordingly.

I also changed the documentation to say that sortby may not be a keyword for special matrix types like SymTridiagonal that may implement their own sorting convention. SymTridiagonal already didn't accept the other eigen keywords like permute and scale, and allowing special matrices to decide on their preferred eigenvalue ordering seemed like the most sensible and conservative approach to me.

[stevengj]

This seems to be working now. (Travis CI failure on MacOS is #30949.)

@andreasnoack, @StefanKarpinski, can we consider merging this?

[StefanKarpinski]

Seems good to me!