Julia: Type instability in eig()

[Linusnie]

There seems to be an issue with getindex for the type Eigen resulting in eig(::Matrix) being type unstable:

julia> @code_warntype eig(rand(5,5))
Variables:
  A::Array{Float64,2}
  args::Tuple{}

Body:
  begin $(Expr(:line, 66, symbol("linalg/eigen.jl"), symbol("")))
      GenSym(8) = (top(ccall))(:jl_alloc_array_1d,$(Expr(:call1, :(top(apply_type)), :Array, Any, 1)),$(Expr(:call1, :(top(svec)), :Any, :Int)),Array{Any,1},0,0,0)::Array{Any,1}
      GenSym(9) = GenSym(8)
      return __eig#179__(GenSym(9),A::Array{Float64,2})::Tuple{Any,Any}
  end::Tuple{Any,Any}

The output is assigned Tuple(Any,Any) instead of Tuple(Number,Number) as you might expect.

[jiahao]

In short, this type instability is harmless. Performance-sensitive code should use eigfact and the Eigen objects directly. Also, the return type of eig is _never_ Tuple{Number, Number} or any subtype thereof.

Longer explanation:

There seems to be an issue with getindex for the type Eigen

I think you misunderstand the origin of the type instability. You can verify for yourself with specific instances of the Eigen object that the return type is as tightly typed as possible:

julia> A=eigfact(randn(5,5)); #Nonsymmetric eigenproblem - returns complex data

julia> @code_warntype(A[:vectors])
Variables:
  A::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}}
  d::Symbol

Body:
  begin  # linalg/eigen.jl, line 21:
      unless d::Symbol === :values::Bool goto 0
      return (top(getfield))(A::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}},:values)::Array{Complex{Float64},1}
      0:  # line 22:
      unless d::Symbol === :vectors::Bool goto 1
      return (top(getfield))(A::Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}},:vectors)::Array{Complex{Float64},2}
      1:  # line 23:
      return (Base.LinAlg.throw)($(Expr(:new, :((top(getfield))(Base,:KeyError)::Type{KeyError}), :(d::Symbol))))::Union{}
  end::Union{Array{Complex{Float64},1},Array{Complex{Float64},2}}

The return type is either a 1-D or 2-D array of complexes, depending on the value of the index being called. Because :value and :vector are not lifted into the type system, the type inference algorithm only knows that getfield(::Eigen, ::Symbol) is being called.

The real cause is that eigfact (called by eig) returns either an Eigen object with complex numeric entries, or real numeric entries, depending on the runtime value of the input matrix:

julia> @code_warntype eigfact(rand(5,5))
Variables:
  A::Array{Float64,2}
  ##S#6437::Type{Float64}

Body:
  begin $(Expr(:line, 56, symbol("linalg/eigen.jl"), symbol("")))
      GenSym(0) = (Base.LinAlg.norm)((Base.box)(Float64,(Base.sitofp)(Float64,1)::Any)::Float64)::Float64
      ##S#6437 = Float64
      return (top(kwcall))((top(getfield))(Base.LinAlg,:call)::F,2,:permute,true,:scale,true,Base.LinAlg.eigfact!,(top(ccall))(:jl_alloc_array_1d,(top(apply_type))(Base.Array,Any,1)::Type{Array{Any,1}},(top(svec))(Base.Any,Base.Int)::SimpleVector,Array{Any,1},0,4,0)::Array{Any,1},(Base.LinAlg.copy)(A::Array{Float64,2})::Array{Float64,2})::Union{Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}},Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}}
  end::Union{Base.LinAlg.Eigen{Complex{Float64},Complex{Float64},Array{Complex{Float64},2},Array{Complex{Float64},1}},Base.LinAlg.Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}}

From this output you can see that when eig returns the eigenpairs, you'd want the return type to be inferred as something like

(F[:values], F[:vectors]) :: Union{Tuple{Vector{Complex128}, Matrix{Complex128}}, Tuple{Vector{Float64}, Matrix{Float64}}}

(not Tuple{Number, Number} as you wrote.)

However the type inference algorithm currently knows to use the runtime value of the Symbol in getindex(::Eigen, ::Symbol), so instead it infers at some point something like

(F[:values], F[:vectors]) :: Tuple{Z, Z} #getindex(::Eigen, ::Symbol) is called twice

Z = Union{
   Tuple{Vector{Complex128}, Matrix{Complex128}}, #complex case
   Tuple{Vector{Float64}, Matrix{Float64}} #real case
}

and Z is sufficiently complicated that the type inference algorithm gives up recording the specifics of Z and widens it to Any.

Summary: Tuple{Any,Any} is returned by widening an excessively complicated inferred type.

[mauro3]

I do find the design of the getindex method https://github.com/JuliaLang/julia/blob/8a7752d7fb0aaad4670f13a6763785cbc3758a27/base/linalg/eigen.jl#L20 a bit odd: The values and vectors of Eigen don't necessarily have the same type, thus getindex is type unstable. Further the getindex method seems only used internally, so why not just access the fields directly?

[jiahao]

Also it is useful to cross-reference mailing list discussions. Ref: https://groups.google.com/forum/#!topic/julia-users/rXSfEQV8WTM

[mschauer]

@mauro3 There is a lack of a syntactically blessed Val{:L}, where do we have an issue for that?

[mauro3]

I guess that goes into #1974, could be an alternative to dot-overloading.

[IainNZ]

This came up on StackoverFlow just now too: http://stackoverflow.com/questions/31633153/peformance-difference-due-to-type-instability

The Tuple{Any,Any} typing propagated to the length of eigenvalue vector being Any, when then propagated further. A bit unfortunate.

[andyferris]

I'm still confused why eigfact shouldn't be type-stable. In the case of sqrt it is quite clear: you want a type stable result.

I feel we should solve this. Some options:

1. Have eigfact() always return complex output, and let people use Symmetric and Hermitian to get more efficiency
2. Have eigfact() throw an error like sqrt() does, for inappropriate (real, non-symmetric) input.
3. Since the overhead of casting a real input to a complex input isn't strictly necessary, define a method like eig_complex for when the output is complex.
4. Optionally let the user define the output eltype. E.g. eigfact(Float64, mat::Matrix{Float64}) would work if and only if the matrix is symmetric. Users could type eigfact(Complex{Float64}, mat::Matrix{Float64}) in the general case. Similarly for eig().

I like 1 and 4 the best. I don't feel these are particularly onerous for the user, and 4 is non-breaking.

[andreasnoack]

I'd be okay with 1. It would be much more robust to ask people to annotate their symmetric matrix with Symmetric. The slightest non-symmetry in a matrix will result in a call to the non-symmetric solver which is slow but the main problem is that you loose the sorting of the eigenvalues. The main problem with 1. is that some non-symmetric matrices have real eigenvalues but they are relatively rare. At least for reasonable sizes matrices.

[andyferris]

I'd be okay with 1.

I've been leaning more and more this way.

The main problem with 1. is that some non-symmetric matrices have real eigenvalues but they are relatively rare.

I had thought of that. Are there specialized LAPACK routines for these cases? If so, maybe a simple wrapper for RealDiagonlizableMatrix would be sufficient to plug the gap in the relatively small cases where it is important.

[nalimilan]

+1 for solution 1 too.

[andreasnoack]

I had thought of that. Are there specialized LAPACK routines for these cases? If so, maybe a simple wrapper for RealDiagonlizableMatrix would be sufficient to plug the gap in the relatively small cases where it is important.

LAPACK returns the result of the non-symmetric real eigenvalue problem in a packed format. We could add a special type for handling this. The point being that when the solution is complex, it still has a lot of structure, i.e. the complex values are conjugate pairs. However, I don't think it is really worth it. If you are serious about your numerical algorithm you'd probably avoid the non-symmetric eigendecomposition anyway and just stick with the real Schur.

[andyferris]

That makes sense, Andreas.

Anyway, we could move ahead with option 1 and if it turns out lots of people want to take advantage of the LAPACK non-symmetric eigenvalue algorithm and its packed format, we can implement that special wrapper type later.

[dlfivefifty]

Have eigfact() throw an error like sqrt() does, for inappropriate (real, non-symmetric) input.

It's not only symmetric matrices that have real eigenvalues, and I don't think it's possible to tell a priori:

julia> A=rand(10,10); A = A + A'; D=diagm(rand(10)); A = D*A*inv(D); issymmetric(A)
false

julia> eigvals(A)
10-element Array{Float64,1}:
 10.0053  
 -1.90995 
  2.24982 
 -0.959746
 -0.749866
 -0.183301
  0.212858
  0.763341
  1.06319 
  1.54392 

[nalimilan]

@dlfivefifty The idea is that you would get a complex result in that case, unless you wrap the input matrix in a Symmetric object.

[dlfivefifty]

It's not a symmetric matrix, so why would you wrap it in Symmetric?

[andreasnoack]

@dlfivefifty The point is that the price of type stability could be that your result becomes complex with zero imaginary part or an error. If you have a symmetric matrix, you probably want real eigenvalues. To get that in a type stable way, you'd have to use the Symmetric wrapper.

[dlfivefifty]

I know you're point. My point is sometimes you know your answer is real even when it's not symmetric. Above someone proposed RealDiagonalizable For this situation. Another solution is to return a special array ComplexConjugateVector <: Vector{Complex128} that uses the optimal LAPACK storage

[andreasnoack]

My judgment based on the matrices I've worked with over the years is that in practice and for non-toy problems, if the matrix is not symmetric, then you might as well just assume that it as complex eigenvalues. The ComplexConjugateVector solution is appealing but I'm not sure it is worth the work required.

[dlfivefifty]

Will have to think if there's a good example... (see 2x deleted comments of bad examples)

Whether it's worth the work depends on whether the extra memory usage is going to impact algorithms. I could see this happening.