Yesterday I asked a question about why my randomly generated distance matrices were leading to Matérn covariance matrices that were not positive definite. The answer there called my attention to the fact that my simplistic way of generating randomly generated distance matrices was not ensuring that such matrices obeyed the triangle inequality. Indeed, enforcing that solved the problem.

However, what I realize now is that my real data is also leading to that same problem. That is, when calculating a distance matrix for the distances between rows in my dataset, I end up getting a distance matrix that... does not obey triangle inequality. And thus, leads to a Matérn covariance matrix that is not positive definite.

Here is a small subset (1000 rows) of my actual data (about 40,000 rows), which does reproduce the issue, uploaded to Filebin:


Here is the illustration of what I am doing:

library(Rfast) # for fast calculation of distances and checking symmetry
library(fossil) # for checking triangle inequality

# Load subset of my actual data and calculate the distance matrix:
X <- read.csv("actualdata.csv")
d <- Rfast:Dist(X)
d[d==0] <- 0.000000000001 # instead of zero to avoid numerical issues

# 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

# Check if covariance matrix is squared:
nrow(m) == ncol(m)

# Check eigenvalues of Marten covariance matrix:
g <- eigen(m, only.values=TRUE, symmetric = TRUE)
print(min(g$values)) # should be greter than zero, but isn't
print(sum(g$values<0)/length(g$values)) # thus, this should be zero, but isn't

Here is a set of checks of whether the d obeys all 4 properties of distance matrices:

# check that we have a correct distance matrix (warning: the triangle inequality checks may take a while):
Rfast::is.symmetric(d) # symmetric matrix? YES
sum(diag(d)!=0.000000000001)==0 # zero diagonal values? YES
min(d[row(d) == (col(d) - 1)])>0 # positive off-diagonal values? YES
fossil::tri.ineq(d) # obeys triangular inequality? NO (warning: this line is quite slow)

The problem is: because here the input dataset X is not simulated, I cannot force the distance matrix to obey triangle inequality by design. Therefore my question: in this case, what can be done to ensure the generation of a positive definite Matérn covariance matrix?

  • 2
    $\begingroup$ I just ran your code (except for d[d==0] <- 0.000000000001, which has nothing to do with computing a distance matrix). Testing non-negativity reveals an issue: all(upper.tri(d) > 0.0) returns FALSE, so there's at least one non-positive element in the upper triangular part of d. Moreover, examining dim(unique(X)) reveals that there are 16 duplicated rows. $\endgroup$
    – Sycorax
    Commented Mar 28 at 4:28
  • 1
    $\begingroup$ Matern covariances apply to Euclidean distances. Thus, to create a distance matrix, (a) create a set of points in a Euclidean space and (b) compute their distances. Here is an example in dimension d using n points uniformly distributed in $[0,1]^d:$ d <- 2; n <- 10; p <- matrix(runif(d*n), n); m <- as.matrix(dist(p)) Besides being correct, it's even simpler than what you're doing. $\endgroup$
    – whuber
    Commented Mar 28 at 14:26
  • $\begingroup$ @Sycorax thanks for your insights. But your test is incorrect, since upper.tri returns TRUE/FALSE. It should have been all(d[upper.tri(d)] > 0.0), not all(upper.tri(d) > 0.0). Using the correct version returns TRUE. So there are none non-positive values in d. You are correct, however, that X has duplicate entries (as the full dataset also has). I tried adding a small random value to those entries, or even to the entire dataset, to force uniqueness, but that didn't help. $\endgroup$
    – hannah
    Commented Mar 28 at 17:31
  • 1
    $\begingroup$ No -- you are non-randomly altering the distances! Get rid of the 0.000000000001 kludge and please provide a reproducible example of the problem you encounter. $\endgroup$
    – whuber
    Commented Mar 28 at 17:42
  • 2
    $\begingroup$ @hannah d[d==0] <- 0.000000000001 is a kludge. If you truly need the duplicated data for some reason then you'll need to change the covariance matrix, not the distances. The answer stats.stackexchange.com/a/643743/22311 explains this in more detail: replace the removable discontinuity at 0. Using a kernel matrix of the form $K + \epsilon I$ is a standard approach for mitigating the numerical issues of identical or nearly-identical data. $\endgroup$
    – Sycorax
    Commented Mar 28 at 18:38

1 Answer 1

  1. Don't fudge the diagonal. Replace the final diagonal in m by sigma**2 to get rid of the (removable) singularity from 0 * Inf in the original expression.
  2. Your parameter v=3 is unusual. Are you sure about this? Matérn-5/2 is with v=5/2.
  3. Your data may be difficult there seem to be many identical coordinates, i.e. points are very close.
  4. Scaling by p=5 decreases distances further.

The items above mean you will run into numerical difficulties for a large data set. Test this by starting with just a few points and increase to see when neg Eigenvalues will appear.

If you really need to fudge, add a small diagonal term. I.e. instead of using m as covariance use $m + \epsilon I$ for small $\epsilon>0$.

I have no idea what tri.ineq is doing. Most likely it runs into numerical difficulties as well. Don't bother with this test too much. A distance matrix is a distance matrix (up to pesky numerical issues).

You may also want to check issues against the simpler expressions for specific choices of v. See Wikipedia In particular, the covariance with v<-1/2 should be pretty robust.

  • 2
    $\begingroup$ Ah and a final remark: Don't call it Marten. $\endgroup$
    – g g
    Commented Mar 28 at 9:50
  • $\begingroup$ Fantastic, this was really helpful. Playing with the parameters proved to be enough, but the other insights were also quite helpful. $\endgroup$
    – hannah
    Commented Mar 28 at 18:46

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.