# Matrix square root of known positive semi-definite matrix does not exist?
**URL:** https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766
**Category:** General Usage
**Tags:** linearalgebra
**Created:** [March 9, 2020, 7:17pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766 "2020-03-09T19:17:03Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@torgo](https://discourse.julialang.org/u/torgo)
#### Post date: [March 9, 2020, 7:17pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/1 "2020-03-09T19:17:03Z")
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](https://en.wikipedia.org/wiki/Definiteness_of_a_matrix#Square_root).
However, Julia returns a matrix of complex `NaN`’s.
Can anyone help me out with what is going on here?
```julia
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:
```julia
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
```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [March 9, 2020, 7:25pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/2 "2020-03-09T19:25:43Z")
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.
---
### Author:  [@torgo](https://discourse.julialang.org/u/torgo)
#### Post date: [March 9, 2020, 7:36pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/3 "2020-03-09T19:36:54Z")
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.
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [March 9, 2020, 8:13pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/4 "2020-03-09T20:13:04Z")
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.
---
### Author:  [@torgo](https://discourse.julialang.org/u/torgo)
#### Post date: [March 9, 2020, 8:22pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/5 "2020-03-09T20:22:46Z")
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
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
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 9, 2020, 9:06pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/6 "2020-03-09T21:06:02Z")
> [@torgo](#):
>
> `C = A * B * A'`
As noted by others above, roundoff errors give you tiny negative eigenvalues. One thing you can do is a slight shift:
```julia
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](https://github.com/JuliaLang/julia/pull/35057)
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 9, 2020, 9:24pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/7 "2020-03-09T21:24:02Z")
> [@Oscar\_Smith](#):
>
> As such, `sqrt(C)` does not exist.
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](https://github.com/JuliaLang/julia/issues/35058)
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [March 9, 2020, 10:03pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/8 "2020-03-09T22:03:56Z")
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.
---
### Author:  [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)
#### Post date: [March 9, 2020, 10:18pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/9 "2020-03-09T22:18:39Z")
It is defined by default; yielding NaNs is a bug. And `sqrt(-1)` gives a `DomainError`, not `NaN`.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [March 9, 2020, 10:24pm UTC](https://discourse.julialang.org/t/matrix-square-root-of-known-positive-semi-definite-matrix-does-not-exist/35766/10 "2020-03-09T22:24:14Z")
Oops. Yeah, this should be fixed then.