Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RFC: sort eigenvalues in a canonical order #21598

Merged
merged 15 commits into from Feb 5, 2019

Conversation

@stevengj
Copy link
Member

commented Apr 27, 2017

As discussed in #15734, 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 #30681

(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 #17093, #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
Copy link
Member

left a comment

But only after 0.6 has been branched.

eigvals(D::Diagonal{<:Number}) = D.diag
eigvals(D::Diagonal) = [eigvals(x) for x in D.diag] #For block matrices, etc.
eigvals(D::Diagonal{<:Number}) = sort(D.diag, by=eigsortby)
eigvals(D::Diagonal) = sort!(vcat((eigvals(x) for x in D.diag)...),by=eigsortby) #For block matrices, etc.

This comment has been minimized.

Copy link
@ararslan

ararslan Apr 27, 2017

Member

Would something like

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

be any more efficient, as it avoids the splatting?

This comment has been minimized.

Copy link
@stevengj

stevengj Apr 28, 2017

Author Member

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

This comment has been minimized.

Copy link
@ararslan

ararslan Apr 28, 2017

Member

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

eigvals(M::Bidiagonal) = M.dv
function eigvecs(M::Bidiagonal{T}) where T
eigvals(M::Bidiagonal) = sort(M.dv, by=eigsortby)
function eigvecs_(M::Bidiagonal{T}) where T # non-sorted eigenvectors

This comment has been minimized.

Copy link
@tkelman

tkelman Apr 27, 2017

Contributor

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

This comment has been minimized.

Copy link
@stevengj

stevengj Apr 28, 2017

Author Member

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

@tkelman

This comment has been minimized.

Copy link
Contributor

commented Apr 28, 2017

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

This comment has been minimized.

Copy link
Member Author

commented Apr 28, 2017

@tkelman, yes.

@tkelman

This comment has been minimized.

Copy link
Contributor

commented Apr 28, 2017

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

This comment has been minimized.

Copy link
Member

commented Apr 28, 2017

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

This comment has been minimized.

Copy link
Member Author

commented Apr 28, 2017

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

@nanosoldier

This comment has been minimized.

Copy link
Collaborator

commented Apr 28, 2017

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

@Sacha0

This comment has been minimized.

Copy link
Member

commented Apr 29, 2017

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

@nanosoldier

This comment has been minimized.

Copy link
Collaborator

commented Apr 29, 2017

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

@stevengj

This comment has been minimized.

Copy link
Member Author

commented Apr 29, 2017

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

@tkelman

This comment has been minimized.

Copy link
Contributor

commented Apr 29, 2017

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

This comment has been minimized.

Copy link
Member

commented May 3, 2017

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

This comment has been minimized.

Copy link
Contributor

commented May 3, 2017

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

This comment has been minimized.

Copy link
Member

commented May 3, 2017

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.

@dalum dalum referenced this pull request Nov 8, 2017

@stevengj stevengj force-pushed the stevengj:eigsort branch from 4bce535 to 726b309 Jan 10, 2019

@stevengj

This comment has been minimized.

Copy link
Member Author

commented Jan 10, 2019

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 added 7 commits Jan 10, 2019
stevengj added 5 commits Feb 1, 2019
@stevengj

This comment has been minimized.

Copy link
Member Author

commented Feb 4, 2019

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

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

@StefanKarpinski

This comment has been minimized.

Copy link
Member

commented Feb 4, 2019

Seems good to me!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
8 participants
You can’t perform that action at this time.