R中随机指数相关矩阵的Cholesky分解答案

【问题标题】：Cholesky Decomposition of a random exponential correlation matrix in RR中随机指数相关矩阵的Cholesky分解
【发布时间】：2017-10-24 10:36:53
【问题描述】：

我有一组使用以下代码创建的指数相关矩阵。

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);
  }
}

现在我想得到他们的 Cholesky 分解。但其中许多是负确定的。我该如何解决这个问题？

【问题讨论】：

欢迎来到 SO！请使用set.seed(..) 作为可复制的示例。您可以通过crossprod(R) 生成半正定矩阵。如果 det(R) > 0 则有偶数个负特征值。

标签： r matrix statistics correlation matrix-factorization

【解决方案1】：

空间或时间建模中使用的指数相关矩阵，有一个因子alpha 控制衰减速度：

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

您已将该因子固定为 1。但实际上，该因子是根据数据估算的。

请注意，alpha 是确保数值正定性所必需的。这个矩阵原则上是正定的，但如果alpha 不足以快速衰减，则数值上不是。

鉴于x <- runif(n, 0, 1)，x[i] 和x[j] 之间的距离聚集在一个短距离[0, 1] 内。看到相关性衰减的范围不是很大，也许您想尝试alpha = 10000。

或者，如果你想留在alpha = 1，你需要让距离更分散。试试x <- runif(n, 0, 100)。衰减非常快，即使使用alpha = 1。

所以我们看到距离和alpha 之间存在对偶性。这也是为什么这种相关矩阵可以稳定地用于统计建模的原因。在估计alpha时，可以使其适应距离，使相关矩阵始终为正定。

例子：

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)

【讨论】：

【解决方案2】：

尝试用它来构造正定矩阵

A<-matrix(runif(n^2),n,n)
dim(A)
A<-A%*%t(A)
chol(A)

【讨论】：