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
Merged

Conversation

stevengj
Copy link
Member

@stevengj stevengj 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.

@stevengj stevengj added the domain:linear algebra Linear algebra label Apr 27, 2017
Copy link
Member

@andreasnoack andreasnoack left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would something like

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

be any more efficient, as it avoids the splatting?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@tkelman
Copy link
Contributor

tkelman 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
Copy link
Member Author

@tkelman, yes.

@tkelman
Copy link
Contributor

tkelman 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
Copy link
Member

andreasnoack 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
Copy link
Member Author

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

@nanosoldier
Copy link
Collaborator

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

@Sacha0
Copy link
Member

Sacha0 commented Apr 29, 2017

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

@nanosoldier
Copy link
Collaborator

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

@stevengj
Copy link
Member Author

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

@tkelman
Copy link
Contributor

tkelman 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
Copy link
Member

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
Copy link
Member

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
Copy link
Sponsor Member

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
Copy link
Member Author

stevengj 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
Copy link
Member Author

stevengj 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
Copy link
Sponsor Member

Seems good to me!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

8 participants