完成正定相关矩阵

Question

我想complete一个相关矩阵，让它变成正定的，我只关心每个变量与第一个变量的相关性，其余的可以是任何东西。因此，例如，我想complete这个相关矩阵，将

--

填入任何一组数字，以便矩阵变为正定矩阵。

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,] 1.00 0.76 0.77 0.78 0.79 0.81 0.82 0.83 0.84
 [2,] 0.76 1.00  --   --   --   --   --   --   -- 
 [3,] 0.77  --  1.00  --   --   --   --   --   -- 
 [4,] 0.78  --   --  1.00  --   --   --   --   -- 
 [5,] 0.79  --   --   --  1.00  --   --   --   -- 
 [6,] 0.81  --   --   --   --  1.00  --   --   -- 
 [7,] 0.82  --   --   --   --   --  1.00  --   -- 
 [8,] 0.83  --   --   --   --   --   --  1.00  -- 
 [9,] 0.84  --   --   --   --   --   --   --  1.00

像 nearcor 这样的函数用于“修复”正定相关矩阵，但它们优化所有变量，而不是填充一些变量，同时尊重其他变量。或者有像 pos_def_limits 这样的函数用于 while-looping 以填充单个下一个相关值，但同时对许多值执行此操作似乎很棘手。

（此相关矩阵的预期用途，是将其提供给 rnorm_multi 以生成满足它的数据。）

Answer 1

可能有一些优雅的线性代数可以解决这个问题，但蛮力似乎可以解决这个问题。

## set up matrix with required form
m <- matrix(NA, 9, 9)
diag(m) <- 1
m[1,-1] <- m[-1, 1] <- c(76:79, 81:84)/100

mk_matrix <- function(p) {
    ## fill in lower triangle
    m[col(m) > 1 & row(m) > col(m)] <- p
    ## symmetrize
    m[upper.tri(m)] <- t(m)[upper.tri(m)]
    m
}

objfun <- function(p = rnorm(8*7/2)) {
    r <- -1*det(mk_matrix(p))
    ## cat(r, "\n")
    r
}

objfun

的正值表示负定矩阵。

似乎我们有大约 50/50 的机会通过填充随机法线值来获得 pos-def 矩阵：

set.seed(101)
rr <- replicate(1000, objfun())
hist(rr, breaks = 50)

在任何情况下，我们都可以使用

optim()

和

abstol

在我们得到负值（== pos def 结果）时立即停止

set.seed(104) ## chosen to get a neg def starting matrix
optim(par = rnorm(8*7/2), 
      fn = objfun, control=list(abstol = 1e-8))

Answer 2

我想你可以从一个正对角线

M

的下三角开始构造一个正定矩阵

L

，例如，

M = L*L^T

，如果

x^T*M*x = x^T*L*L^T*x = |L^T*x|^2

.一个实现例如：

x>0

我们可以看到

n <- 9
L <- diag(n)
L[-1, 1] <- c(76:79, 81:84) / 100
v <- 1 - L[, 1]^2
for (k in 2:length(v)) {
  L[k, 2:k] <- sqrt((r <- runif(k - 1)) / sum(r) * v[k])
}
M <- tcrossprod(L)

我们可以检查所有特征值是否大于

> M
      [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
 [1,] 1.00 0.7600000 0.7700000 0.7800000 0.7900000 0.8100000 0.8200000
 [2,] 0.76 1.0000000 0.8883572 0.7988221 0.7328434 0.7003986 0.7966054
 [3,] 0.77 0.8883572 1.0000000 0.8972778 0.7937081 0.8500339 0.8200438
 [4,] 0.78 0.7988221 0.8972778 1.0000000 0.9093460 0.8651427 0.8668928
 [5,] 0.79 0.7328434 0.7937081 0.9093460 1.0000000 0.8581035 0.9336535
 [6,] 0.81 0.7003986 0.8500339 0.8651427 0.8581035 1.0000000 0.9105348
 [7,] 0.82 0.7966054 0.8200438 0.8668928 0.9336535 0.9105348 1.0000000
 [8,] 0.83 0.6964236 0.7789939 0.7994304 0.8621129 0.9136291 0.9507678
 [9,] 0.84 0.7344271 0.8091267 0.7860827 0.8387232 0.9097534 0.9200662
           [,8]      [,9]
 [1,] 0.8300000 0.8400000
 [2,] 0.6964236 0.7344271
 [3,] 0.7789939 0.8091267
 [4,] 0.7994304 0.7860827
 [5,] 0.8621129 0.8387232
 [6,] 0.9136291 0.9097534
 [7,] 0.9507678 0.9200662
 [8,] 1.0000000 0.9652907
 [9,] 0.9652907 1.0000000

0

Answer 3

> all(eigen(M)$values > 0)
[1] TRUE