# 附录 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 1

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 <- 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 对数、指数与幂

expm::logm(A)
           [,1]      [,2]
[1,] -0.4485660 0.8050907
[2,]  0.8050907 0.8932519

$\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

(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

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

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$

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^{\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

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$$ - 范数

$$\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
.. ..$: 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

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

$$|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 分解

(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) {
return es.eigenvalues();
}
// [[Rcpp::export]]
MatrixXd getEigenVectors(Map<MatrixXd> M) {
return es.eigenvectors();
}

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

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

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

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

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

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

# 计算特征值
getEigenValues(A)
[1] 0.4379501 3.5620499

# 计算特征向量
getEigenVectors(A)
           [,1]       [,2]
[1,] -0.9055894 -0.4241554
[2,]  0.4241554 -0.9055894

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}

\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}

\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)

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))
}

$\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)

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
}