利用QR分解、SVD (和Cholesky分解？)计算投影/ hat矩阵。

Question

我试图在R中计算任意N矩阵P的投影矩阵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分解来解决问题，或者将现有的代码付诸实施。我还尝试了所建议的代码，即将解决方案编写为一个系统：

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

但还是不起作用。如有任何建议，将不胜感激。

我很肯定，我的矩阵应该是可逆的，没有任何协线性，如果我试过用一个正交的虚拟变量矩阵来测试它，它仍然不起作用。

另外，我想把它应用到相当大的矩阵中，所以我正在寻找一个简洁的通用解决方案。

Answer 1

虽然OP已经一年多没有激活了，但我仍然决定发布一个答复。我将使用X代替S，因为在统计学中，我们通常需要线性回归上下文中的投影矩阵，其中X是模型矩阵，y是响应向量，而H = X(X'X)^{-1}X'是hat /投影矩阵，因此Hy给出了预测值。

这个答案假定为普通最小二乘的上下文。有关加权最小二乘，请参见加权最小二乘回归QR分解的hat矩阵。

概述

solve是基于一般平方矩阵的LU分解。对于对称的X'X (应该由crossprod(X)而不是t(X) %*% X在R中计算)，可以使用基于Choleksy分解的chol2inv。

然而，三角分解比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

同样，我们得到了同样的投影矩阵。

使用旋转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

同样，我们得到了同样的投影矩阵。