LinearAlgebra eigen errors when given Diagonal

[KronosTheLate]

When I pass matrices of type Diagonal to eigen, I get errors. This should work, right?

Demonstrating the issue

julia> using LinearAlgebra

julia> A = Diagonal([1, 2])
2×2 Diagonal{Int64, Vector{Int64}}:
 1  ⋅
 ⋅  2

julia> eigen(A, A)
ERROR: MethodError: no method matching eigen!(::Diagonal{Float64, Vector{Float64}}, ::Diagonal{Float64, Vector{Float64}})
Closest candidates are:
  eigen!(::Union{Hermitian{T, S}, Hermitian{Complex{T}, S}, Symmetric{T, S}}, ::AbstractMatrix{T}; sortby) where {T<:Number, S<:(StridedMatrix{T} where T)} at C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\symmetriceigen.jl:170
Stacktrace:
 [1] eigen(A::Diagonal{Int64, Vector{Int64}}, B::Diagonal{Int64, Vector{Int64}}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ LinearAlgebra C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:511
 [2] eigen(A::Diagonal{Int64, Vector{Int64}}, B::Diagonal{Int64, Vector{Int64}})
   @ LinearAlgebra C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:509
 [3] top-level scope
   @ REPL[3]:1

julia> eigen(A|>Matrix, A)
ERROR: MethodError: no method matching eigen!(::Matrix{Float64}, ::Diagonal{Float64, Vector{Float64}})
Closest candidates are:
  eigen!(::StridedMatrix{T}; permute, scale, sortby) where T<:Union{Float32, Float64} at C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:149
  eigen!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64} at C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:436
  eigen!(::Union{Hermitian{T, S}, Hermitian{Complex{T}, S}, Symmetric{T, S}}, ::AbstractMatrix{T}; sortby) where {T<:Number, S<:(StridedMatrix{T} where T)} at C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\symmetriceigen.jl:170
Stacktrace:
 [1] eigen(A::Matrix{Int64}, B::Diagonal{Int64, Vector{Int64}}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ LinearAlgebra C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:511
 [2] eigen(A::Matrix{Int64}, B::Diagonal{Int64, Vector{Int64}})
   @ LinearAlgebra C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:509
 [3] top-level scope
   @ REPL[4]:1

julia> eigen(A, A|>Matrix)
ERROR: MethodError: no method matching eigen!(::Diagonal{Float64, Vector{Float64}}, ::Matrix{Float64})
Closest candidates are:
  eigen!(::StridedMatrix{T}, ::StridedMatrix{T}; sortby) where T<:Union{Float32, Float64} at C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:436
  eigen!(::Union{Hermitian{T, S}, Hermitian{Complex{T}, S}, Symmetric{T, S}}, ::AbstractMatrix{T}; sortby) where {T<:Number, S<:(StridedMatrix{T} where T)} at C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\symmetriceigen.jl:170
Stacktrace:
 [1] eigen(A::Diagonal{Int64, Vector{Int64}}, B::Matrix{Int64}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ LinearAlgebra C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:511
 [2] eigen(A::Diagonal{Int64, Vector{Int64}}, B::Matrix{Int64})
   @ LinearAlgebra C:\Users\Dennis Bal\.julia\juliaup\julia-1.8.0+0.x64\share\julia\stdlib\v1.8\LinearAlgebra\src\eigen.jl:509
 [3] top-level scope
   @ REPL[5]:1

julia> eigen(A|>Matrix, A|>Matrix)
GeneralizedEigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 0.9999999999999998
 1.0
vectors:
2×2 Matrix{Float64}:
 0.0       1.0
 0.707107  0.0

julia> versioninfo()
Julia Version 1.8.0
Commit 5544a0fab7 (2022-08-17 13:38 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × Intel(R) Core(TM) i5-8265U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
  Threads: 8 on 8 virtual cores

[pcjentsch]

I think this should work yes.

The easy fix would be to define

eigen!(A::AbstractMatrix,B::Diagonal; sortby::Union{Function,Nothing}=nothing) = eigen!(A, Matrix(B);sortby)
eigen(A::Diagonal,B::AbstractMatrix; sortby::Union{Function,Nothing}=nothing) = eigen(Matrix(A), B;sortby)
eigen(A::Diagonal,B::Diagonal;sortby::Union{Function,Nothing}=nothing) = eigen(Matrix(A), Matrix(B);sortby)

in LinearAlgebra. eigen!(A::Diagonal, B::AbstractMatrix) isn't defined because we can't use Diagonal as a output array in LAPACK.

Ideally, we would have more efficient methods for diagonal matrices. In theory, we should just be able to do something like eigen(A::AbstractMatrix,B::Diagonal) = eigen(inv(B)*A), which would be efficient because inv and * are fast for diagonal matrices.

However, I feel like eigen(A::Diagonal,B::AbstractMatrix) == eigen(Matrix(A),B) should hold, and I think the old issue #230 is making that harder than it should be.

[dkarrasch]

Creating dense arrays seems a bit wasteful. I think one can do much better in the case when the Diagonal is the second argument/rhs, as you say @pcjentsch. I think I have a draft PR, just need some tests.

[stevengj]

Duplicate of #594?

[KronosTheLate]

Sure looks like a duplicate to me.

[dkarrasch]

No, it's not a duplicate. JuliaLang/julia#46797 fixes this issue here, but not the other one. The old issue is about standard eigendecompositions of matrices of matrices, whereas this one is about generalized eigendecomposition of abstract matrices.