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
```