Eigen decomposition of Diagonal matrix returns Dense eigenvectors

[jecs]

Hello. I noticed today that if you have a matrix of Diagonal type, then the eigenvector matrix (from eigvecs or eigen) is not of type Diagonal. Wouldn't it make more sense for the eigenvectors to match the type of the input matrix, when possible?

For example, if I have A = Diagonal(rand(10)), then eigvecs(A;kwarg...) should return a Diagonal (e.g. one(A)), but eigvecs(Matrix(A);kwarg...) should return a Matrix (e.g. one(Matrix(A))), and so on. The same changes would be needed for eigen(...).

[fredrikekre]

This would be incompatible with the sortby kwarg since the return type would have to depend on the value of sortby then.

[jecs]

Maybe for Diagonal, using sortby could return view(one(A),:,p) where p is the permutation? Or would this introduce type instability?

[eschnett]

Isn't eigenvector calculation type unstable anyway? eigvecs checks whether the matrix is symmetric (by checking the values in the matrix), and if so, returns a real-valued result. Thus a dependence on a keyword argument shouldn't matter.

julia> A = randn(4,4);

julia> typeof(eigvecs(A))
Matrix{ComplexF64} (alias for Array{Complex{Float64}, 2})

julia> typeof(eigvecs(A+A'))
Matrix{Float64} (alias for Array{Float64, 2})