# 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