在正半确定矩阵的Cholesky分解中正确使用枢轴

Question

我不明白如何使用cholR中的函数来计算正半正定矩阵.(或者我这样做,并且有一个错误.)文档说明:

如果pivot = TRUE,则可以计算正半正定x的Choleski分解.x的等级作为attr(Q,"rank")返回,受数字误差的影响.枢轴以attr(Q,"pivot")返回.不再是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)

有些人称之为'bug' chol,但实际上它是底层LAPACK例程的一个特性dpstrf.分解进行到第一对角线元素,其低于公差,使得尾随矩阵在退出时简单地不受影响.

感谢Ian对以下观察:

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