通过QR分解,SVD(和Cholesky分解？)计算投影/帽子矩阵

Question

我试图在R中计算P任意N x J矩阵的投影矩阵S:

P = S (S'S) ^ -1 S'

我一直试图用以下功能执行此操作:

P <- function(S){
  output <- S %*% solve(t(S) %*% S) %*% t(S)
  return(output)
}

但是当我使用它时,我得到的错误看起来像这样:

# Error in solve.default(t(S) %*% S, t(S), tol = 1e-07) : 
#  system is computationally singular: reciprocal condition number = 2.26005e-28

我认为这是数值下溢和/或不稳定的结果在许多地方一样讨论R-帮助和这里,但我没有足够的经验使用SVD或QR分解来解决这个问题,要不然就把这个现有的代码到行动.我也尝试了建议的代码,即将solve写成一个系统:

output <- S %*% solve (t(S) %*% S, t(S), tol=1e-7)

但它仍然无效.任何建议,将不胜感激.

我很确定我的矩阵应该是可逆的并且没有任何共线性,只是因为我尝试用正交虚拟变量矩阵进行测试,但它仍然不起作用.

另外,我想将它应用于相当大的矩阵,所以我正在寻找一个简洁的通用解决方案.

Answer 1

虽然OP已经活动了一年多,但我仍然决定发布答案.我将使用X而不是S像统计中那样,我们经常需要在线性回归上下文中使用投影矩阵,其中X模型矩阵y是响应向量,而H = X(X'X)^{-1}X'帽子/投影矩阵则是Hy预测值.

这个答案假设普通最小二乘的背景.对于加权最小二乘,请参阅从QR分解中获取帽子矩阵以进行加权最小二乘回归.

概述

solve基于一般方阵的LU分解.对于X'X(应该通过crossprod(X)而不是t(X) %*% X在R中计算,读取?crossprod更多)是对称的,我们可以使用chol2inv基于Choleksy分解的.

然而,三角分解比QR分解更不稳定.这不难理解.如果X有条件数kappa,X'X则会有条件数kappa ^ 2.这可能导致很大的数值困难.您收到的错误消息:

# system is computationally singular: reciprocal condition number = 2.26005e-28

只是告诉这个.kappa ^ 2大约是e-28在比周围机器精度远小得多e-16.随着宽容tol = .Machine$double.eps,X'X将被视为排名不足,因此LU和Cholesky分解将被打破.

一般情况下,我们切换到在这种情况下SVD或QR,但摆动 Cholesky分解是另一种选择.

• SVD是最稳定的方法,但太贵了;
• QR具有令人满意的稳定性,计算成本适中,并且在实践中常用;
• 透视Cholesky速度快,稳定性可接受.对于大型基质,这是优选的.

在下文中,我将解释所有三种方法.

使用QR分解

注意,投影矩阵是排列独立的,即,无论是否使用旋转进行QR分解都无关紧要.

在R中,qr.default可以调用LINPACK例程DQRDC进行非旋转QR分解,并使用LAPACK例程DGEQP3进行块旋转QR分解.让我们生成一个玩具矩阵并测试两个选项:

set.seed(0); X <- matrix(rnorm(50), 10, 5)
qr_linpack <- qr.default(X)
qr_lapack <- qr.default(X, LAPACK = TRUE)

str(qr_linpack)
# List of 4
# $ qr   : num [1:10, 1:5] -3.79 -0.0861 0.3509 0.3357 0.1094 ...
# $ rank : int 5
# $ qraux: num [1:5] 1.33 1.37 1.03 1.01 1.15
# $ pivot: int [1:5] 1 2 3 4 5
# - attr(*, "class")= chr "qr"

str(qr_lapack)
# List of 4
# $ qr   : num [1:10, 1:5] -3.79 -0.0646 0.2632 0.2518 0.0821 ...
# $ rank : int 5
# $ qraux: num [1:5] 1.33 1.21 1.56 1.36 1.09
# $ pivot: int [1:5] 1 5 2 4 3
# - attr(*, "useLAPACK")= logi TRUE
# - attr(*, "class")= chr "qr"

请注意$pivot,两个对象不同.

现在,我们定义一个包装函数来计算QQ':

f <- function (QR) {
  ## thin Q-factor
  Q <- qr.qy(QR, diag(1, nrow = nrow(QR$qr), ncol = QR$rank))
  ## QQ'
  tcrossprod(Q)
  }

我们将看到qr_linpack并qr_lapack给出相同的投影矩阵:

H1 <- f(qr_linpack)
H2 <- f(qr_lapack)

mean(abs(H1 - H2))
# [1] 9.530571e-17

使用奇异值分解

在R中,svd计算奇异值分解.我们仍然使用上面的示例矩阵X:

SVD <- svd(X)

str(SVD)
# List of 3
# $ d: num [1:5] 4.321 3.667 2.158 1.904 0.876
# $ u: num [1:10, 1:5] -0.4108 -0.0646 -0.2643 -0.1734 0.1007 ...
# $ v: num [1:5, 1:5] -0.766 0.164 0.176 0.383 -0.457 ...

H3 <- tcrossprod(SVD$u)

mean(abs(H1 - H3))
# [1] 1.311668e-16

同样,我们得到相同的投影矩阵.

使用Pivoted Cholesky分解

为了演示,我们仍然使用X上面的示例.

## pivoted Chol for `X'X`; we want lower triangular factor `L = R'`:
## we also suppress possible rank-deficient warnings (no harm at all!)
L <- t(suppressWarnings(chol(crossprod(X), pivot = TRUE)))

str(L)
# num [1:5, 1:5] 3.79 0.552 -0.82 -1.179 -0.182 ...
# - attr(*, "pivot")= int [1:5] 1 5 2 4 3
# - attr(*, "rank")= int 5

## compute `Q'`
r <- attr(L, "rank")
piv <- attr(L, "pivot")
Qt <- forwardsolve(L, t(X[, piv]), r)

## P = QQ'
H4 <- crossprod(Qt)

## compare
mean(abs(H1 - H4))
# [1] 6.983997e-17

同样,我们得到相同的投影矩阵.