r语言 - 在正半定矩阵的乔列斯基分解中正确使用枢轴

Question

我不明白如何使用 R 中的 chol 函数来分解正半定矩阵。(或者我这样做，并且有一个错误。文档指出：

如果枢轴 = TRUE，则可以计算正半定 x 的乔列斯基分解。x 的秩以 attr(Q， "rank"( 的形式返回，受数值误差的影响。枢轴返回为 attr(Q， "枢轴"(。t(Q( %*% Q 等于 x 的情况已不再如此。但是，设置 pivot <- attr(Q， "pivot"( 和 oo <- order(pivot(，t(Q[， oo]( %*% Q[， oo] 等于 x ...

下面的例子似乎与这种描述相悖。

> x <- matrix(1, nrow=3, ncol=3)
> Q <- chol(x, pivot=TRUE)
> oo <- order(attr(Q, 'pivot'))
> t(Q[, oo]) %*% Q[, oo]
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    3

结果不是x.我是否错误地使用了枢轴？

Answer 1

对于全秩输入，即正定矩阵x，我们需要

Q <- chol(x, TRUE)
oo <- order(attr(Q, 'pivot'))
unpivQ <- Q[, oo]
all.equal(crossprod(unpivQ), x)

对于有效的秩缺陷输入，即正半定矩阵x(具有负特征值的不定矩阵是非法的，但不在chol选中(，请记住将缺陷尾随对角线块归零：

Q <- chol(x, TRUE)
r <- attr(Q, 'rank')
if (r < nrow(x)) Q[(r+1):nrow(x), (r+1):nrow(x)] <- 0
oo <- order(attr(Q, 'pivot'))
unpivQ <- Q[, oo]
all.equal(crossprod(unpivQ), x)

有些人称之为chol的"错误"，但实际上它是底层LAPACK例程dpstrf的一个功能。因式分解一直持续到第一个低于容差的对角线元素，在退出时保持尾随矩阵不变。

<小时 />

感谢伊恩的以下观察：

您可以在 Q[-(1:r), -(1:r)] <- 0 中使用 R 的负索引来跳过 if 语句。