R 中线性回归线的 QR 分解（使用 qr.solve）

Question

所以我试图弄清楚如何使用 R 中的命令

qr.solve

来根据特定数据集求解最佳拟合线（斜率和截距）的分量。尽管我已经知道如何使用

lm

函数来执行此操作，但据我了解，矩阵的 qr 分解可能会产生相同的结果。有人能够解释这样的事情是如何成立的吗？也许它比简单的线性模型命令有好处？

qr.solve

函数的输入可以在 R 中实现吗？还是我必须先自己解决它，然后再将其插入 R？

我尝试将数据集（2 列，每行代表图表上的点）输入到矩阵中，并将其用作函数

qr.solve

中的 x 参数。但是，我不太确定要为 b 插入什么。

Answer 1

QR分解可以求解Ax = b。因此，当最小二乘解达到这一点时：XTy = XTXB，我们可以使用 A = XTX 和 b = XTy 来求解 x = B。

?qr.solve

mtcars
lm(mpg ~ cyl, data=mtcars)
with(mtcars, {
  X <<- cbind(1, cyl)
  y <- mpg
  A = t(X) %*% X
  b = t(X) %*% y
  qr.solve(A, b)
})

Answer 2

稍微详细一点，

线性回归旨在用正规方程求解最小二乘问题：

beta = inv(X'X)(X'y) = solve(X'X, X'y)

这里，

X

表示具有独立预测变量的矩阵，

y

是因响应变量，

beta

是要估计的系数。另外，

X'

表示矩阵

X

的转置。

通过

QR

分解，

X

被分解为正交矩阵

Q

和上三角矩阵

R

，如下所示：

X = QR

，这样

Q'Q = I

然后

beta

可以估计如下：

beta = inv(X'X)(X'y) = inv(R'Q'QR)(R'Q'y) = inv(R'R)(R'Q'y) = inv(R)inv(R')(R')(Q'y)

即，

beta = inv(R)(Q'y) = solve(R, Q'y)

上面显示了

QR

分解如何使线性回归变得更便宜，因为

R

是上三角矩阵，所以可以更有效地求解线性系统。

现在，让我们比较一下求解线性系统的不同方法的性能。一些随机生成的数据：

m <- 100
n <- 6
set.seed(1)
X <- matrix(rnorm(m*n), nrow=m)
beta <- matrix(sample(1:(n+1),n+1), ncol=1)
y <- X1%*% beta + rnorm(m)
X1 <- cbind(rep(1,m), X) # add intercept

首先让我们确认以下所有方法都会产生相同的输出（计算相同的系数）：

lm(y~X)$coeff   
# (Intercept)          X1          X2          X3          X4          X5          X6 
#   6.5233774  -0.2827844   0.3469484  -0.9083762   0.4903771  -0.3217761  -0.0326980 

as.numeric(solve(t(X1)%*%X1, t(X1)%*%y))
# [1]  6.5233774 -0.2827844  0.3469484 -0.9083762  0.4903771 -0.3217761 -0.0326980

q <- qr(X1)
solve(qr.R(q), t(qr.Q(q))%*%y)

as.numeric(solve(qr.R(q), t(qr.Q(q))%*%y))
# [1]  6.5233774 -0.2827844  0.3469484 -0.9083762  0.4903771 -0.3217761 -0.0326980

as.numeric(qr.solve(X1, y))
# [1]  6.5233774 -0.2827844  0.3469484 -0.9083762  0.4903771 -0.3217761 -0.0326980

正如预期的那样，我们在上述方法中获得了完全相同的系数值。

现在，让我们比较一下使用正规方程和使用

QR

分解求解线性系统的运行时间（假设

Q

和

R

已经预先计算）。

library(microbenchmark)
library(ggplot2)
q <- qr(X1)
Q <- qr.Q(q)
R <- qr.R(q)
autoplot(
  microbenchmark(solve(t(X1)%*%X1, t(X1)%*%y),
                 solve(R, t(Q)%*%y), times=1000L)
) 

再次正如预期的那样，通过

QR

分解获得的矩阵可以更快地获得解。