Type instability and/or behavior instability of `sqrt` on nightly

Internals & Design
,
1

I noticed that since this commit for LinearAlgebra.jl, the default implementation of sqrt will now check if a matrix is diagonal to improve performance.
There is however a change in behavior associated to this which might be surprising to users, related to the difference in implementation for non-positive definite diagonal arrays or regular matrices.

In particular, the behavior in Julia 1.11 and before is that sqrt(::Matrix{Real}) returns a real matrix whenever the input is positive definite, and a complex matrix whenever it is not. On the other hand, sqrt(::Diagonal{Real}) always returns a real diagonal, and thus throws an error when some elements are negative.

As a result, after the change in the mentioned commit, now sqrt(::Matrix{Real}) will return a complex matrix when it is non-positive definite but not diagonal, throw an error for diagonal non-positive definite, and return a real matrix for positive definite matrices. This seems to be slightly undesirable.

I’m aware there have been previous discussions and concerns about the type-instability of sqrt, and I do agree that the current design choices are not great, but I would argue to keep consistency between Diagonal and Matrix, especially if the Matrix implementation will dispatch through to a Diagonal one whenever possible.

(v1.11)

Summary 
julia> d = Diagonal([1., -1.])
2×2 Diagonal{Float64, Vector{Float64}}:
 1.0    ⋅
  ⋅   -1.0

julia> sqrt(d)
ERROR: DomainError with -1.0:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

julia> sqrt(Array(d))
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+1.0im

julia> sqrt(Array(d) + [0. 0.1; 0. 0.]) # array but not strictly diagonal
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.05-0.05im
 0.0+0.0im   0.0+1.0im

(v1.12)

Summary 
julia> using LinearAlgebra

julia> d = Diagonal([1., -1.])
2×2 Diagonal{Float64, Vector{Float64}}:
 1.0    ⋅
  ⋅   -1.0

julia> sqrt(d)
ERROR: DomainError with -1.0:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

julia> sqrt(Array(d))
ERROR: DomainError with -1.0:
sqrt was called with a negative real argument but will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).

julia> sqrt(Array(d) + [0. 0.1; 0. 0.]) # array but not strictly diagonal
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.05-0.05im
 0.0+0.0im   0.0+1.0im

(Not sure if this had to be an issue or a discussion post here…)

2 Likes
2

I was looking at sqrt the other day (for triangular matrices, in my case) and also noticed the instability. My initial opinion (which is not the result of extensive thought) was that it should do the thing we do for sqrt(::Real). I.e., it should throw for a real matrix with nonreal sqrt and demand that the input matrix be converted to complex instead.

1 Like
3

I noticed that time-to-load became super slow in a package or two of mine. It was due to this. I assumed it was a known issue.

4

Just as a small aside, it’s not necessarily the type instability that is the biggest problem here to me, while I have opinions about that as well, at least I know what’s going on. The thing that is a bit problematic here is that I think this is really a breaking change, since I can now no longer assume that sqrt(::Array) will always return (even though it probably shouldn’t). I’ll open an issue in LinearAlgebra.jl to alert people about the breaking change.

1 Like