通过QR分解，SVD(和Cholesky分解?)计算投影/帽子矩阵 [英] Compute projection / hat matrix via QR factorization, SVD (and Cholesky factorization?)

问题描述

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

I'm trying to calculate in R a projection matrix P of an arbitrary N x J matrix S:

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

我一直在尝试通过以下功能执行此操作:

I've been trying to perform this with the following function:

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

但是当我使用它时，会出现如下错误:

But when I use this I get errors that look like this:

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

我认为，这是数字下溢和/或不稳定性的结果，正如许多地方在

I think that this is a result of numerical underflow and/or instability as discussed in numerous places like r-help and here, but I'm not experienced enough using SVD or QR decomposition to fix the problem, or else put this existing code into action. I've also tried the suggested code, which is to write solve as a system:

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

但是仍然不起作用.任何建议，将不胜感激.

But still it doesn't work. Any suggestions would be appreciated.

我很确定我的矩阵应该是可逆的，并且没有任何共线性，仅是因为我尝试使用正交虚拟变量的矩阵对其进行测试，但仍然无法正常工作.

I'm pretty sure that my matrix should be invertible and does not have any co-linearities, if only because I have tried testing this with a matrix of orthogonal dummy variables, and it still doesn't work.

此外，我想将其应用于相当大的矩阵，所以我正在寻找一种简洁的通用解决方案.

Also, I'd like to apply this to fairly large matrices, so I'm looking for a neat general solution.

推荐答案

尽管OP已有超过一年没有活跃，但我仍然决定发布答案.我将使用X而不是S，因为在统计中，我们经常希望在线性回归上下文中使用投影矩阵，其中X是模型矩阵，y是响应矢量，而H = X(X'X)^{-1}X'是hat/投影矩阵，以便Hy给出预测值.

Although OP has not been active for more than a year, I still decide to post an answer. I would use X instead of S, as in statistics, we often want projection matrix in linear regression context, where X is the model matrix, y is the response vector, while H = X(X'X)^{-1}X' is hat / projection matrix so that Hy gives predictive values.

此答案假设使用普通最小二乘的上下文.有关加权最小二乘，请参见从QR分解中获取帽子矩阵以进行加权最小二乘回归.

This answer assumes the context of ordinary least squares. For weighted least squares, see Get hat matrix from QR decomposition for weighted least square regression.

概述

solve基于一般方阵的LU分解.对于对称的X'X(应该由R中的crossprod(X)而不是t(X) %*% X计算，更多的内容请阅读?crossprod)，我们可以使用基于Choleksy分解的chol2inv.

solve is based on LU factorization of a general square matrix. For X'X (should be computed by crossprod(X) rather than t(X) %*% X in R, read ?crossprod for more) which is symmetric, we can use chol2inv which is based on Choleksy factorization.

但是，三角分解不如QR分解稳定.这并不难理解.如果X具有条件编号kappa，则X'X将具有条件编号kappa ^ 2.这可能会造成很大的数值困难.您收到的错误消息:

However, triangular factorization is less stable than QR factorization. This is not hard to understand. If X has conditional number kappa, X'X will have conditional number kappa ^ 2. This can cause big numerical difficulty. The error message you get:

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

只是在告诉这个. kappa ^ 2约为e-28，远小于e-16左右的机器精度.在容差为tol = .Machine$double.eps的情况下，X'X将被视为秩不足，因此LU和Cholesky分解将崩溃.

is just telling this. kappa ^ 2 is about e-28, much much smaller than machine precision at around e-16. With tolerance tol = .Machine$double.eps, X'X will be seen as rank deficient, thus LU and Cholesky factorization will break down.

通常，在这种情况下，我们会切换到SVD或QR，但是枢纽性乔尔斯基分解是另一种选择.

Generally, we switch to SVD or QR in this situation, but pivoted Cholesky factorization is another choice.

• SVD是最稳定的方法，但是太昂贵了；
• QR具有令人满意的稳定性，且计算成本适中，并且在实践中经常使用；
• 带枢轴的Cholesky速度快，并且具有可接受的稳定性.对于大型矩阵，此为首选.

下面，我将解释所有三种方法.

In the following, I will explain all three methods.

使用QR因式分解

请注意，投影矩阵是独立于排列的，也就是说，无论是否进行透视都可以执行QR因式分解.

Note that the projection matrix is permutation independent, i.e., it does not matter whether we perform QR factorization with or without pivoting.

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

In R, qr.default can call LINPACK routine DQRDC for non-pivoted QR factorization, and LAPACK routine DGEQP3 for block pivoted QR factorization. Let's generate a toy matrix and test both options:

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对于两个对象是不同的.

Note the $pivot is different for two objects.

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

Now, we define a wrapper function to compute 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给出相同的投影矩阵:

We will see that qr_linpack and qr_lapack give the same projection matrix:

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

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

---

使用奇异值分解

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

In R, svd computes singular value decomposition. We still use the above example matrix 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

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

Again, we get the same projection matrix.

使用透视Cholesky分解

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

For demonstration, we still use the example X above.

## 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

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

Again, we get the same projection matrix.

这篇关于通过QR分解，SVD(和Cholesky分解?)计算投影/帽子矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！