Creating a covariance matrix from a distance matrix—normalizing matrix with infinity

I'm trying to simulate spatially autocorrelated data by creating a covariance matrix from a distance matrix and using said covariance matrix as the $$\Sigma$$ parameter of a multivariate normal.

I have no issues doing this with exponential and linear decay across space, but I'm having trouble with inverse power $$\frac{\delta}{(d_{i,j})^k}$$ for $$\delta, k>0$$ since as $$d_{i,j} \to 0$$, the function goes to $$\infty$$ as opposed to 1.

Suppose you have units located on the x-y plane like so:

set.seed(08192020)
N <- 5
data <- data.frame(lat = runif(N, 0, sqrt(N)), long = runif(N, 0, sqrt(N)))


I can create a distance matrix for them:

distance_matrix <- as.matrix(dist(data, diag=T, upper=T))


Using distance_matrix, I create the first version of $$\Sigma$$:

delta <- sqrt(N)
k <- 2
sigma_inv <- (delta)/(distance_matrix^k)


This yields the following matrix:

╔═══╦═════════════╦═════════════╦═════════════╦═════════════╦═════════════╗
║   ║ 1           ║ 2           ║ 3           ║ 4           ║ 5           ║
║ 1 ║ Inf         ║ 1.373964592 ║ 4.475779735 ║ 1.485955293 ║ 1.61804696  ║
║ 2 ║ 1.373964592 ║ Inf         ║ 1.622451243 ║ 17.37799124 ║ 221.7449699 ║
║ 3 ║ 4.475779735 ║ 1.622451243 ║ Inf         ║ 1.272707732 ║ 1.891004307 ║
║ 4 ║ 1.485955293 ║ 17.37799124 ║ 1.272707732 ║ Inf         ║ 18.05266302 ║
║ 5 ║ 1.61804696  ║ 221.7449699 ║ 1.891004307 ║ 18.05266302 ║ Inf         ║
╚═══╩═════════════╩═════════════╩═════════════╩═════════════╩═════════════╝


Now, for this to be a valid covariance matrix, I need the matrix to be positive semi-definite and the elements along the diagonal to be 1. The infinite values along the diagonal, as well as the values that are > 1 on the off-diagonal are giving me issues.

I've tried:

1. Normalizing the data, e.g. (sigma_inv - min(sigma_inv)) / (max(sigma_inv) - min(sigma_inv)) but that just yields 0s and NaNs.
2. Normalizing the off-diagonal elements—but then the max value of the off-diagonals is 1, when I need the max value of the off-diagonals to not be equal to 1.
3. Changing the values of the diagonal elements to 1 and then normalizing, but then of course some off-diagonal elements will be larger than the diagonal elements.

What's the best way to handle this issue?

• There's no reason to suppose this approach will ever create a valid covariance matrix out of a distance matrix. The functions of distance that will succeed are called correlograms. Because it's generally difficult to determine whether a function is a correlogram, people usually construct them using basic operations applied to a small but fairly rich family of standard correlogram models. These include linear and exponential decay (which is why it has worked for you before). The valid models depend on the dimension of the space, BTW. – whuber Aug 19 at 15:37
• @whuber Thank you! That makes a lot of sense. I'm working just in two-dimensions, so I'll look into the correlograms for that space. – Julian Aug 19 at 15:54