Why are there type-unstable matrix operations in `Base`?

Internals & Design
1

In other domains, type stability is given more priority than generalizability, e.g. for sqrt with numbers. But not for matrices:

julia> code_warntype(sqrt, (Float64,)) # this is type-stable
MethodInstance for sqrt(::Float64)
  from sqrt(x::Union{Float32, Float64}) in Base.Math at math.jl:566
Arguments
  #self#::Core.Const(sqrt)
  x::Float64
Body::Float64
1 ─      nothing
│   %2 = Base.Math.zero(x)::Core.Const(0.0)
│   %3 = (x < %2)::Bool
└──      goto #3 if not %3
2 ─      Base.Math.throw_complex_domainerror(:sqrt, x)
└──      Core.Const(:(goto %7))
3 ┄ %7 = Base.Math.sqrt_llvm(x)::Float64
└──      return %7


julia> code_warntype(sqrt, (Matrix{Float64},)) # but this is not
MethodInstance for sqrt(::Matrix{Float64})
  from sqrt(A::StridedMatrix{T}) where T<:Union{Real, Complex} in LinearAlgebra at D:\.julia\juliaup\julia-1.7.0+0~x64\share\julia\stdlib\v1.7\LinearAlgebra\src\dense.jl:836
Static Parameters
  T = Float64
Arguments
  #self#::Core.Const(sqrt)
  A::Matrix{Float64}
Locals
  #33::LinearAlgebra.var"#33#34"
  SchurS::Schur{ComplexF64, Matrix{ComplexF64}}
  is_sqrt_real::Union{Missing, Bool}
  sqrtA::Union{Matrix{ComplexF64}, Matrix{Float64}}
  SchurF::Schur{Float64, Matrix{Float64}}
  sqrtHermA::Union{Symmetric{ComplexF64, Matrix{ComplexF64}}, Symmetric{Float64, Matrix{Float64}}}
Body::Union{Matrix{ComplexF64}, Matrix{Float64}}

What is the reason for that difference? And, can type-unstability be avoided in such operations with matrices (e.g. with alternative functions that throw an error if the arguments can’t yield results in the target domain)?

1 Like
2

See https://github.com/JuliaLang/julia/issues/4006 and linked discussions. I think the basic reasoning was that matrix square roots are so expensive that the cost of dynamic dispatch is not so big a concern in this case.

The workaround is to add a type assertion at the call site, I guess?

5 Likes
3

This is probably true (at least for large systems) but the instability may have an undesired effect on compile times as well I guess? In particular if it propagates far.

eigvals sometimes returning real values and sometimes complex have caused me a lot of headache when trying to improve compile times

2 Likes
4

It will always return real values if you give it an appropriate matrix type, e.g. eigvals(Hermitian(A)). If your matrix is not Hermitian or related, you might want to be cautious about using eigenvalues in the first place, since eigenvalues and eigenvectors can be very badly conditioned if the matrix is nearly defective.

In general, the main strategy for type stability is to indicate the type of the output via the type of the input. Unfortunately, with matrix square roots we don’t have a matrix type that guarantees that the output will be real.

2 Likes
5

In my case, most uses of eigvals are to calculate the poles of linear dynamical systems, they are often real without the matrix being Hermitian, and are central to a lot of concepts in control theory. The systems are often small as well, making the dynamic dispatch rather expensive in many cases.

Maybe using a separate definition of eigvals that always returns complex values in ControlSystems.jl would solve some of this headache.

3 Likes
6

You don’t necessarily need a type that guarantees real sqrt, in the same way that Float64 does not guarantee real sqrt: it could just throw an error if the return value is complex and advise the user to use sqrt(complex(A)) instead.

7

The problem with this is that there are faster algorithms for complex square roots of real matrices than there are for complex square roots of complex matrices.

Maybe we could introduce sqrt(A, Real) and sqrt(A, Complex) (and similarly for eigvals) to force a real or complex result, respectively (and throw an exception if this isn’t possible).

(Eventually we could contemplate defaulting to sqrt(A::AbstractMatrix{<:Real}) = sqrt(A, Real), but since that is a breaking change it might have to wait for Julia 2.)

4 Likes