-
-
Notifications
You must be signed in to change notification settings - Fork 5.5k
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
Implement sorting for Eigen #24536
Implement sorting for Eigen #24536
Changes from 4 commits
a7340a0
468709f
e5aa710
5af6876
0594363
447ea8e
c65c90f
53b0be2
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -71,6 +71,35 @@ isperm(p::Tuple{}) = true | |
isperm(p::Tuple{Int}) = p[1] == 1 | ||
isperm(p::Tuple{Int,Int}) = ((p[1] == 1) & (p[2] == 2)) | ((p[1] == 2) & (p[2] == 1)) | ||
|
||
#swap columns i and j of a, in-place | ||
function swapcols!(a::AbstractMatrix, i, j) | ||
i == j && return | ||
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] | ||
end | ||
end | ||
# like permute!! applied to each row of a, in-place in a (overwriting p). | ||
function permutecols!!(a::AbstractMatrix, p::AbstractVector{<:Integer}) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. might be generally useful to have and export, along with corresponding |
||
count = 0 | ||
start = 0 | ||
while count < length(p) | ||
ptr = start = findnext(!iszero, p, start+1) | ||
next = p[start] | ||
count += 1 | ||
while next != start | ||
swapcols!(a, ptr, next) | ||
p[ptr] = 0 | ||
ptr = next | ||
next = p[next] | ||
count += 1 | ||
end | ||
p[ptr] = 0 | ||
end | ||
a | ||
end | ||
|
||
function permute!!(a, p::AbstractVector{<:Integer}) | ||
count = 0 | ||
start = 0 | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -27,6 +27,45 @@ function getindex(A::Union{Eigen,GeneralizedEigen}, d::Symbol) | |
throw(KeyError(d)) | ||
end | ||
|
||
eigsortby(λ::Real) = λ | ||
eigsortby(λ::Complex) = reim(λ) | ||
|
||
""" | ||
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)`. | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Why point to There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Unless you do something like There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Even then, it's not clear what this reference means: does this function support all keyword arguments accepted by There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. 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 |
||
|
||
# Examples | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Should the example perhaps be attached to the There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I think that makes sense. But for There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. 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 Also, why the reference to |
||
```jldoctest | ||
julia> F = eigfact([4.0 1.0; 0 1.0]) | ||
Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([4.0, 1.0], [1.0 -0.316228; 0.0 0.948683]) | ||
|
||
julia> F[:values] | ||
2-element Array{Float64,1}: | ||
4.0 | ||
1.0 | ||
|
||
julia> sort!(F) | ||
Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}([1.0, 4.0], [-0.316228 1.0; 0.948683 0.0]) | ||
|
||
julia> F[:values] | ||
2-element Array{Float64,1}: | ||
1.0 | ||
4.0 | ||
``` | ||
""" | ||
function sort!(F::Eigen; by=eigsortby, kw...) | ||
if !issorted(F[:values], by=by) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. need to pass |
||
perm = sortperm(F[:values]; by=by, kw...) | ||
permute!(F[:values], perm) | ||
Base.permutecols!!(F[:vectors], perm) | ||
end | ||
return F | ||
end | ||
|
||
sort(F::Eigen; kw...) = sort!(deepcopy(F), kw...) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Could be worthwhile to define a Base.copy(F::Eigen) = Eigen(copy(F.values), copy(F.vectors)) |
||
|
||
isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values) | ||
|
||
""" | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -105,6 +105,13 @@ end | |
end | ||
end | ||
|
||
let a = rand(10, 10) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Perhaps add a |
||
f = eigfact(a) | ||
sort!(f) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. probably |
||
@test a ≈ f[:vectors] * Diagonal(f[:values]) / f[:vectors] | ||
end | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Since the purpose is to sort, should we test that too? E.g. |
||
|
||
|
||
# test a matrix larger than 140-by-140 for #14174 | ||
let aa = rand(200, 200) | ||
for a in (aa, view(aa, 1:n, 1:n)) | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
would be clearer with spaces after the outer commas