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.

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

Given a covariance matrix, show that it is positive semidefinite.

However, I am not sure

What properties should a matrix have to be a covariance matrix.

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.

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

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.

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.

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

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.

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.