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

General Usage
1

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
2

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
3

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.

4

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
5

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
6

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
7

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: spurious NaN values from matrix square root · Issue #35058 · JuliaLang/julia · GitHub

1 Like
8

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.

9

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

10

Oops. Yeah, this should be fixed then.