在正半定矩阵的Cholesky分解中正确使用枢轴 [英] Correct use of pivot in Cholesky decomposition of positive semi-definite matrix
问题描述
我不明白如何在R中使用chol
函数来分解一个正半定矩阵. (或者我这样做,并且有一个错误.)
I don't understand how to use the chol
function in R to factor a positive semi-definite matrix. (Or I do, and there's a bug.) The documentation states:
如果ivot = TRUE,则可以计算正半定x的Choleski分解. x的等级以attr(Q,"rank")的形式返回,这会受到数值误差的影响.枢轴返回为attr(Q,"pivot"). t(Q)%*%Q不再等于x.但是,设置枢轴<-attr(Q,"pivot")和oo <-order(pivot),则确实是t(Q [,oo])%*%Q [,oo]等于x ...
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 ...
以下示例似乎掩盖了这一描述.
The following example seems to belie this description.
> 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
.我是否使用了不正确的枢轴?
The result is not x
. Am I using the pivot incorrectly?
推荐答案
对于全秩输入,即正定矩阵x
,我们需要
For full rank input, i.e., a positive-definite matrix x
, we need
Q <- chol(x, TRUE)
oo <- order(attr(Q, 'pivot'))
unpivQ <- Q[, oo]
all.equal(crossprod(unpivQ), x)
对于有效秩不足的输入,即正半定矩阵x
(具有负本征值的不确定矩阵是非法的,但未在chol
中进行检查),请记住零缺陷尾随对角线块:
For a valid rank-deficient input, i.e., a positive semi-definite matrix x
(indefinite matrix with negative eigen values are illegal, but not checked in chol
), remember to zero deficient trailing diagonal block:
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
的功能.分解继续进行,直到第一个对角线元素低于公差为止,而尾随矩阵在出口处保持不变.
Some people call this a 'bug' of chol
, but actually it is a feature of the underlying LAPACK routine dpstrf
. The factorization proceeds till the first diagonal element which is below a tolerance, leaving the trailing matrix simply untouched on exit.
感谢伊恩(Ian)的以下观察:
Thanks to Ian for the following observation:
您可以在Q[-(1:r): -(1:r)] <- 0
中使用R的负索引来跳过if
语句.
You could use R's negative indexing in Q[-(1:r): -(1:r)] <- 0
to skip the if
statement.
这篇关于在正半定矩阵的Cholesky分解中正确使用枢轴的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!