# Verify a matrix is positive semi-definite

I am trying to numerically verify that

A symmetric matrix \mathbf{A} is positive semidefinite if and only if it is a covariance matrix.

Then I need to verify in both directions, i.e.

1. Given a positive semidefinite matrix \mathbf{A}, show that it is a covariance matrix.
2. Given a covariance matrix, show that it is positive semidefinite.

However, I am not sure

1. What properties should a matrix have to be a covariance matrix.
2. I know I could generate a covariance matrix using the following and I know that cov is positive semidefinite if and only if all of its eigenvalues are non-negative. But it turns out that minimum(eigvals(cov)) is a negative number close to 0 (on the order \sim 10^{-15}), I am not sure if I could conclude that cov is positive semidefinite since numerical reasons.
n = 100
u = randn(n);
cov = u * u'


Any input will be appreciated.

The issue is that what you are creating is not a covariance matrix. Consider this

julia> n = 100

julia> x = randn(n,3); # 3 random variables with sample size n

julia> demean_x = x .- mean(x,dims=1)

julia> cov_x = demean_x'*demean_x./(n-1) # or simply use cov function from Statistics

julia> minimum(eigvals(cov_x))>=0
true


Julia also has a function isposdef under LinearAlgebra, it checks if the matrix is positive definite.

julia> using LinearAlgebra

julia> isposdef(cov_x)
true

2 Likes

I am not sure I understand the whole context, but this is a well-known property (covariance matrices are psd, and psd matrices are covariance matrices). I don’t know how you could verify a statement about an uncountably infinite set numerically.

It is actually a homework question. I think the point of this question is that professor wants us to get familiar with basic Julia programming.

To verify numerically 1, you may find https://github.com/mcreel/Econometrics/blob/master/Examples/GLS/cholesky.jl useful.

To verify 2, just check the eigenvalues. Your code is missing the step of averaging to compute the sample covariance, and is not keeping the dimensionality fixed.

1 Like

A Cholesky decomposition requires positive definiteness though; semidefiniteness is not enough.

1 Like

Sorry for reviving this topic but I think that my question is quite similar to the title of this one. Let me know if I should have created a new one.

I have a Gram matrix A constructed from a Matérn covariance, so mathematically A is posdef. The thing is that for small such matrices, posdef(A) is true, but for large such matrices, posdef(A) becomes false. I guess that this arises from rounding errors but do you know any way to deal with this numeric issue ?

This is probably due to floating point calculations.
Just add to it scaled identity matrix. Scaling factor should have the value of the minimum eigen value of your matrix (Maybe additional eps()).

1 Like

Oh, that’s the so-called jitter, right ? Let’s say that one obtains a matrix Ã =A + tau^2*I
But then, does a random number generation rand(MvNormal(0,Ã)) have the same properties than rand(MvNormal(0,A)) ?

Well, since A isn’t a PSD it is not a valid covariance matrix, so the questions is not well defined :-).
But generally speaking, yes, it will be quite similar.

Another method is to threshold negative eigen values.
So you calculate the Eigen Decomposition of A , and any negative value you set to 0 and then rebuild the matrix. It might be closer to what you need.

2 Likes

Rebuilding may not be necessary, if factors are needed in the later steps anyway. Which is the case for multivariate normals.

Generally, if the problem permits, it is best not to form A at all, but work with a representation that either precludes impossible values, or makes it very easy to regularize them. I am not familiar with Matérn covariance, but for plain vanilla covariance one could just proceed from the QR decomposition of the data.

1 Like

I come back to say that was a perfect solution to solve my problem, thank you !

1 Like

I am happy to hear I could assist.

Enjoy…

2 Likes