6

我不明白如何使用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。我是否错误地使用了枢轴？

4

1 回答 1

5

对于满秩输入，即正定矩阵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。分解一直进行到第一个低于容差的对角线元素，在退出时简单地保持尾随矩阵不变。

感谢 Ian 的以下观察：

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

于 2017-03-30T14:53:19.977 回答