`sqrt(::Matrix{Float64})` is not type stable

[LilithHafner]

julia> using StableRNGs

julia> rng = StableRNG(0)
StableRNGs.LehmerRNG(state=0x00000000000000000000000000000001)

julia> sqrt(rand(rng,3,3))
3×3 Matrix{ComplexF64}:
    0.410548+0.369007im      0.49482-0.209855im   0.394519-0.231836im
 -0.00928873+0.0595217im    0.935186-0.0338502im   0.11389-0.0373957im
    0.661866-0.458654im   -0.0589303+0.260838im   0.524852+0.288158im

julia> sqrt(rand(rng,3,3))
3×3 Matrix{Float64}:
 0.301226   0.392135   -0.452243
 0.183647   0.25732     1.14379
 0.185353  -0.0785732   0.578105

For scalars, we throw when sqrting a negative unless the negative is of a complex type in which case we return a complex result. IMO that is a much better solution. Switching to that would be breaking in the case of square rooting real matrices with negative real eigenvalues.

[jishnub]

Previous discussion in #21

[giordano]

Shouldn't this be in the LinearAlgebra repo?

[stevengj]

A non-breaking change could be to pass an additional argument, e.g. sqrt(A, Real) if you want a real result (throwing an exception if not real) and sqrt(A, Complex) if you want a complex result (even if the result is real).

Note that eigen (and eigvals and eigvecs) have the same issue — for a real matrix, they return real matrices/vectors if possible and complex matrices/vectors if not. So any change to sqrt should probably be mirrored there.

[araujoms]

And also some functions defined in symmetric.jl. Namely ^, log, acos, asin, and acosh.

[aravindh-krishnamoorthy]

For scalars, we throw when sqrting a negative unless the negative is of a complex type in which case we return a complex result. IMO that is a much better solution. Switching to that would be breaking in the case of square rooting real matrices with negative real eigenvalues.

I might be in the minority here, but I'm not a fan of sqrt(-1) returning an error. IMHO it seems like a case of prioritising language features over functionality.

In fact, in this case, returning an error would be worse because it's not obvious (or easily computable) whether the eigenvalues of a real-valued input matrix are such that the square root will be real-valued. So, we'll always need to code-in two calls.

[stevengj]

but I'm not a fan of sqrt(-1) returning an error. IMHO it seems like a case of prioritising language features over functionality.

It was a necessary sacrifice: having a basic scalar function like that be type-unstable would completely kill performance. Whereas with matrices the cost of dynamic dispatch is not as significant.

In fact, in this case, returning an error would be worse because it's not obvious (or easily computable) whether the eigenvalues of a real-valued input matrix are such that the square root will be real-valued.

I suspect that in most real-world cases where you are taking the matrix square root, and in which you care about performance, you will know whether the square root is necessarily real — typically because the matrix was computed via $A^T A$ or similar and hence is guaranteed to be positive (semi)definite. (Though in many such cases you could use Cholesky instead.)

Whereas if the matrix square root might be complex and you are okay with that, you might be fine with just returning a complex matrix all the time, i.e. you might call sqrt(A, Complex) or similar.

In any case, because of backwards compatibility, we can't change the default sqrt(A) — it must remain type-unstable. But we could add an extra argument to indicate the desired return type.

It would be helpful to have a use-case here, though. Are there any good examples of performance-sensitive applications using a matrix square root?

[LilithHafner]

We should re-open this iff someone provides such an example.