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