# Calculating Matérn covariance in R returns matrices that are not positive definite

In R, I am trying to calculate Matérn covariance matrices whose inputs are randomly created distance matrices. However, I often end up getting covariance matrices that are not positive definite, which makes little sense.

At first, I had asked a version of this questions at StackOverflow, but a few comments made there reassured me that the code was not at fault. Additionally, I was told that while the Matérn covariance matrix obviously is, in theory, always positive definite, in practice it may not be.

At that point, I realized that this became a question more suited to CrossValidated. So, what I am trying to understand is how to ensure the generation of a positive definite Matérn covariance matrix in practice.

Perhaps there are specific relations between the Matérn equation parameters that ensure that? Perhaps there are any other tricks - for example, more general questions about ensuring positive definiteness in practice are usually answered by simply adding a tiny constant to the otherwise zero diagonals of the distance matrices - which does not solve my problem here.

For reference, here's how I generate a random distance matrix:

nrows <- 100
ncols <- 100
d <- matrix(runif(nrows*ncols, 0, 1), ncols, nrows)
# enforce symmetry:
d[lower.tri(d)] <- t(d)[lower.tri(d)]
diag(d) <- 0.000000000001 # instead of zero to avoid numerical issues
# check that we have a correct distance matrix:
range(d)


Then I follow this Wikipedia formula for Matérn covariances:

# Calculate Marten covariance matrix:
sigma <- 1
v <- 3
p <- 5
term1 <- (2**(1-v))/(gamma(v))
term2 <- (sqrt(2*v)*(abs(d)/p))**v
term3 <- besselK(sqrt(2*v)*(abs(d)/p), nu = v)
m <- (sigma**2)*term1*term2*term3


However, when I check the eigenvalues of the Matérn covariance matrix, I see that it's not a positive definite matrix, as it should be:

g <- eigen(m, only.values=TRUE)
print(min(g$$values)) # should be greater than zero, but instead is approx. -0.15 print(sum(g$$values<0)/length(g\$values)) # thus, this should be zero but instead is approx. 0.46


It surprises me that it is so easy to arrive at a non-positive definite Matérn covariance matrix, and I am then puzzled about how one may ensure positive definiteness.

• What are the actual numerical values that are negative? Is it -100 or -0.00000000000000001?
– Sycorax
Commented Mar 26 at 23:27
• @Sycorax running the code snippets above as-is shows that around 47% of the eigen values are negative each time, with the lowest eigen values usually around -0.15. Commented Mar 26 at 23:29
• Marten is a repeated typo for Matérn. I have fixed your text but left your code untouched, Commented Mar 27 at 9:20

Your problem is not in the Matérn part. It is already in the "distance" part, i.e. with your matrix $$d$$. In your code, $$d$$ is just a (slightly fudged) symmetric random matrix with zeros on the diagonal and non-negative entries. But this is not enough to qualify the matrix as a distance matrix. In addition to the three properties (symmetric, zero diagonal and non-negative) for distances, you have the crucial triangle inequality. This is missing from your code. See the relevant section on Wikipedia's article on distance matrices.

• +1. But the code has numerous errors anyway. m is not even a square matrix; when it is converted to one, floating point imprecision will likely create some tiny imaginary components in the eigenvalues (causing an error when calling min); the calculation of m for larger distances might be numerically unstable; and the eigenvalues should be checked either with svd or using symmetric = TRUE in the call to eigen.
– whuber
Commented Mar 27 at 15:01
• @whuber I am really puzzled by your comment. In my code above m is a square matrix. And using eigen(m, symmetric = TRUE) changes nothing (as should be the case, since internally eigen() checks for symmetry anyways, hence making the parameter performance-oriented only). Commented Mar 27 at 17:09
• This seems unlikely, why should your iid generated entries satisfy such a relation? Anyway, why don't you just sample some points, calculate their distances using dist and use this matrix for a test?
– g g
Commented Mar 27 at 17:30
• @hannah I ran your code before making that comment and tested my recommendations. I warmly suggest you do the same! Remember, you need to create an actual distance matrix... . Use the dist function to do that.
– whuber
Commented Mar 27 at 19:56
• @whuber thanks for doing so, and indeed g g's answer proved to be correct; the problem was in d disobeying triangle inequality, which I solved by using the dist function over vectors of data, as you both recommended. That said, for the sake of clarity, what puzzles me is your comment about m not being even a square matrix. I don't know what you tried exactly, but I did run my code, of course, and a check ncol(m) == nrow(m) shows the code does generate a square m. Also, I did try your eigen(m, symmetric = TRUE), which does nothing. See: onecompiler.com/r/428j7yjmp Commented Mar 27 at 20:16