RFC: sort eigenvalues in a canonical order (#21598)

* sort eigenvalues in a canonical order

NEWS.md

@@ -37,6 +37,7 @@ Standard library changes
 
 * Added keyword arguments `rtol`, `atol` to `pinv` and `nullspace` ([#29998]).
 * `UniformScaling` instances are now callable such that e.g. `I(3)` will produce a `Diagonal` matrix ([#30298]).
+* Eigenvalues λ of general matrices are now sorted lexicographically by (Re λ, Im λ) ([#21598]).
 
 #### SparseArrays
 

base/combinatorics.jl

@@ -63,6 +63,37 @@ 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 = axes(a,2)
+    @boundscheck i in cols || throw(BoundsError(a, (:,i)))
+    @boundscheck j in cols || throw(BoundsError(a, (:,j)))
+    for k in axes(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})
+    require_one_based_indexing(a, p)
+    count = 0
+    start = 0
+    while count < length(p)
+        ptr = start = findnext(!iszero, p, start+1)::Int
+        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})
     require_one_based_indexing(a, p)
     count = 0

stdlib/LinearAlgebra/docs/src/index.md

@@ -38,13 +38,13 @@ julia> A = [-4. -17.; 2. 2.]
 
 julia> eigvals(A)
 2-element Array{Complex{Float64},1}:
- -1.0 + 5.0im
  -1.0 - 5.0im
+ -1.0 + 5.0im
 
 julia> eigvecs(A)
 2×2 Array{Complex{Float64},2}:
-  0.945905+0.0im        0.945905-0.0im
- -0.166924-0.278207im  -0.166924+0.278207im
+  0.945905-0.0im        0.945905+0.0im
+ -0.166924+0.278207im  -0.166924-0.278207im
 ```
 
 In addition, Julia provides many [factorizations](@ref man-linalg-factorizations) which can be used to

stdlib/LinearAlgebra/src/diagonal.jl

@@ -528,11 +528,11 @@ eigvals(D::Diagonal{<:Number}; permute::Bool=true, scale::Bool=true) = D.diag
 eigvals(D::Diagonal; permute::Bool=true, scale::Bool=true) =
     [eigvals(x) for x in D.diag] #For block matrices, etc.
 eigvecs(D::Diagonal) = Matrix{eltype(D)}(I, size(D))
-function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true)
+function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=nothing)
     if any(!isfinite, D.diag)
         throw(ArgumentError("matrix contains Infs or NaNs"))
     end
-    Eigen(eigvals(D), eigvecs(D))
+    Eigen(sorteig!(eigvals(D), eigvecs(D), sortby)...)
 end
 
 #Singular system

stdlib/LinearAlgebra/src/eigen.jl

@@ -27,18 +27,32 @@ Base.iterate(S::Union{Eigen,GeneralizedEigen}, ::Val{:done}) = nothing
 
 isposdef(A::Union{Eigen,GeneralizedEigen}) = isreal(A.values) && all(x -> x > 0, A.values)
 
+# pick a canonical ordering to avoid returning eigenvalues in "random" order
+# as is the LAPACK default (for complex λ — LAPACK sorts by λ for the Hermitian/Symmetric case)
+eigsortby(λ::Real) = λ
+eigsortby(λ::Complex) = (real(λ),imag(λ))
+function sorteig!(λ::AbstractVector, X::AbstractMatrix, sortby::Union{Function,Nothing}=eigsortby)
+    if sortby !== nothing && !issorted(λ, by=sortby)
+        p = sortperm(λ; alg=QuickSort, by=sortby)
+        permute!(λ, p)
+        Base.permutecols!!(X, p)
+    end
+    return λ, X
+end
+sorteig!(λ::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby)
+
 """
-    eigen!(A, [B])
+    eigen!(A, [B]; permute, scale, sortby)
 
 Same as [`eigen`](@ref), but saves space by overwriting the input `A` (and
 `B`), instead of creating a copy.
 """
-function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasReal
+function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
     n = size(A, 2)
     n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
-    issymmetric(A) && return eigen!(Symmetric(A))
+    issymmetric(A) && return eigen!(Symmetric(A), sortby=sortby)
     A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
-    iszero(WI) && return Eigen(WR, VR)
+    iszero(WI) && return Eigen(sorteig!(WR, VR, sortby)...)
     evec = zeros(Complex{T}, n, n)
     j = 1
     while j <= n
@@ -53,18 +67,19 @@ function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where
         end
         j += 1
     end
-    return Eigen(complex.(WR, WI), evec)
+    return Eigen(sorteig!(complex.(WR, WI), evec, sortby)...)
 end
 
-function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasComplex
+function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
     n = size(A, 2)
     n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
-    ishermitian(A) && return eigen!(Hermitian(A))
-    return Eigen(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]...)
+    ishermitian(A) && return eigen!(Hermitian(A), sortby=sortby)
+    eval, evec = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]
+    return Eigen(sorteig!(eval, evec, sortby)...)
 end
 
 """
-    eigen(A; permute::Bool=true, scale::Bool=true) -> Eigen
+    eigen(A; permute::Bool=true, scale::Bool=true, sortby) -> Eigen
 
 Computes the eigenvalue decomposition of `A`, returning an `Eigen` factorization object `F`
 which contains the eigenvalues in `F.values` and the eigenvectors in the columns of the
@@ -79,6 +94,12 @@ before the eigenvector calculation. The option `permute=true` permutes the matri
 closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to
 make rows and columns more equal in norm. The default is `true` for both options.
 
+By default, the eigenvalues and vectors are sorted lexicographically by `(real(λ),imag(λ))`.
+A different comparison function `by(λ)` can be passed to `sortby`, or you can pass
+`sortby=nothing` to leave the eigenvalues in an arbitrary order.   Some special matrix types
+(e.g. `Diagonal` or `SymTridiagonal`) may implement their own sorting convention and not
+accept a `sortby` keyword.
+
 # Examples
 ```jldoctest
 julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
@@ -112,18 +133,18 @@ julia> vals == F.values && vecs == F.vectors
 true
 ```
 """
-function eigen(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T
+function eigen(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
     AA = copy_oftype(A, eigtype(T))
-    isdiag(AA) && return eigen(Diagonal(AA), permute = permute, scale = scale)
-    return eigen!(AA, permute = permute, scale = scale)
+    isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
+    return eigen!(AA; permute=permute, scale=scale, sortby=sortby)
 end
 eigen(x::Number) = Eigen([x], fill(one(x), 1, 1))
 
 """
-    eigvecs(A; permute::Bool=true, scale::Bool=true) -> Matrix
+    eigvecs(A; permute::Bool=true, scale::Bool=true, `sortby`) -> Matrix
 
 Return a matrix `M` whose columns are the eigenvectors of `A`. (The `k`th eigenvector can
-be obtained from the slice `M[:, k]`.) The `permute` and `scale` keywords are the same as
+be obtained from the slice `M[:, k]`.) The `permute`, `scale`, and `sortby` keywords are the same as
 for [`eigen`](@ref).
 
 # Examples
@@ -135,19 +156,17 @@ julia> eigvecs([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
  0.0  0.0  1.0
 ```
 """
-eigvecs(A::Union{Number, AbstractMatrix}; permute::Bool=true, scale::Bool=true) =
-    eigvecs(eigen(A, permute=permute, scale=scale))
+eigvecs(A::Union{Number, AbstractMatrix}; kws...) =
+    eigvecs(eigen(A; kws...))
 eigvecs(F::Union{Eigen, GeneralizedEigen}) = F.vectors
 
 eigvals(F::Union{Eigen, GeneralizedEigen}) = F.values
 
 """
-    eigvals!(A; permute::Bool=true, scale::Bool=true) -> values
+    eigvals!(A; permute::Bool=true, scale::Bool=true, sortby) -> values
 
 Same as [`eigvals`](@ref), but saves space by overwriting the input `A`, instead of creating a copy.
-The option `permute=true` permutes the matrix to become
-closer to upper triangular, and `scale=true` scales the matrix by its diagonal elements to
-make rows and columns more equal in norm.
+The `permute`, `scale`, and `sortby` keywords are the same as for [`eigen`](@ref).
 
 !!! note
     The input matrix `A` will not contain its eigenvalues after `eigvals!` is
@@ -171,29 +190,27 @@ julia> A
   0.0        5.37228
 ```
 """
-function eigvals!(A::StridedMatrix{<:BlasReal}; permute::Bool=true, scale::Bool=true)
-    issymmetric(A) && return eigvals!(Symmetric(A))
+function eigvals!(A::StridedMatrix{<:BlasReal}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
+    issymmetric(A) && return sorteig!(eigvals!(Symmetric(A)), sortby)
     _, valsre, valsim, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)
-    return iszero(valsim) ? valsre : complex.(valsre, valsim)
+    return sorteig!(iszero(valsim) ? valsre : complex.(valsre, valsim), sortby)
 end
-function eigvals!(A::StridedMatrix{<:BlasComplex}; permute::Bool=true, scale::Bool=true)
-    ishermitian(A) && return eigvals(Hermitian(A))
-    return LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)[2]
+function eigvals!(A::StridedMatrix{<:BlasComplex}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby)
+    ishermitian(A) && return sorteig!(eigvals(Hermitian(A)), sortby)
+    return sorteig!(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'N', 'N', A)[2], sortby)
 end
 
 # promotion type to use for eigenvalues of a Matrix{T}
 eigtype(T) = promote_type(Float32, typeof(zero(T)/sqrt(abs2(one(T)))))
 
 """
-    eigvals(A; permute::Bool=true, scale::Bool=true) -> values
+    eigvals(A; permute::Bool=true, scale::Bool=true, sortby) -> values
 
 Return the eigenvalues of `A`.
 
 For general non-symmetric matrices it is possible to specify how the matrix is balanced
-before the eigenvalue calculation. The option `permute=true` permutes the matrix to
-become closer to upper triangular, and `scale=true` scales the matrix by its diagonal
-elements to make rows and columns more equal in norm. The default is `true` for both
-options.
+before the eigenvalue calculation. The `permute`, `scale`, and `sortby` keywords are
+the same as for [`eigen!`](@ref).
 
 # Examples
 ```jldoctest
@@ -208,8 +225,8 @@ julia> eigvals(diag_matrix)
  4.0
 ```
 """
-eigvals(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T =
-    eigvals!(copy_oftype(A, eigtype(T)), permute = permute, scale = scale)
+eigvals(A::StridedMatrix{T}; kws...) where T =
+    eigvals!(copy_oftype(A, eigtype(T)); kws...)
 
 """
 For a scalar input, `eigvals` will return a scalar.
@@ -309,11 +326,11 @@ inv(A::Eigen) = A.vectors * inv(Diagonal(A.values)) / A.vectors
 det(A::Eigen) = prod(A.values)
 
 # Generalized eigenproblem
-function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasReal
-    issymmetric(A) && isposdef(B) && return eigen!(Symmetric(A), Symmetric(B))
+function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
+    issymmetric(A) && isposdef(B) && return eigen!(Symmetric(A), Symmetric(B), sortby=sortby)
     n = size(A, 1)
     alphar, alphai, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
-    iszero(alphai) && return GeneralizedEigen(alphar ./ beta, vr)
+    iszero(alphai) && return GeneralizedEigen(sorteig!(alphar ./ beta, vr, sortby)...)
 
     vecs = zeros(Complex{T}, n, n)
     j = 1
@@ -329,13 +346,13 @@ function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasReal
         end
         j += 1
     end
-    return GeneralizedEigen(complex.(alphar, alphai)./beta, vecs)
+    return GeneralizedEigen(sorteig!(complex.(alphar, alphai)./beta, vecs, sortby)...)
 end
 
-function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasComplex
-    ishermitian(A) && isposdef(B) && return eigen!(Hermitian(A), Hermitian(B))
+function eigen!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
+    ishermitian(A) && isposdef(B) && return eigen!(Hermitian(A), Hermitian(B), sortby=sortby)
     alpha, beta, _, vr = LAPACK.ggev!('N', 'V', A, B)
-    return GeneralizedEigen(alpha./beta, vr)
+    return GeneralizedEigen(sorteig!(alpha./beta, vr, sortby)...)
 end
 
 """
@@ -348,6 +365,9 @@ Computes the generalized eigenvalue decomposition of `A` and `B`, returning a
 
 Iterating the decomposition produces the components `F.values` and `F.vectors`.
 
+Any keyword arguments passed to `eigen` are passed through to the lower-level
+[`eigen!`](@ref) function.
+
 # Examples
 ```jldoctest
 julia> A = [1 0; 0 -1]
@@ -364,29 +384,29 @@ julia> F = eigen(A, B);
 
 julia> F.values
 2-element Array{Complex{Float64},1}:
- 0.0 + 1.0im
  0.0 - 1.0im
+ 0.0 + 1.0im
 
 julia> F.vectors
 2×2 Array{Complex{Float64},2}:
-  0.0-1.0im   0.0+1.0im
- -1.0-0.0im  -1.0+0.0im
+  0.0+1.0im   0.0-1.0im
+ -1.0+0.0im  -1.0-0.0im
 
 julia> vals, vecs = F; # destructuring via iteration
 
 julia> vals == F.values && vecs == F.vectors
 true
 ```
 """
-function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB}
+function eigen(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
     S = promote_type(eigtype(TA),TB)
-    return eigen!(copy_oftype(A, S), copy_oftype(B, S))
+    eigen!(copy_oftype(A, S), copy_oftype(B, S); kws...)
 end
 
 eigen(A::Number, B::Number) = eigen(fill(A,1,1), fill(B,1,1))
 
 """
-    eigvals!(A, B) -> values
+    eigvals!(A, B; sortby) -> values
 
 Same as [`eigvals`](@ref), but saves space by overwriting the input `A` (and `B`),
 instead of creating copies.
@@ -409,8 +429,8 @@ julia> B = [0. 1.; 1. 0.]
 
 julia> eigvals!(A, B)
 2-element Array{Complex{Float64},1}:
- 0.0 + 1.0im
  0.0 - 1.0im
+ 0.0 + 1.0im
 
 julia> A
 2×2 Array{Float64,2}:
@@ -423,15 +443,15 @@ julia> B
  0.0  1.0
 ```
 """
-function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasReal
-    issymmetric(A) && isposdef(B) && return eigvals!(Symmetric(A), Symmetric(B))
+function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
+    issymmetric(A) && isposdef(B) && return sorteig!(eigvals!(Symmetric(A), Symmetric(B)), sortby)
     alphar, alphai, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
-    return (iszero(alphai) ? alphar : complex.(alphar, alphai))./beta
+    return sorteig!((iszero(alphai) ? alphar : complex.(alphar, alphai))./beta, sortby)
 end
-function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}) where T<:BlasComplex
-    ishermitian(A) && isposdef(B) && return eigvals!(Hermitian(A), Hermitian(B))
+function eigvals!(A::StridedMatrix{T}, B::StridedMatrix{T}; sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
+    ishermitian(A) && isposdef(B) && return sorteig!(eigvals!(Hermitian(A), Hermitian(B)), sortby)
     alpha, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
-    alpha./beta
+    return sorteig!(alpha./beta, sortby)
 end
 
 """
@@ -453,13 +473,13 @@ julia> B = [0 1; 1 0]
 
 julia> eigvals(A,B)
 2-element Array{Complex{Float64},1}:
- 0.0 + 1.0im
  0.0 - 1.0im
+ 0.0 + 1.0im
 ```
 """
-function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) where {TA,TB}
+function eigvals(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}; kws...) where {TA,TB}
     S = promote_type(eigtype(TA),TB)
-    return eigvals!(copy_oftype(A, S), copy_oftype(B, S))
+    return eigvals!(copy_oftype(A, S), copy_oftype(B, S); kws...)
 end
 
 """
@@ -482,11 +502,11 @@ julia> B = [0 1; 1 0]
 
 julia> eigvecs(A, B)
 2×2 Array{Complex{Float64},2}:
-  0.0-1.0im   0.0+1.0im
- -1.0-0.0im  -1.0+0.0im
+  0.0+1.0im   0.0-1.0im
+ -1.0+0.0im  -1.0-0.0im
 ```
 """
-eigvecs(A::AbstractMatrix, B::AbstractMatrix) = eigvecs(eigen(A, B))
+eigvecs(A::AbstractMatrix, B::AbstractMatrix; kws...) = eigvecs(eigen(A, B; kws...))
 
 function show(io::IO, mime::MIME{Symbol("text/plain")}, F::Union{Eigen,GeneralizedEigen})
     summary(io, F); println(io)