附录 B — 矩阵运算

$$$$

There’s probably some examples, but there are some examples of people using solve(t(X) %*% W %*% X) %*% W %*% Y to compute regression coefficients, too.

— Thomas Lumley ¹

本文主要介绍 Base R 提供的矩阵运算，包括加、减、乘等基础矩阵运算和常用的矩阵分解方法，总结 Base R 、Matrix 包和 Eigen 库对应的矩阵运算函数，分别对应基础、进阶和高阶的读者。最后，介绍矩阵运算在线性回归中的应用。

library(Matrix)

B.1 基础运算

约定符号

$$A=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$$

B.1.1 加、减、乘

矩阵 $A$

A <- matrix(c(1, 1.2, 1.2, 3), nrow = 2)
A

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

B <- matrix(c(1, 2, 3, 4), nrow =2)
B

     [,1] [,2]
[1,]    1    3
[2,]    2    4

A + A # 对应元素相加

     [,1] [,2]
[1,]  2.0  2.4
[2,]  2.4  6.0

A - A # 对应元素相减

     [,1] [,2]
[1,]    0    0
[2,]    0    0

A %*% A # 矩阵乘法

     [,1]  [,2]
[1,] 2.44  4.80
[2,] 4.80 10.44

B.1.2 对数、指数与幂

矩阵 $A$ 的对数 $\log A$ ，就是找一个矩阵 $L$ 使得 $A={\mathrm{e}}^L$

expm::logm(A)

           [,1]      [,2]
[1,] -0.4485660 0.8050907
[2,]  0.8050907 0.8932519

矩阵 $A$ 的指数 ${\mathrm{e}}^A$ 的定义

$${\mathrm{e}}^A=\sum _{k=1}^{\mathrm{\infty }}\frac{A^k}{k!}$$

expm 包可以计算矩阵的指数、开方、对数等。

expm::expm(A)

         [,1]     [,2]
[1,]  7.60987 12.93908
[2,] 12.93908 29.17501

或者使用奇异值分解 $A=UD{V}^{\mathrm{\top }}$ ，则 ${\mathrm{e}}^A=U{\mathrm{e}}^D{V}^{\mathrm{\top }}$ ，其中，D 是对角矩阵。

(res <- svd(A))

$d
[1] 3.5620499 0.4379501

$u
           [,1]       [,2]
[1,] -0.4241554 -0.9055894
[2,] -0.9055894  0.4241554

$v
           [,1]       [,2]
[1,] -0.4241554 -0.9055894
[2,] -0.9055894  0.4241554

res$u %*% diag(exp(res$d)) %*% res$v

         [,1]     [,2]
[1,]  7.60987 12.93908
[2,] 12.93908 29.17501

矩阵 $A$ 的 $n$ 次幂 $A^n$ ，利用奇异值分解 $A=UD{V}^{\mathrm{\top }}$

$$\begin{array}{rl}{A}^n& =A\times A\times \cdots \times A\\ & =UD{V}^{\mathrm{\top }}UD{V}^{\mathrm{\top }}\cdots UD{V}^{\mathrm{\top }}\end{array}$$

计算 $A^3$

res$u %*% (diag(res$d)^3) %*% res$v

       [,1]   [,2]
[1,]  8.200 17.328
[2,] 17.328 37.080

B.1.3 迹、秩、条件数

矩阵 $A$ 的迹 $\mathrm{tr}(A)=\sum _{i=1}^n{a}_{ii}$

sum(diag(A))

[1] 4

qr(A)$rank

[1] 2

kappa(A)

[1] 10.41469

B.1.4 求逆与广义逆

Moore-Penrose Generalized Inverse 摩尔广义逆 $A^{-}$。

$$A^{-}=(A^{\mathrm{\top }}A{)}^{-1}A$$

如果 A 可逆，则广义逆就是逆。

solve(A) # 逆

           [,1]       [,2]
[1,]  1.9230769 -0.7692308
[2,] -0.7692308  0.6410256

MASS::ginv(A) # 广义逆

           [,1]       [,2]
[1,]  1.9230769 -0.7692308
[2,] -0.7692308  0.6410256

B.1.5 行列式与伴随

矩阵必须是方阵

伴随矩阵 $A\ast A^{\star }=A^{\star }\ast A=|A|\ast I,A^{\star }=|A|\ast A^{-1}$

• $|A^{\star }|=|A{|}^{n-1},A\in {\mathbb{R}}^{n\times n},n\ge 2$
• $(A^{\star }{)}^{\star }=|A{|}^{n-2}A,A\in {\mathbb{R}}^{n\times n},n\ge 2$
• $(A^{\star }{)}^{\star }$ A 的 n 次伴随是？

det(A)

[1] 1.56

det(A) * solve(A)

     [,1] [,2]
[1,]  3.0 -1.2
[2,] -1.2  1.0

B.1.6 外积、直积与交叉积

通常的矩阵乘法也叫矩阵内积

A %*% B

     [,1] [,2]
[1,]  3.4  7.8
[2,]  7.2 15.6

外积

A %o% B # outer(A, B, FUN = "*")

, , 1, 1

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

, , 2, 1

     [,1] [,2]
[1,]  2.0  2.4
[2,]  2.4  6.0

, , 1, 2

     [,1] [,2]
[1,]  3.0  3.6
[2,]  3.6  9.0

, , 2, 2

     [,1] [,2]
[1,]  4.0  4.8
[2,]  4.8 12.0

直积/克罗内克积

A %x% B # kronecker(A, B, FUN = "*")

     [,1] [,2] [,3] [,4]
[1,]  1.0  3.0  1.2  3.6
[2,]  2.0  4.0  2.4  4.8
[3,]  1.2  3.6  3.0  9.0
[4,]  2.4  4.8  6.0 12.0

交叉积 $A^{\mathrm{\top }}A$

crossprod(A, A)  #  t(x) %*% y

     [,1]  [,2]
[1,] 2.44  4.80
[2,] 4.80 10.44

tcrossprod(A, A) #  x %*% t(y)

     [,1]  [,2]
[1,] 2.44  4.80
[2,] 4.80 10.44

B.1.7 Hadamard 积

Hadamard 积（法国数学家 Jacques Hadamard）也叫 Schur 积（德国数学家 Issai Schur ）或 entrywise 积是两个维数相同的矩阵对应元素相乘，特别地，$A^2$ 表示将矩阵 $A$ 的每个元素平方

$$(A\circ B{)}_{ij}=(A{)}_{ij}(B{)}_{ij}$$

$$\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\circ \left[\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{21}& b_{22}& b_{23}\\ b_{31}& b_{32}& b_{33}\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}{b}_{11}& a_{12}{b}_{12}& a_{13}{b}_{13}\\ a_{21}{b}_{21}& a_{22}{b}_{22}& a_{23}{b}_{23}\\ a_{31}{b}_{31}& a_{32}{b}_{32}& a_{33}{b}_{33}\end{array}\right]$$

fastmatrix::hadamard(A, B)

     [,1] [,2]
[1,]  1.0  3.6
[2,]  2.4 12.0

A^2     # 每个元素平方 a_ij ^ 2

     [,1] [,2]
[1,] 1.00 1.44
[2,] 1.44 9.00

A ** A  # 每个元素的幂 a_ij ^ a_ij

         [,1]      [,2]
[1,] 1.000000  1.244565
[2,] 1.244565 27.000000

2^A     # 每个元素的指数 2 ^ a_ij

         [,1]     [,2]
[1,] 2.000000 2.297397
[2,] 2.297397 8.000000

exp(A)  # 每个元素的指数 exp(a_ij)

         [,1]      [,2]
[1,] 2.718282  3.320117
[2,] 3.320117 20.085537

B.1.8 矩阵范数

矩阵的范数，包括 1，2，无穷范数

$1$-范数

列和绝对值最大的

$2$ - 范数

又称谱范数，矩阵最大的奇异值，如果是方阵，就是最大的特征值

$\mathrm{\infty }$ - 范数

行和绝对值最大的

Frobenius - 范数

Euclidean 范数

$M$ - 范数

矩阵里模最大的元素，矩阵里面的元素可能含有复数，所以取模最大

norm(A, type = "1") # max(abs(colSums(A)))

[1] 4.2

norm(A, type = "I") # max(abs(rowSums(A)))

[1] 4.2

norm(A, type = "F")

[1] 3.588872

norm(A, type = "M") #

[1] 3

norm(A, type = "2") # max(svd(A)$d)

[1] 3.56205

B.1.9 转置与旋转

矩阵 $A$

t(A) # 转置

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

B.1.10 正交与投影

矩阵 $A$ 的投影

$$I-A(A^{\mathrm{\top }}A{)}^{-1}{A}^{\mathrm{\top }}$$

diag(rep(1, 2)) - A %*% solve(t(A) %*% A) %*% t(A)

              [,1]          [,2]
[1,] -6.661338e-16  3.330669e-16
[2,]  3.837386e-16 -2.220446e-16

B.1.11 Givens 变换(*)

• Givens 旋转
• 帽子矩阵在统计中的应用，回归与方差分析 (Hoaglin 和 Welsch 1978)

B.1.12 Householder 变换(*)

Householder 变换是平面反射的一般情况： 要计算 $N\times P$ 维矩阵 $X$ 的 QR 分解，我们采用 Householder 变换

$${\mathbf{H}}_u=\mathbf{I}-2\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}$$

其中 $I$ 是 $N\times N$ 维的单位矩阵，$u$ 是 $N$ 维单位向量，即 $\parallel \mathbf{u}\parallel =\sqrt{\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}}=1$。则 $H_u$ 是对称正交的，因为

$${\mathbf{H}}_u^{\mathrm{\top }}={\mathbf{I}}^{\mathrm{\top }}-2\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}={\mathbf{H}}_u$$

并且

$${\mathbf{H}}_u^{\mathrm{\top }}{\mathbf{H}}_u=\mathbf{I}-4\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}+4\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}=\mathbf{I}$$

让 ${\mathbf{H}}_u$ 乘以向量 $\mathbf{y}$，即

$${\mathbf{H}}_u\mathbf{y}=\mathbf{y}-2\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}\mathbf{y}$$

它是 $y$ 关于垂直于过原点的 $u$ 的直线的反射，只要

$$\begin{array}{}\text{(B.1)}& \begin{array}{r}\mathbf{u}=\frac{\mathbf{y}-\parallel \mathbf{y}\parallel {\mathbf{e}}_1}{\parallel \mathbf{y}-\parallel \mathbf{y}\parallel {\mathbf{e}}_1\parallel }\end{array}\end{array}$$

或者

$$\begin{array}{}\text{(B.2)}& \begin{array}{r}\mathbf{u}=\frac{\mathbf{y}+\parallel \mathbf{y}\parallel {\mathbf{e}}_1}{\parallel \mathbf{y}+\parallel \mathbf{y}\parallel {\mathbf{e}}_1\parallel }\end{array}\end{array}$$

其中 ${\mathbf{e}}_1=(1,0,\dots ,0{)}^{\mathrm{\top }}$，Householder 变换使得向量 $y$ 成为 $x$ 轴，在新的坐标系统中，向量 $H_{u}y$ 的坐标为 $(\pm \parallel y\parallel ,0,\dots ,0{)}^{\mathrm{\top }}$

举个例子

借助 Householder 变换做 QR 分解的优势：

1. 更快、数值更稳定比直接构造 Q，特别当 N 大于 P 的时候
2. 相比于存储矩阵 Q 的 $N^2$ 个元素，Householder 变换只存储 P 个向量 $u_1,\dots ,u_P$
3. QR 分解的真实实现，比如在 LINPACK 中，定义 $u$ 的时候， 式 B.1 或 式 B.2 的选择基于 $y$ 的第一个坐标的符号。如果坐标是负的，使用 式 B.1 ，如果是正的，使用 式 B.2 ， 这个做法可以使得数值计算更加稳定。

用 Householder 变换做 QR 分解 (Bates 和 Watts 1988) 及其 R 语言、Eigen 实现。

B.1.13 单位矩阵

矩阵对角线上全是1，其余位置都是0

$$A=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$$

diag(rep(3))

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

而全1矩阵是所有元素都是1的矩阵，可以借助外积运算构造，如3阶全1矩阵

rep(1,3) %o% rep(1,3) 

     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1

B.1.14 对角矩阵

diag(A)       # 矩阵的对角

[1] 1 3

diag(x = c(1, 2, 3)) # 构造对角矩阵

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    2    0
[3,]    0    0    3

B.1.15 稀疏矩阵

稀疏矩阵的典型构造方式是通过三元组。

i <- c(1, 3:8) # 行指标
j <- c(2, 9, 6:10) # 列指标
x <- 7 * (1:7) # 数据
Matrix::sparseMatrix(i, j, x = x)

8 x 10 sparse Matrix of class "dgCMatrix"
                             
[1,] . 7 . . .  .  .  .  .  .
[2,] . . . . .  .  .  .  .  .
[3,] . . . . .  .  .  . 14  .
[4,] . . . . . 21  .  .  .  .
[5,] . . . . .  . 28  .  .  .
[6,] . . . . .  .  . 35  .  .
[7,] . . . . .  .  .  . 42  .
[8,] . . . . .  .  .  .  . 49

B.1.16 上、下三角矩阵

m <- A
m

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

upper.tri(m) # 矩阵上三角

      [,1]  [,2]
[1,] FALSE  TRUE
[2,] FALSE FALSE

m[upper.tri(m)]

[1] 1.2

m[lower.tri(m)] <- 0 # 获得上三角矩阵
m

     [,1] [,2]
[1,]    1  1.2
[2,]    0  3.0

矩阵 A 的下三角矩阵

m <- matrix(c(1, 2, 2, 3), nrow = 2)
m[row(m) < col(m)] <- 0
m

     [,1] [,2]
[1,]    1    0
[2,]    2    3

B.2 矩阵分解

B.2.1 LU 分解

矩阵 $A$ 的 LU 分解 $A=LU$ ， $L$ 是下三角矩阵，$U$ 是上三角矩阵

Matrix::lu(A)

LU factorization of Formal class 'denseLU' [package "Matrix"] with 4 slots
  ..@ x       : num [1:4] 1.2 0.833 3 -1.3
  ..@ perm    : int [1:2] 2 2
  ..@ Dim     : int [1:2] 2 2
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL

B.2.2 Schur 分解

矩阵 $A$ 的 Schur 分解 $A=QT{Q}^{\mathrm{\top }}$

(res <- Matrix::Schur(A))

$Q
           [,1]       [,2]
[1,] -0.9055894 -0.4241554
[2,]  0.4241554 -0.9055894

$T
          [,1]    [,2]
[1,] 0.4379501 0.00000
[2,] 0.0000000 3.56205

$EValues
[1] 0.4379501 3.5620499

其中 $Q$ 是一个正交矩阵 $QQ=I$ ，$T$ 是一个分块上三角矩阵

res$Q %*% t(res$Q)

             [,1]         [,2]
[1,] 1.000000e+00 8.052102e-18
[2,] 8.052102e-18 1.000000e+00

res$Q %*% res$T %*% t(res$Q)

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

B.2.3 QR 分解

矩阵 $A$ 的 QR 分解 $A=QR$

(res <- qr(A))

$qr
           [,1]       [,2]
[1,] -1.5620499 -3.0728851
[2,]  0.7682213  0.9986877

$rank
[1] 2

$qraux
[1] 1.6401844 0.9986877

$pivot
[1] 1 2

attr(,"class")
[1] "qr"

QR 分解结果中的 Q

qr.Q(res)

           [,1]       [,2]
[1,] -0.6401844 -0.7682213
[2,] -0.7682213  0.6401844

QR 分解结果中的 R

qr.R(res)

         [,1]       [,2]
[1,] -1.56205 -3.0728851
[2,]  0.00000  0.9986877

恢复矩阵 $A$

qr.Q(res) %*% qr.R(res)

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

B.2.4 Cholesky 分解

矩阵 $A$ 的 Cholesky 分解 $A=L^{\mathrm{\top }}L$ ，其中 $L$ 是上三角矩阵

(res <- chol(A))

     [,1]  [,2]
[1,]    1 1.200
[2,]    0 1.249

t(res) %*% res

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

B.2.5 特征值分解

特征值分解（Eigenvalues Decomposition）也叫谱分解（Spectral Decomposition）

矩阵 $A$ 的特征值分解 $A=V\mathrm{\Lambda }{V}^{-1}$

(res <- eigen(A))

eigen() decomposition
$values
[1] 3.5620499 0.4379501

$vectors
          [,1]       [,2]
[1,] 0.4241554 -0.9055894
[2,] 0.9055894  0.4241554

返回值列表中的元素 vectors 就是 $V$

res$vectors %*% diag(res$values) %*% solve(res$vectors)

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

计算特征值，即求解如下一元 $n$ 次方程

$|A-\lambda I|=0$

rootSolve::uniroot.all(
  f = function(x) (x - 1) * (x - 3) - 1.2^2,
  lower = -10, upper = 10
)

[1] 0.4379747 3.5620253

B.2.6 SVD 分解

矩阵 $A$ 的 SVD 分解 $A=UD{V}^{\mathrm{\top }}$ ，矩阵 U 和 V 是正交的，矩阵 D 是对角的，矩阵 D 的对角元素是按降序排列的奇异值。

当矩阵是对称矩阵时，SVD 分解和特征值分解结果是一样的。

(res <- svd(A))

$d
[1] 3.5620499 0.4379501

$u
           [,1]       [,2]
[1,] -0.4241554 -0.9055894
[2,] -0.9055894  0.4241554

$v
           [,1]       [,2]
[1,] -0.4241554 -0.9055894
[2,] -0.9055894  0.4241554

# A = U D V'
res$u %*% diag(res$d) %*% t(res$v)

     [,1] [,2]
[1,]  1.0  1.2
[2,]  1.2  3.0

# D = U'AV
t(res$u) %*% A %*% res$v

             [,1]          [,2]
[1,] 3.562050e+00 -3.276838e-16
[2,] 3.602566e-17  4.379501e-01

# I = VV'
res$v %*% t(res$v)

              [,1]          [,2]
[1,]  1.000000e+00 -1.647255e-17
[2,] -1.647255e-17  1.000000e+00

# I = UU'
res$u %*% t(res$u)

              [,1]          [,2]
[1,]  1.000000e+00 -7.722454e-17
[2,] -7.722454e-17  1.000000e+00

B.3 Eigen 库

Eigen 是一个高性能的线性代数计算库，基于 C++ 编写，有 R 语言接口 RcppEigen 包。示例来自 RcppEigen 包，本文增加了特征向量，下面介绍如何借助 RcppEigen 包调用 Eigen 库做 SVD 矩阵分解。

#include <RcppEigen.h>

// [[Rcpp::depends(RcppEigen)]]

using Eigen::Map;                       // 'maps' rather than copies
using Eigen::MatrixXd;                  // variable size matrix, double precision
using Eigen::VectorXd;                  // variable size vector, double precision
using Eigen::SelfAdjointEigenSolver;    // one of the eigenvalue solvers

// [[Rcpp::export]]
VectorXd getEigenValues(Map<MatrixXd> M) {
  SelfAdjointEigenSolver<MatrixXd> es(M);
  return es.eigenvalues();
}
// [[Rcpp::export]]
MatrixXd getEigenVectors(Map<MatrixXd> M) {
  SelfAdjointEigenSolver<MatrixXd> es(M);
  return es.eigenvectors();
}

对上面的代码做几点说明：

1. // [[Rcpp::depends(RcppEigen)]] 可以看作一种标记，表示依赖 RcppEigen 包提供的 C++ 头文件，并导入到 C++ 命名空间中。// [[Rcpp::export]] 也可以看作一种标记，表示下面的函数需要导出到 R 语言环境中，这样 C++ 中定义的函数可以在 R 语言环境中使用。
2. MatrixXd 和 VectorXd 分别是 Eigen 库中定义的可变大小的双精度矩阵、向量类型。
3. SelfAdjointEigenSolver 是 Eigen 库中关于特征值分解方法中的一个求解器，特征值分解的结果有两个部分：一个是由特征值构成的向量，一个是特征向量构成的矩阵。求解器 SelfAdjointEigenSolver 名称中 SelfAdjoint 是伴随的意思，它是做矩阵 $A$ 的伴随矩阵 $A^{\star }$ 的特征值分解。
4. getEigenValues 和 getEigenVectors 是用户自定义的两个函数名称，分别计算特征值和特征向量。

伴随矩阵的特征值分解和原矩阵的特征值分解有何关系？为什么不直接求原矩阵的特征值分解呢？

1. 伴随矩阵的特征值与原矩阵是一样的。
2. 伴随矩阵的特征向量有一个符号差异。

RcppEigen 包封装了 Eigen 库，它在 RcppEigen 包的源码路径为

RcppEigen/inst/include/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

在 Eigen 库的源码路径如下：

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h 。

如何使用 RcppEigen 包加速计算？还是要看 Eigen 库的文档和源码，通过阅读源码，可以知道有哪些求解器，比如名称 SelfAdjointEigenSolver ，以及求解器包含的方法，比如 eigenvalues() 和 eigenvectors()，还有参数和返回值类型等。以特征值分解器 SelfAdjointEigenSolver 为例，编译上面的 C++ 代码，获得在 R 语言环境中可直接使用的函数 getEigenValues() 。

# 编译代码
Rcpp::sourceCpp(file = "code/rcpp_eigen.cpp")

然后，函数 getEigenValues() 计算特征值，返回一个向量。

# 计算特征值
getEigenValues(A)

[1] 0.4379501 3.5620499

返回一个矩阵，列是特征向量。

# 计算特征向量
getEigenVectors(A)

           [,1]       [,2]
[1,] -0.9055894 -0.4241554
[2,]  0.4241554 -0.9055894

根据上述分解结果计算矩阵 A 的伴随矩阵 $A^{\star }$ 。

t(getEigenVectors(A)) %*% diag(getEigenValues(A)) %*% getEigenVectors(A)

     [,1] [,2]
[1,]  1.0 -1.2
[2,] -1.2  3.0

B.4 应用

以线性模型为例讲述一些初步的计算性能提升办法。回顾一下线性回归的矩阵表示。

$$\begin{array}{rl}& \mathit{y}=X\mathit{\beta }+\mathit{ϵ}\\ & \mathit{ϵ}\sim \mathrm{MVN}(\mathbf{0},{\sigma }^{2}I)\end{array}$$

模型中 $\mathit{\beta },{\sigma }^2$ 是待估的参数，它们的最小二乘估计分别记为 $\hat{\mathit{\beta }},\widehat{{\sigma }^2}$ 。

$$\begin{array}{rl}\hat{\mathit{\beta }}& =(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}\\ \widehat{{\sigma }^2}& =\frac{{\mathit{y}}^{\mathrm{\top }}(I-X(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }})\mathit{y}}{n-\mathrm{rank}(X)}\end{array}$$

在获得参数的估计后，响应变量 $\mathit{y}$ 的预测 $\hat{\mathit{y}}$ 及其预测方差 $\mathsf{Var}(\hat{\mathit{y}})$ 如下。

$$\begin{array}{rl}\hat{\mathit{y}}& =X(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}\\ \mathsf{Var}(\hat{\mathit{y}})& ={\sigma }^{2}X(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}\end{array}$$

set.seed(2023)
n <- 200
p <- 50
x <- matrix(rnorm(n * p), n)
y <- rnorm(n)
fit_lm <- lm(y ~ x + 0)

下面不同的方法来计算预测值 $\hat{\mathit{y}}$ ，从慢到快地优化。教科书版就是从左至右依次计算。

fit_base = function(x, y) {
  x %*% solve(t(x) %*% x) %*% t(x) %*% y
}

矩阵乘向量比矩阵乘矩阵快。虽然矩阵乘法没有交换律，但是有结合律。先向量计算，然后矩阵计算。

$$\hat{\mathit{y}}=X(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}$$

fit_vector <- function(x, y) {
  x %*% (solve(t(x) %*% x) %*% (t(x) %*% y))
}

解线性方程组比求逆快。 $X^{\mathrm{\top }}X$ 是对称的，通过解线性方程组来避免求逆。

$$\hat{\mathit{y}}=X(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}$$

fit_inv <- function(x, y) {
  x %*% solve(crossprod(x), crossprod(x, y))
}

QR 分解。 $X_{n\times p}=Q_{n\times p}{R}_{p\times p}$，$n>p$，$Q^{\mathrm{\top }}Q=I$，$R$ 是上三角矩阵。

$$\begin{array}{rl}\hat{\mathit{y}}& =X(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}\\ & =QR((QR{)}^{\mathrm{\top }}QR{)}^{-1}(QR{)}^{\mathrm{\top }}\mathit{y}\\ & =QR(R^{\mathrm{\top }}R{)}^{-1}{R}^{\mathrm{\top }}{Q}^{\mathrm{\top }}\mathit{y}\\ & =Q{Q}^{\mathrm{\top }}\mathit{y}\end{array}$$

fit_qr <- function(x, y) {
  decomp <- qr(x)
  qr.qy(decomp, qr.qty(decomp, y))
}
fit_qr2 <- lm.fit(x, y)

其中，函数 qr.qy(decomp, y) 表示 Q %*% y ，函数 qr.qty(decomp, y) 表示 t(Q) %*% y 。实际上，Base R 提供的线性回归拟合函数 lm() 就采用 QR 分解。

Cholesky 分解。记 $A=X^{\mathrm{\top }}X$ ，若 $A$ 是正定矩阵，则 $A$ 可做 Cholesky 分解。不妨设$A=L^{\mathrm{\top }}L$，其中 $L$ 是上三角矩阵。

$$\begin{array}{rl}\hat{\mathit{y}}& =X(X^{\mathrm{\top }}X{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}\\ & =X(L^{\mathrm{\top }}L{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}\\ & =X{L}^{-1}(L^{\mathrm{\top }}{)}^{-1}{X}^{\mathrm{\top }}\mathit{y}\end{array}$$

fit_chol <- function(x, y) {
  decomp <- chol(crossprod(x))
  lxy <- backsolve(decomp, crossprod(x, y), transpose = TRUE)
  b <- backsolve(decomp, lxy)
  x %*% b
}

函数 backsolve() 求解上三角线性方程组。

Bates, Douglas M., 和 Donald G. Watts. 1988. Nonlinear Regression Analysis and Its Applications. New York, NY: John Wiley & Sons. https://doi.org/10.1002/9780470316757.app2.

Hoaglin, David C., 和 Roy E. Welsch. 1978. 《The Hat Matrix in Regression and ANOVA》. The American Statistician 32 (1): 17–22. https://www.jstor.org/stable/2683469.

---

1. https://stat.ethz.ch/pipermail/r-help/2006-March/101596.html