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

**URL:** <https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999>\
**Category:** Internals & Design\
**Tags:** linearalgebra\
**Created:** [April 4, 2022, 10:41am UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999 "2022-04-04T10:41:50Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![heliosdrm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/heliosdrm/32/3851_2.png) [@heliosdrm](https://discourse.julialang.org/u/heliosdrm)\
**Post date:** [April 4, 2022, 10:41am UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999/1 "2022-04-04T10:41:50Z")

</div>

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

```julia-repl
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)?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 4, 2022, 11:26am UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999/2 "2022-04-04T11:26:40Z")

</div>

See [https://github.com/JuliaLang/julia/issues/4006](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?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 4, 2022, 2:04pm UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999/3 "2022-04-04T14:04:47Z")

</div>

> [@stevengj](#):
>
> 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.

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 4, 2022, 2:49pm UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999/4 "2022-04-04T14:49:16Z")

</div>

> [@baggepinnen](#):
>
> `eigvals` sometimes returning real values and sometimes complex have caused me a lot of headache when trying to improve compile times

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.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 4, 2022, 2:54pm UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999/5 "2022-04-04T14:54:56Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [April 5, 2022, 7:15am UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999/6 "2022-04-05T07:15:12Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [April 5, 2022, 12:06pm UTC](https://discourse.julialang.org/t/why-are-there-type-unstable-matrix-operations-in-base/78999/7 "2022-04-05T12:06:56Z")

</div>

> [@dlfivefifty](#):
>
> advise the user to use `sqrt(complex(A))` instead.

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.)