Implement sorting for Eigen #24536

dalum · 2017-11-08T16:48:59Z

The helper functions were taken verbatim from @stevengj's PR: #21598 (with a single minor deprecation fix). This PR implements the (probably) least radical solution to sorting eigenfactorisations proposed in #21598, with sorting of eigenvalues being opt-in rather than opt-out, and not enforcing any particular canonical sorting.

stevengj · 2017-11-08T17:19:17Z

base/linalg/eigen.jl

+```
+"""
+function sort!(F::Eigen; kw...)
+ perm = sortperm(F[:values]; kw...)


Seems like we need a default by=eigsortby, if by wasn't supplied in kw, using e.g. the eigsortby function from #21598.

Otherwise it will fail with no method matching isless for complex eigenvalues.

This was actually intentional, but maybe the error message is confusing? I wanted to avoid choosing any particular sorting order for complex eigenvalues, since I generally need by=angle and others might want real, imag or eigsortby.

I can also see a good argument for having some default, and eigsortby is a pretty good candidate.

This could also be done by dispatching on the eltype of the eigenvalues, defaulting to reim, without the need to define eigsortby.

Definitely the caller should be able to supply their own by, but I think there should be a default, and reim seems like a good a choice for complex values (since it will match the sorting of real values in the case of complex numbers that are purely real).

stevengj · 2017-11-08T18:04:49Z

base/linalg/eigen.jl

- permute!(F[:values], perm)
- Base.permutecols!!(F[:vectors], perm)
+function sort!(F::Eigen; by=eigsortby, kw...)
+ if !issorted(F[:values], by=by)


need to pass kw... to issorted too.

stevengj · 2017-11-08T18:08:14Z

test/linalg/eigen.jl

@@ -105,6 +105,13 @@ end
 end
 end

+let a = rand(10, 10)
+ f = eigfact(a)
+ sort!(f)


probably f = sort(eigfact(a)) instead to make sure that both sort and sort! are tested.

fredrikekre · 2017-11-09T12:21:03Z

base/linalg/eigen.jl

+ return F
+end
+
+sort(F::Eigen; kw...) = sort!(deepcopy(F), kw...)


Could be worthwhile to define a copy method for Eigen. Such methods exists for e.g. Cholesky and LU factorizations. Something like:

Base.copy(F::Eigen) = Eigen(copy(F.values), copy(F.vectors))

fredrikekre · 2017-11-09T12:23:36Z

base/linalg/eigen.jl

+Sort the eigenvectors and eigenvalues of an eigenfactorization `F` in place using the eigenvalues
+for comparison. See also `sort!(::AbstractVector)`.
+
+# Examples


Should the example perhaps be attached to the sort method instead? It is a fairly common pattern throughout the documentation to have a more extensive docstring for the non-inplace method, and just reference that there is an inplace equivalent.

I think that makes sense. But for sort(::Vector), it actually references sort! as the extensive docstring. Maybe this should be changed?

Hmm okay, I guess we don't have an official policy, in any case it is unnecessary to repeat information, so one or the other should have a more extensive docstring. Perhaps you can simply add a short docstring for the non-inplace version instead then, and reference this. I guess most of the time you want to use sort! anyway.

Also, why the reference to sort!(::AbstractVector)? For docs about the keywords?

fredrikekre · 2017-11-09T12:25:48Z

test/linalg/eigen.jl

@@ -105,6 +105,12 @@ end
 end
 end

+let a = rand(10, 10)


Perhaps add a b = a'a and test the case with real eigenvalues too?

fredrikekre · 2017-11-09T12:26:26Z

test/linalg/eigen.jl

+let a = rand(10, 10)
+ f = sort(eigfact(a))
+ @test a ≈ f[:vectors] * Diagonal(f[:values]) / f[:vectors]
+end


Since the purpose is to sort, should we test that too? E.g. @test issorted(...)

fredrikekre · 2017-11-09T12:31:31Z

Should probably also add sort[!](::Eigen) to appropriate manual section, otherwise it will show up at some random place with the rest of the sort[!] methods.

nalimilan · 2017-11-09T12:53:09Z

base/linalg/eigen.jl

+ sort!(F::Base.LinAlg.Eigen; kw...)
+
+Sort the eigenvectors and eigenvalues of an eigenfactorization `F` in place using the eigenvalues
+for comparison. See also `sort!(::AbstractVector)`.


Why point to sort!(::AbstractVector)? Other docstrings for sort! will be listed in the same place anyway.

Unless you do something like ?> sort!(F) to filter out all the noise. 🙂 Or in the manual.

Even then, it's not clear what this reference means: does this function support all keyword arguments accepted by sort!(::AbstractVector)? Is it just that they are similar? Also, the generic docstring is defined for any v, not just ::AbstractVector.

You are right. It is because it accepts the same keyword arguments. I will change the wording to make it more clear. Also, AFAIU, the "generic" docstring is a fallback to all docstrings, but the one currently being shown is only for AbstractVector, despite the missing ::AbstractVector in the signature of the docstring.

StefanKarpinski

Overall, this is great to have a nice interface for, but I'm not convinced that adding this as a method to sort and sort! for the Eigen type is the right way to go. Those functions are for sorting collections of values, which I'm not sure we can legitimately consider factorization objects to be. However, permuting the rows/columns of factorizations according to some sort order seems to be a more general thing, so perhaps it would make sense to have sortfact and sortfact! functions?

StefanKarpinski · 2017-11-09T15:23:51Z

base/combinatorics.jl

+ cols = indices(a,2)
+ (i in cols && j in cols) || throw(BoundsError())
+ for k in indices(a,1)
+ @inbounds a[k,i],a[k,j] = a[k,j],a[k,i]


would be clearer with spaces after the outer commas

StefanKarpinski · 2017-11-09T15:26:43Z

base/combinatorics.jl

+ end
+end
+# like permute!! applied to each row of a, in-place in a (overwriting p).
+function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer})


might be generally useful to have and export, along with corresponding permutecols!, permutecols, permuterows!!, permuterows! and permuterows (yikes, that's a lot variations).

nalimilan · 2017-11-09T17:14:28Z

How about adding a sort::Bool keyword argument to the corresponding factorization functions? That sounds simpler than adding a new function which will be harder to discover.

StefanKarpinski · 2017-11-09T17:22:28Z

Or sort::Union{Bool,Function}=false where if you ask for true you get a "standard" sort order, false is unsorted and if you pass a function, that's the function whose value is used to sort the eigenvalues (in this case). Is it always that eigenvalues that are used for sorting?

stevengj · 2017-11-09T17:39:48Z

See @andreasnoack's comment about a sort keyword for eig. But yes, I think is is always eigenvalues that are used for sorting. Any more complicated filtering based on the eigevectors is something you'd want to implement on your own as a postprocessing step.

Evey Dee added 3 commits November 8, 2017 17:12

Implement sort!(::Eigen)

a7340a0

Add tests

468709f

Remove redundant Base.LinAlg.

e5aa710

stevengj reviewed Nov 8, 2017

View reviewed changes

Let default sorting be eigsortby

5af6876

stevengj reviewed Nov 8, 2017

View reviewed changes

kshyatt added the domain:linear algebra Linear algebra label Nov 8, 2017

Add kw..., call sort in test

0594363

fredrikekre reviewed Nov 9, 2017

View reviewed changes

nalimilan reviewed Nov 9, 2017

View reviewed changes

Evey Dee added 3 commits November 9, 2017 14:22

Define copy and refine docstrings for Eigen

447ea8e

Add (@ref)s and manual entries

c65c90f

More tests

53b0be2

StefanKarpinski reviewed Nov 9, 2017

View reviewed changes

dalum closed this Feb 25, 2019

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement sorting for Eigen #24536

Implement sorting for Eigen #24536

dalum commented Nov 8, 2017

stevengj Nov 8, 2017 •

edited

Loading

dalum Nov 8, 2017 •

edited

Loading

dalum Nov 8, 2017 •

edited

Loading

dalum Nov 8, 2017

stevengj Nov 8, 2017

stevengj Nov 8, 2017

stevengj Nov 8, 2017

fredrikekre Nov 9, 2017

fredrikekre Nov 9, 2017

dalum Nov 9, 2017

fredrikekre Nov 9, 2017

fredrikekre Nov 9, 2017

fredrikekre Nov 9, 2017

fredrikekre commented Nov 9, 2017

nalimilan Nov 9, 2017

dalum Nov 9, 2017

nalimilan Nov 9, 2017

dalum Nov 9, 2017 •

edited

Loading

StefanKarpinski left a comment

StefanKarpinski Nov 9, 2017

StefanKarpinski Nov 9, 2017

nalimilan commented Nov 9, 2017

StefanKarpinski commented Nov 9, 2017

stevengj commented Nov 9, 2017

Implement sorting for Eigen #24536

Implement sorting for Eigen #24536

Conversation

dalum commented Nov 8, 2017

stevengj Nov 8, 2017 • edited Loading

Choose a reason for hiding this comment

dalum Nov 8, 2017 • edited Loading

Choose a reason for hiding this comment

dalum Nov 8, 2017 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fredrikekre commented Nov 9, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dalum Nov 9, 2017 • edited Loading

Choose a reason for hiding this comment

StefanKarpinski left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nalimilan commented Nov 9, 2017

StefanKarpinski commented Nov 9, 2017

stevengj commented Nov 9, 2017

stevengj Nov 8, 2017 •

edited

Loading

dalum Nov 8, 2017 •

edited

Loading

dalum Nov 8, 2017 •

edited

Loading

dalum Nov 9, 2017 •

edited

Loading