Square root of a matrix that has no square root

General Usage
1

Is this expected?

julia> sqrt([0 1; 0 0])
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

I thought it would give an error.

2

Throwing an error here would be hard. The sqrt is defined as long as the eigenvalues are unique, and numerically determining this would be really unreliable.

3

WolframAlpha produces the same result happily and quietly:

MatrixPower[{{0, 1}, {0, 0}}, 1/2]
#result:
{{0, 0}, {0, 0}}
4

It’s fine to have repeated eigenvalues, even for a defective matrix:

julia> A = [4 1; 0 4]
2×2 Matrix{Int64}:
 4  1
 0  4

julia> S = sqrt(A)
2×2 Matrix{Float64}:
 2.0  0.25
 0.0  2.0

julia> S^2
2×2 Matrix{Float64}:
 4.0  1.0
 0.0  4.0

The problem only arises for repeated eigenvalues that are exactly zero and defective, I think?

It seems like we could check for this case in _sqrt_quasitriu_diag_block! — perhaps you could just set a flag on this line when iszero(A[i,i]), and throw an error if you hit that condition twice in a row (assuming the diagonal entries are sorted, as usual for Schur).

Of course, roundoff errors may defeat this too if the matrix is not in Schur-like form to start with:

julia> A = [0 1; 0 0]
2×2 Matrix{Int64}:
 0  1
 0  0

julia> sqrt(A) # fail
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

julia> B = rand(2,2); # random change of basis

julia> A2 = B * A / B # similar to A, should also not have a square root
2×2 Matrix{Float64}:
 0.68983  -0.154736
 3.07534  -0.68983

julia> S = B \ sqrt(A2) * B # change basis from √A2 back to original basis
2×2 Matrix{ComplexF64}:
   7.63631e-5+1.68025e-20im     6547.66+5.55286e-13im
 -1.34717e-12-3.78013e-21im  7.63631e-5-4.38691e-13im

julia> S^2 # should == A if square root existed
2×2 Matrix{ComplexF64}:
   -2.9895e-9-2.4751e-17im          1.0-2.8724e-9im
 -2.05748e-16+1.36679e-26im  -2.9895e-9-9.17507e-17im

julia> norm(S^2 - A) / norm(A) # loses half of the significant digits due to ill-conditioning
3.147765149816803e-8
3 Likes