make eigvals(::Diagonal) accept sortby

[araujoms]

Solves the Diagonal part of #1471 and #1450.

I thought it would be better to start with Diagonal, as it was the only part that required a bit of thinking. I think it also makes sense in itself, since eigen accepted sortby, just eigvals and eigvecs did not. Even though the docstring said none did.

If this is welcome I'll also do SymTridiagonal, and then we can document that eigenvalues are sorted by default.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 93.92%. Comparing base (7e11b5e) to head (3b0cb4f).
⚠️ Report is 11 commits behind head on master.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1477   +/-   ##
=======================================
  Coverage   93.92%   93.92%           
=======================================
  Files          34       34           
  Lines       15933    15934    +1     
=======================================
+ Hits        14965    14966    +1     
  Misses        968      968           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ViralBShah]

cc @stevengj

[dkarrasch]

Just a minor comment. I like it in terms of consistency, but I'm a little afraid due to changing the default behavior.

[araujoms]

I think one can hardly rely on eigenvalues not being sorted, making them sorted by default is harmless.

Another possibility would be to add the sortby keyword argument to eigvals(::Diagonal) with default value nothing. It would make the interface inconsistent, but at least it would be possible to pass eigsortby and be certain that the eigenvalues will be sorted.

[dkarrasch]

Do we have a test for the eigenvectors with sorted eigenvalues? Currently, the default eigvecs is a Diagonal, so we can't sort inplace, and need to make sure this works correctly.

About the default sorting, @stevengj's comment is relevant here: #1450 (comment) Before v.12, the eigenvectors were returned as a dense matrix, for which it doesn't seem to make any difference whether the eigenvectors are the identity matrix or a permuted version thereof. With the new Diagonal return type, it starts to make a difference, so, as you say, we may consider making the default for Diagonals no sorting, even though that would be inconsistent with the symmetric case. On the other hand, it would save memory and computational effort by default.

In any case, making sorting accessible is a clear (consistency) win!

[araujoms]

There is already a test, line 422 of test/diagonal.jl.

Note that I didn't change the behaviour of eigen(::Diagonal), just the default, so it will still return a Diagonal matrix of eigenvectors. I don't think that's a good idea, I think it should always return a dense matrix even when the matrix of eigenvectors is identity, for two reasons: 1 - it made the function type unstable, and 2 - the optimization is not really useful, as nobody calls eigen specifically on a diagonal matrix, this is done only in generic code. and for generic code it's better to avoid special cases.

[dkarrasch]

Now I see clearer. If you change the default, it won't return a Diagonal anymore, since that is only if the sortby function is nothing. In all other cases, it will allocate a matrix and potentially permute the ones. Type instability is not an issue as long as the sortby function is a constant, like hard-coded or not given at all and then filled by the default value. Given the older discussions, I'd personally prefer to keep this "less controversial" and keep the default, at least for now, but make eigvals accept the sorting function as the title says.

[araujoms]

Ah, I see, functions are types, so it is not actually type unstable. Just a trip hazard.

Is it less controversial, though? People already complained about it in the original PR JuliaLang/julia#21598, and another person besides me filed an issue about it.

[dkarrasch]

I guess differing opinions is what makes it controversial. 😂 I don't have a personal preference TBH.

If we change the default behavior, this needs to be listed in some sort of NEWS.md, together with the advice that if the user wants to benefit from low-memory eigvecs, the sortby=nothing keyword argument needs to be provided explicitly.

[araujoms]

Sorry, I thought you meant keeping the default sortby=nothing is less controversial than changing it to sortby=eigsortby. Now I see you meant that that adding sortby to eigvals(::Diagonal) is less controversial than changing the default.

Fair enough, I want to do another PR for SymTridiagonal anyway, I can reverse the change here and leave it for there.

[dkarrasch]

Yes, adding the keyword is nice, and an "uncontroversial" new feature if we, for the time being, use/leave the current default sorting function nothing. Adding the keyword also to SymTridiagonal is also nice.

Orthogonally to that, perhaps in an issue, we should discuss if we shift towards having all eig*, which accept a sorting function, towards using eigsortby. But I guess we should try to reach more people and get them involved.

[araujoms]

I did open an issue to discuss that before writing the PR, that was #1471, which was closed as a duplicate of #1450, which got crickets.