R 的 `chol` 不同于 MATLAB 的 `cholcov`。如何进行类似 Cholesky 的协方差分解？

Question

我一直在尝试在 R 中重现类似 cholesky 的协方差分解 - 就像在 Matlab 中使用 cholcov() 完成的那样。示例取自 https://uk.mathworks.com/help/stats/cholcov.html.

他们示例中原始 cholcov() 函数的结果：

T =
   -0.2113    0.7887   -0.5774         0
    0.7887   -0.2113   -0.5774         0
    1.1547    1.1547    1.1547    1.7321

我正在尝试在 R 中复制这个 T。我试过：

C1 <- cbind(c(2,1,1,2), c(1,2,1,2), c(1,1,2,2), c(2,2,2,3))
T1 <- chol(C1)
C2 <- t(T1) %*% T1

我的结果：

         [,1]      [,2]      [,3]         [,4]
[1,] 1.414214 0.7071068 0.7071068 1.414214e+00
[2,] 0.000000 1.2247449 0.4082483 8.164966e-01
[3,] 0.000000 0.0000000 1.1547005 5.773503e-01
[4,] 0.000000 0.0000000 0.0000000 1.290478e-08

C2恢复了C1，但是T1和MATLAB的解法大不相同。然后我想也许这将是协方差矩阵的 Cholesky 组合：

T1 <- chol(cov(C1))

但我明白了

          [,1]      [,2]      [,3]         [,4]
[1,] 0.5773503 0.0000000 0.0000000 2.886751e-01
[2,] 0.0000000 0.5773503 0.0000000 2.886751e-01
[3,] 0.0000000 0.0000000 0.5773503 2.886751e-01
[4,] 0.0000000 0.0000000 0.0000000 3.725290e-09

这也不对。

谁能告诉我如何在 Matlab 中计算 cholcov()，以便我可以在 R 中复制它？

Answer 1

在这种情况下，您实际上是在滥用 R 函数 chol。来自 MATLAB 的 cholcov 函数是一个复合函数。

如果协方差为正，则进行Cholesky分解，返回满秩上三角Cholesky因子；
如果协方差是半正定的，则进行特征分解，返回一个矩形矩阵。

另一方面，R 中的 chol 只进行 Choleksy 分解。你举的例子C1属于第二种情况。所以，我们应该求助于 R 中的 eigen 函数。

E <- eigen(C1, symmetric = TRUE)
#$values
#[1] 7.000000e+00 1.000000e+00 1.000000e+00 2.975357e-17
#
#$vectors
#           [,1]          [,2]          [,3]       [,4]
#[1,] -0.4364358  0.000000e+00  8.164966e-01 -0.3779645
#[2,] -0.4364358 -7.071068e-01 -4.082483e-01 -0.3779645
#[3,] -0.4364358  7.071068e-01 -4.082483e-01 -0.3779645
#[4,] -0.6546537  8.967707e-16 -2.410452e-16  0.7559289

V <- E$vectors
D <- sqrt(E$values)  ## root eigen values

由于数值秩为 3，我们舍弃最后一个特征值和特征向量：

V1 <- V[, 1:3]
D1 <- D[1:3]

因此你想要的因素是：

R <- D1 * t(V1)  ## diag(D1) %*% t(V1)
#           [,1]       [,2]       [,3]          [,4]
#[1,] -1.1547005 -1.1547005 -1.1547005 -1.732051e+00
#[2,]  0.0000000 -0.7071068  0.7071068  8.967707e-16
#[3,]  0.8164966 -0.4082483 -0.4082483 -2.410452e-16

我们可以验证：

crossprod(R)  ## t(R) %*% R

#     [,1] [,2] [,3] [,4]
#[1,]    2    1    1    2
#[2,]    1    2    1    2
#[3,]    1    1    2    2
#[4,]    2    2    2    3

上面的R因子与cholcov返回的因子不同，因为特征分解使用的算法不同。 R 使用 LAPACK 例程 DSYVER，其中进行了一些旋转，以便本征值不增加。 MATLAB 的 cholcov 不是开源的，所以我不确定它使用的是什么算法。但是很容易证明它没有按非递增顺序排列特征值。

考虑 cholcov 返回的因子 T:

T <- structure(c(-0.2113, 0.7887, 1.1547, 0.7887, -0.2113, 1.1547, 
-0.5774, -0.5774, 1.1547, 0, 0, 1.7321), .Dim = 3:4)

我们可以通过

得到特征值

rowSums(T ^ 2)
# [1] 1.000086 1.000086 7.000167

因为T不精确，所以有一些舍入误差，但我们可以清楚地看到特征值是1, 1, 7。另一方面，我们有来自 R 的 7, 1, 1（回想一下 D1）。

R 的 `chol` 不同于 MATLAB 的 `cholcov`。如何进行类似 Cholesky 的协方差分解？

R's `chol` differs from MATLAB's `cholcov`. How to do a Cholesky-alike covariance decomposition?

matlab

r

matrix

covariance

matrix-factorization