附录 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 = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \]

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}^{\infty}\frac{A^k}{k!} \]

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

expm::expm(A)

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

或者使用奇异值分解 \(A = UDV^{\top}\) ，则 \(\mathrm{e}^A = U\mathrm{e}^DV^{\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 = UDV^{\top}\)

\[ \begin{aligned} A^n &= A \times A \times \cdots \times A \\ & = UDV^{\top} UDV^{\top} \cdots UDV^{\top} \end{aligned} \]

计算 \(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\) 的迹 \(\operatorname{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^{\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*A^{\star} = A^{\star} *A = |A|*I, A^{\star} = |A|*A^{-1}\)

\(|A^{\star}| = |A|^{n-1}, A \in \mathbb{R}^{n\times n},n \geq 2\)
\((A^{\star})^{\star} = |A|^{n-2}A, A \in \mathbb{R}^{n\times n},n \geq 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^{\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} \]

\[ \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \circ \begin{bmatrix} b_{11} & b_{12} & b_{13} \\ b_{21} & b_{22} & b_{23} \\ b_{31} & b_{32} & b_{33} \end{bmatrix} = \begin{bmatrix} 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{bmatrix} \]

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\) - 范数: 又称谱范数，矩阵最大的奇异值，如果是方阵，就是最大的特征值
\(\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^{\top}A)^{-1}A^{\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}^{\top} \]

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

\[ \mathbf{H}_{u}^{\top} = \mathbf{I}^{\top} - 2\mathbf{u}\mathbf{u}^{\top} = \mathbf{H}_{u} \]

并且

\[ \mathbf{H}_{u}^{\top}\mathbf{H}_{u} = \mathbf{I} -4\mathbf{u}\mathbf{u}^{\top} + 4\mathbf{u}\mathbf{u}^{\top}\mathbf{u}\mathbf{u}^{\top} = \mathbf{I} \]

让 \(\mathbf{H}_{u}\) 乘以向量 \(\mathbf{y}\)，即

\[ \mathbf{H}_{u}\mathbf{y} = \mathbf{y} - 2\mathbf{u}\mathbf{u}^{\top}\mathbf{y} \]

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

\[ \begin{aligned} \mathbf{u} = \frac{\mathbf{y} - \| \mathbf{y} \|\mathbf{e}_{1}}{\| \mathbf{y} - \| \mathbf{y} \|\mathbf{e}_{1}\|} \end{aligned} \tag{B.1}\]

或者

\[ \begin{aligned} \mathbf{u} = \frac{\mathbf{y} + \| \mathbf{y} \|\mathbf{e}_{1}}{\| \mathbf{y} + \| \mathbf{y} \|\mathbf{e}_{1}\|} \end{aligned} \tag{B.2}\]

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

举个例子

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

更快、数值更稳定比直接构造 Q，特别当 N 大于 P 的时候
相比于存储矩阵 Q 的 \(N^2\) 个元素，Householder 变换只存储 P 个向量 \(u_1,\ldots,u_P\)
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 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \]

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 = QTQ^{\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^{\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\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 = UDV^{\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();
}

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

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

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

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

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{aligned} &\boldsymbol{y} = X\boldsymbol{\beta} + \boldsymbol{\epsilon} \\ &\boldsymbol{\epsilon} \sim \mathrm{MVN}(\boldsymbol{0}, \sigma^2I) \end{aligned} \]

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

\[ \begin{aligned} \hat{\boldsymbol{\beta}} &= (X^{\top}X)^{-1}X^{\top}\boldsymbol{y} \\ \hat{\sigma^2} &= \frac{\boldsymbol{y}^{\top}(I - X(X^{\top}X)^{-1}X^{\top})\boldsymbol{y}}{n - \mathrm{rank}(X)} \end{aligned} \]

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

\[ \begin{aligned} \hat{\boldsymbol{y}} &= X(X^{\top}X)^{-1}X^{\top}\boldsymbol{y} \\ \mathsf{Var}(\hat{\boldsymbol{y}}) & = \sigma^2 X(X^{\top}X)^{-1}X^{\top} \end{aligned} \]

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

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

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

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

\[ \hat{\boldsymbol{y}} = X(X^{\top}X)^{-1}X^{\top}\boldsymbol{y} \]

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

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

\[ \hat{\boldsymbol{y}} = X(X^{\top}X)^{-1}X^{\top}\boldsymbol{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^{\top}Q = I\)，\(R\) 是上三角矩阵。

\[ \begin{aligned} \hat{\boldsymbol{y}} &= X(X^{\top}X)^{-1}X^{\top}\boldsymbol{y} \\ & = QR \big((QR)^{\top}QR\big)^{-1}(QR)^{\top}\boldsymbol{y} \\ & = QR(R^{\top}R)^{-1}R^{\top}Q^{\top}\boldsymbol{y} \\ & = QQ^{\top}\boldsymbol{y} \end{aligned} \]

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^{\top}X\) ，若 \(A\) 是正定矩阵，则 \(A\) 可做 Cholesky 分解。不妨设\(A = L^{\top}L\)，其中 \(L\) 是上三角矩阵。

\[ \begin{aligned} \hat{\boldsymbol{y}} &= X(X^{\top}X)^{-1}X^{\top}\boldsymbol{y} \\ & = X\big(L^{\top}L\big)^{-1}X^{\top}\boldsymbol{y} \\ & = XL^{-1}(L^{\top})^{-1}X^{\top}\boldsymbol{y} \end{aligned} \]

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

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

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