Generalized eigenvalues with Diagonal B

Looking to solve a generalized eigenproblem. All is well here:

julia> A = rand(4,4);
julia> B = rand(4,4);
julia> eigen(A,B)
4-element Array{Complex{Float64},1}:
  -1.7837186962342586 + 0.0im               
 -0.24734999477221975 - 0.3007940273651078im

But then I ran into this:

julia> B = Diagonal(1:4);
julia> eigen(A,B)
ERROR: MethodError: no method matching eigen!(::Array{Float64,2}, ::Diagonal{Float64,Array{Float64,1}})
Closest candidates are:

An explicit cast to a matrix makes it all better:

julia> eigen(A,Matrix(B))
4-element Array{Float64,1}:

But I expected the original to work, given that

julia> B isa AbstractMatrix

and eigen(B) is fine, too. So is this the intended behavior?

Our generalized eigensolvers are currently just calling LAPACK so the supported signatures correspond to what LAPACK offers. It wouldn’t be hard to add a specialized method for diagonal right hand side but I think you are the first one to ask for it. It would be great if you could make a PR that adds support for diagonal right hand side.


You can also call C = sqrt.(B); eigen(inv(C) * A * inv(C)), and then multiply the eigenvectors by inv(C) to get the generalized eigenvectors. The eigenvalues are the same.


I’m willing to look into that, but what surprised me is that there is no dispatch on an AbstractMatrix for B that would supply a default fallback behavior. I’m sure there are aspects to doing that I don’t yet appreciate.

Basically it is because generic dense eigensolver algorithms are very complicated to implement well. There is one in the GenericLinearAlgebra package, but it will be much slower than the LAPACK version for Matrix{Float64} and (AFAIK) doesn’t take advantage of special sparsity structures like diagonal B.