Breaking change in sqrt for Matrix which is secretly Diagonal

[lkdvos]

Copied here of a discussion I opened on Discourse:

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

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

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