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

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.

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

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.