Cholesky Decomposition of a random exponential correlation matrix in R

Question

I have a set of exponential correlation matrices created using the following code.

x=runif(n)
R=matrix(0,n,n)  
for (j in 1:n)
{
  for(k in 1:n)
  {
    R[j,k]=exp(-(x[j]-x[k])^2);
  }
}

and now I want to get their Cholesky decomposition. But many of these are negative definite. How could I resolve this?

Answer 1

The exponential correlation matrix used in spatial or temporal modeling, has a factor alpha that controls the speed of decay:

exp(- alpha * (x[i] - x[j]) ^ 2))

You have fixed such factor at 1. But in practice, such factor is estimated from data.

Note that alpha is necessary to ensure numerical positive definiteness. This matrix is in principle positive definite, but numerically not if alpha is not large enough for a fast decay.

Given that x <- runif(n, 0, 1), the distance between x[i] and x[j] is clustered in a short range [0, 1]. This is not a big range to see a decay in correlation, and maybe you want to try alpha = 10000.

Alternatively if you want to stay with alpha = 1, you need to make distance more spread out. Try x <- runif(n, 0, 100). The decay is very fast, even with alpha = 1.

So we see a duality between distance and alpha. This is also the reason why such correlation matrix can be used stably in statistical modeling. When alpha is to be estimated, it can be made adaptive to distance, so that the correlation matrix is always positive definite.

---

Examples:

f <- function (xi, xj, alpha) exp(- alpha * (xi - xj) ^ 2)

n <- 100

# large alpha, small distance
x <- runif(n, 0, 1)
A <- outer(x, x, f, alpha = 10000)
R <- chol(A)

# small alpha, large distance
x <- runif(n, 0, 100)
A <- outer(x, x, f, alpha = 1)
R <- chol(A)