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:

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

enter image description here

# 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.

  • $\begingroup$ What are the actual numerical values that are negative? Is it -100 or -0.00000000000000001? $\endgroup$
    – Sycorax
    Commented Mar 26 at 23:27
  • $\begingroup$ @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. $\endgroup$
    – hannah
    Commented Mar 26 at 23:29
  • $\begingroup$ Marten is a repeated typo for Matérn. I have fixed your text but left your code untouched, $\endgroup$
    – Nick Cox
    Commented Mar 27 at 9:20

1 Answer 1


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.

  • 2
    $\begingroup$ +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. $\endgroup$
    – whuber
    Commented Mar 27 at 15:01
  • $\begingroup$ @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). $\endgroup$
    – hannah
    Commented Mar 27 at 17:09
  • 1
    $\begingroup$ 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? $\endgroup$
    – g g
    Commented Mar 27 at 17:30
  • $\begingroup$ @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. $\endgroup$
    – whuber
    Commented Mar 27 at 19:56
  • $\begingroup$ @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 $\endgroup$
    – hannah
    Commented Mar 27 at 20:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.