broken eig functions for Diagonal{Matrix}

[stevengj]

Bug noted in JuliaLang/julia#21598, filing a separate issue so that it doesn't get lost.

julia> using LinearAlgebra

julia> I2 = Matrix(I, 2,2); D = Diagonal([2.0*I2, 3.0*I2])
2×2 Diagonal{Array{Float64,2},Array{Array{Float64,2},1}}:
 [2.0 0.0; 0.0 2.0]          ⋅         
         ⋅           [3.0 0.0; 0.0 3.0]

julia> eigvals(Diagonal([2.0*I2, 3.0*I2]))  # SHOULD RETURN ARRAY OF SCALARS!
2-element Array{Array{Float64,1},1}:
 [2.0, 2.0]
 [3.0, 3.0]

julia> eigvecs(D)
ERROR: MethodError: no method matching zero(::Type{Array{Float64,2}})
Closest candidates are:
  zero(::Type{LibGit2.GitHash}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LibGit2/src/oid.jl:220
  zero(::Type{Pkg.Resolve.VersionWeights.VersionWeight}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/Pkg/src/resolve/VersionWeights.jl:19
  zero(::Type{Pkg.Resolve.MaxSum.FieldValues.FieldValue}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/Pkg/src/resolve/FieldValues.jl:44
  ...
Stacktrace:
 [1] zeros(::Type{Array{Float64,2}}, ::Tuple{Int64,Int64}) at ./array.jl:467
 [2] Array{Array{Float64,2},2}(::UniformScaling{Bool}, ::Tuple{Int64,Int64}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/uniformscaling.jl:317
 [3] eigvecs(::Diagonal{Array{Float64,2},Array{Array{Float64,2},1}}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/diagonal.jl:486
 [4] top-level scope at none:0

julia> eigen(D)
ERROR: MethodError: no method matching isfinite(::Array{Float64,2})
Closest candidates are:
  isfinite(::BigFloat) at mpfr.jl:858
  isfinite(::Missing) at missing.jl:79
  isfinite(::Float16) at float.jl:545
  ...
Stacktrace:
 [1] (::getfield(Base, Symbol("##54#55")){typeof(isfinite)})(::Array{Float64,2}) at ./operators.jl:853
 [2] _any at ./reduce.jl:614 [inlined]
 [3] #any#548 at ./reducedim.jl:655 [inlined]
 [4] any at ./reducedim.jl:655 [inlined]
 [5] #eigen#114(::Bool, ::Bool, ::Function, ::Diagonal{Array{Float64,2},Array{Array{Float64,2},1}}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/diagonal.jl:488
 [6] eigen(::Diagonal{Array{Float64,2},Array{Array{Float64,2},1}}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/diagonal.jl:488
 [7] top-level scope at none:0

Should be quite easy to fix.

[AzamatB]

To clarify what is the general problem here, in your example D is a matrix whose elements are matrices, i.e. we are dealing with matrix over a field of diagonal nxn matrices.
Are you saying that the output of eigen(D) should return the same as

B = [2.0*I2 0I;
     0I 3.0*I2]
eigen(B)

?

[stevengj]

Yes, because we can also view it as a matrix over the scalar field (acting on blocked vectors of scalars).

[AzamatB]

Ok, I'm going to attempt this. For

julia> begin
           I2 = Matrix(I, 2,2)
           a = Matrix{Matrix{Float64}}(undef, 2,2)
           for i ∈ 1:2, j ∈ 1:2
               a[i,j] = (i+2j)I2
           end
           b = [a[1,1] a[1,2];
                a[2,1] a[2,2]]
       end
4×4 Array{Float64,2}:
 3.0  0.0  5.0  0.0
 0.0  3.0  0.0  5.0
 4.0  0.0  6.0  0.0
 0.0  4.0  0.0  6.0

julia> a
2×2 Array{Array{Float64,2},2}:
 [3.0 0.0; 0.0 3.0]  [5.0 0.0; 0.0 5.0]
 [4.0 0.0; 0.0 4.0]  [6.0 0.0; 0.0 6.0]

do we want eigen(a) to return the same result as eigen(b)?

[stevengj]

It seems like the eigenvectors should be "block" vectors (i.e. vectors of vectors).

[DShivansh]

@stevengj what i understand from your bug is that since the diagonal matrix in your example consists of two matrix so what you want is whenever we get matrix of matrices we should output the array consisting of the eigen vectors of those matrices.

So you want the output like
1.0 0.0
0.0 1.0

1.0 0.0
0.0 1.0

that is the two vectors that is the eigen vector of each matrix inside that of diagonal matrix.

@stevengj can you please conform whether I am interpreting the bug correctly or not.

thanks

[karthikrangasai]

@stevengj , I have submitted a PR, could you please review it