半正定矩阵 Cholesky 分解中主元的正确使用

Question

我不明白如何使用 R 中的 chol 函数对半正定矩阵进行因式分解。（或者我这样做了，但有一个错误。）documentation 状态：

If pivot = TRUE, then the Choleski decomposition of a positive semi-definite x can be computed. The rank of x is returned as attr(Q, "rank"), subject to numerical errors. The pivot is returned as attr(Q, "pivot"). It is no longer the case that t(Q) %*% Q equals x. However, setting pivot <- attr(Q, "pivot") and oo <- order(pivot), it is true that t(Q[, oo]) %*% Q[, oo] equals 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 的 'bug'，但实际上它是底层 LAPACK 例程的一个特征 dpstrf。分解一直进行到低于公差的第一个对角线元素，在退出时留下尾随矩阵不受影响。

感谢 Ian 的以下观察：

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

半正定矩阵 Cholesky 分解中主元的正确使用

Correct use of pivot in Cholesky decomposition of positive semi-definite matrix

r

matrix

matrix-decomposition

matrix-factorization