# Matrix square root of known positive semi-definite matrix does not exist?

In the following example, I believe that the matrix `C` should be positive semi-definite by construction. (Note that the matrix `B` is positive semi-definite.) As a consequence, it should have a unique positive semi-definite square root matrix.
However, Julia returns a matrix of complex `NaN`'s.
Can anyone help me out with what is going on here?

``````using LinearAlgebra

A = [
0.9999999999999998	4.649058915617843e-16	-1.3149405273715513e-16	9.9959579317056e-17
-8.326672684688674e-16	1.0000000000000004	2.9280733590254494e-16	-2.9993900031619594e-16
9.43689570931383e-16	-1.339206523454095e-15	1.0000000000000007	-8.550505126287743e-16
-6.245004513516506e-16	-2.0122792321330962e-16	1.183061278035052e-16	1.0000000000000002
]

B = [
0.09648289218436859	0.023497875751503007	0.0	0.0
0.023497875751503007	0.045787575150300804	0.0	0.0
0.0	0.0	0.0	0.0
0.0	0.0	0.0	0.0
]

display(sqrt(B)) # no problem

C = A * B * A' # should still be positive semi-definite
display(sqrt(C)) # does not exist???
``````

The output:

``````4×4 Array{Float64,2}:
0.307265   0.0455076  0.0  0.0
0.0455076  0.209085   0.0  0.0
0.0        0.0        0.0  0.0
0.0        0.0        0.0  0.0
4×4 Array{Complex{Float64},2}:
NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im
NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im
NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im
NaN+NaN*im  NaN+NaN*im  NaN+NaN*im  NaN+NaN*im
``````

The problem is that `A` and `B` have a lot of detail that is not captured by `Float64`'s. As such, when you take `A*B*A'`The result is not positive semi-definite.

``````C=A*B*A'
4×4 Array{Float64,2}:
0.0964829     0.0234979     5.95814e-17  -6.4982e-17
0.0234979     0.0457876    -3.91443e-17  -2.38882e-17
5.95814e-17  -3.91443e-17   1.08649e-31  -2.93317e-32
-6.4982e-17   -2.38882e-17  -2.93317e-32   4.53883e-32
``````

As such, `sqrt(C)` does not exist. The basic problem here is that all of your matrices suck from a stability perspective.

1 Like

Ok, stability suckiness is duly noted.
Any general suggestions on how to improve it?

Maybe some high-level background on where these matrices come from…
B is an estimated variance-covariance matrix.
`A` is obtained as `D*pinv(D)` where `D` is an input matrix that contains nice, well-rounded entries.
(I would include it, but it is 4 x 102)

My initial thought was that it is the `pinv` that is creating stability problems.
But this only happens on very rare draws, and `D` is staying fixed over draws.
Only `B` is changing.

If you know by construction that it is symmetric and positive definite, you can try `real(Symmetric(C))`.

With these properties you also have the option to do the naive square root via eigenvalue factorization.

1 Like

That seems to solve it.
I get a complex result, albeit with tiny imaginary parts.
Presumably there’s nothing wrong with converting it to real afterwards?

``````julia> using LinearAlgebra

julia> A = [
0.9999999999999998        4.649058915617843e-16        -1.3149405273715513e-16        9.9959579317056e-17
-8.326672684688674e-16        1.0000000000000004        2.9280733590254494e-16        -2.9993900031619594e-16
9.43689570931383e-16        -1.339206523454095e-15        1.0000000000000007        -8.550505126287743e-16
-6.245004513516506e-16        -2.0122792321330962e-16        1.183061278035052e-16        1.0000000000000002
];

julia> B = [
0.09648289218436859        0.023497875751503007        0.0        0.0
0.023497875751503007        0.045787575150300804        0.0        0.0
0.0        0.0        0.0        0.0
0.0        0.0        0.0        0.0
];

julia> display(sqrt(B)) # no problem
4×4 Array{Float64,2}:
0.307265   0.0455076  0.0  0.0
0.0455076  0.209085   0.0  0.0
0.0        0.0        0.0  0.0
0.0        0.0        0.0  0.0

julia> C = A * B * A'; # should still be positive semi-definite

julia> C = real(Symmetric(C)); # enforce positive semi-def

julia> Crt = sqrt(C);

julia> display(Crt) # now it exists
4×4 Symmetric{Complex{Float64},Array{Complex{Float64},2}}:
0.307265+1.85967e-39im     0.0455076-2.89282e-39im   2.88901e-16-1.86116e-24im  -1.98203e-16+5.50328e-39im
0.0455076-2.89282e-39im      0.209085+4.49994e-39im  -2.42324e-16+2.89513e-24im  -7.77374e-17-8.56065e-39im
2.88901e-16-1.86116e-24im  -2.42324e-16+2.89513e-24im   6.73366e-31+1.86265e-9im   -1.18185e-31-5.50768e-24im
-1.98203e-16+5.50328e-39im  -7.77374e-17-8.56065e-39im  -1.18185e-31-5.50768e-24im   1.39421e-31+1.62857e-38im

julia> display(real(Crt))
4×4 Array{Float64,2}:
0.307265      0.0455076     2.88901e-16  -1.98203e-16
0.0455076     0.209085     -2.42324e-16  -7.77374e-17
2.88901e-16  -2.42324e-16   6.73366e-31  -1.18185e-31
-1.98203e-16  -7.77374e-17  -1.18185e-31   1.39421e-31

``````

As noted by others above, roundoff errors give you tiny negative eigenvalues. One thing you can do is a slight shift:

``````julia> sqrt(Hermitian(C) + 1e-15*norm(C)*I)
4×4 Symmetric{Float64,Array{Float64,2}}:
0.307265      0.0455076     2.96653e-16  -1.98203e-16
0.0455076     0.209085     -2.62089e-16  -7.77374e-17
2.96653e-16  -2.62089e-16   1.04541e-8    3.22541e-25
-1.98203e-16  -7.77374e-17   3.22541e-25   1.05758e-8
``````

It might be nice to do this automatically for Hermitian matrices when a tiny negative eigenvalue is detected. Here is a draft PR if people want to bang on it: RFC: treat small negative λ as 0 for sqrt(::Hermitian) by stevengj · Pull Request #35057 · JuliaLang/julia · GitHub

1 Like

It always exists but may be complex.

At the very worst, this should give a complex result — returning `NaN` values here seems like a bug: https://github.com/JuliaLang/julia/issues/35058

1 Like

I can see the argument for the current behavior which is `sqrt(-1)` is also Nan despite being defined in complex arithmetic. Otoh, I can see how linear algebra tends to be more complex number heavy, so maybe this should be defined by default.

It is defined by default; yielding NaNs is a bug. And `sqrt(-1)` gives a `DomainError`, not `NaN`.

Oops. Yeah, this should be fixed then.