# Square root of a matrix that has no square root

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.

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.

WolframAlpha produces the same result happily and quietly:

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

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