matrix-operations.qmd

# 矩阵运算 {#sec-matrix-operations}

::: hidden
$$
 \def\bm#1{{\boldsymbol #1}}
$$
:::

> 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 [^matrix-operations-1]

[^matrix-operations-1]: <https://stat.ethz.ch/pipermail/r-help/2006-March/101596.html>

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

```{r}
library(Matrix)
```

## 基础运算 {#sec-basic-matrix-operations}

约定符号

$$
A = \begin{bmatrix}
 a_{11} & a_{12}  & a_{13} \\ 
 a_{21} & a_{22}  & a_{23} \\ 
 a_{31} & a_{32}  & a_{33} 
\end{bmatrix}
$$

### 加、减、乘

矩阵 $A$

```{r}
A <- matrix(c(1, 1.2, 1.2, 3), nrow = 2)
A
B <- matrix(c(1, 2, 3, 4), nrow =2)
B
```

```{r}
A + A # 对应元素相加
A - A # 对应元素相减
A %*% A # 矩阵乘法
```

### 对数、指数与幂 {#sec-log-exp}

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

```{r}
expm::logm(A)
```

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

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

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

```{r}
expm::expm(A)
```

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

```{r}
(res <- svd(A))
res$u %*% diag(exp(res$d)) %*% res$v
```

矩阵 $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$

```{r}
res$u %*% (diag(res$d)^3) %*% res$v
```

### 迹、秩、条件数

矩阵 $A$ 的迹 $\operatorname{tr}(A) = \sum_{i=1}^{n}a_{ii}$

```{r}
sum(diag(A))
qr(A)$rank
kappa(A)
```

### 求逆与广义逆

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

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

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

```{r}
solve(A) # 逆
MASS::ginv(A) # 广义逆
```

### 行列式与伴随 {#sec-det-adjust}

矩阵必须是方阵

伴随矩阵 $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 次伴随是？

```{r}
det(A)
det(A) * solve(A)
```

### 外积、直积与交叉积 {#sec-crossproduct}

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

```{r}
A %*% B
```

外积

```{r}
A %o% B # outer(A, B, FUN = "*")
```

直积/克罗内克积

```{r}
A %x% B # kronecker(A, B, FUN = "*")
```

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

```{r}
crossprod(A, A)  #  t(x) %*% y
tcrossprod(A, A) #  x %*% t(y)
```

### Hadamard 积 {#subsec-hadamard-product}

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

```{r}
fastmatrix::hadamard(A, B)
```

```{r}
A^2     # 每个元素平方 a_ij ^ 2
A ** A  # 每个元素的幂 a_ij ^ a_ij
2^A     # 每个元素的指数 2 ^ a_ij
exp(A)  # 每个元素的指数 exp(a_ij)
```

### 矩阵范数 {#subsec-matrix-norm}

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

$1$-范数

:   列和绝对值最大的

$2$ - 范数

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

$\infty$ - 范数

:   行和绝对值最大的

Frobenius - 范数

:   Euclidean 范数

$M$ - 范数

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

```{r}
norm(A, type = "1") # max(abs(colSums(A)))
norm(A, type = "I") # max(abs(rowSums(A)))
norm(A, type = "F")
norm(A, type = "M") #
norm(A, type = "2") # max(svd(A)$d)
```

### 转置与旋转 {#sec-transpose-ratate}

矩阵 $A$

```{r}
t(A) # 转置
```

### 正交与投影 {#sec-orthogonal-projection}

矩阵 $A$ 的投影

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

```{r}
diag(rep(1, 2)) - A %*% solve(t(A) %*% A) %*% t(A)
```

### Givens 变换(\*) {#sec-matrix-givens}

-   [Givens 旋转](https://www.wikiwand.com/en/Givens_rotation)
-   帽子矩阵在统计中的应用，回归与方差分析 [@David1978]

### Householder 变换(\*) {#sec-matrix-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} 
$$ {#eq-householder-negative}

或者

$$
\begin{aligned}
\mathbf{u} = \frac{\mathbf{y} + \| \mathbf{y} \|\mathbf{e}_{1}}{\| \mathbf{y} + \| \mathbf{y} \|\mathbf{e}_{1}\|}
\end{aligned} 
$$ {#eq-householder-positive}

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

举个例子

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

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

用 Householder 变换做 QR 分解 [@Bates1988] 及其 [R 语言](https://rpubs.com/aaronsc32/qr-decomposition-householder)、Eigen 实现。

### 单位矩阵 {#sec-identity-matrix}

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

$$
A = \begin{bmatrix}
 1 & 0  & 0 \\ 
 0 & 1  & 0 \\ 
 0 & 0  & 1 
\end{bmatrix}
$$

```{r}
diag(rep(3))
```

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

```{r}
rep(1,3) %o% rep(1,3) 
```

### 对角矩阵 {#sec-matrix-diagonals}

```{r}
diag(A)       # 矩阵的对角
diag(x = c(1, 2, 3)) # 构造对角矩阵
```

### 稀疏矩阵 {#sec-sparse-matrix}

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

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

### 上、下三角矩阵 {#sec-upper-matrix}

```{r}
m <- A
m
upper.tri(m) # 矩阵上三角
m[upper.tri(m)]
m[lower.tri(m)] <- 0 # 获得上三角矩阵
m
```

矩阵 A 的下三角矩阵

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

## 矩阵分解 {#sec-matrix-decomposition}

### LU 分解 {#sec-lu}

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

```{r}
Matrix::lu(A)
```

### Schur 分解 {#sec-schur}

矩阵 $A$ 的 Schur 分解 $A = QTQ^{\top}$

```{r}
(res <- Matrix::Schur(A))
```

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

```{r}
res$Q %*% t(res$Q)
```

```{r}
res$Q %*% res$T %*% t(res$Q)
```

### QR 分解 {#sec-qr}

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

```{r}
(res <- qr(A))
```

QR 分解结果中的 Q

```{r}
qr.Q(res)
```

QR 分解结果中的 R

```{r}
qr.R(res)
```

恢复矩阵 $A$

```{r}
qr.Q(res) %*% qr.R(res)
```

### Cholesky 分解 {#sec-cholesky}

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

```{r}
(res <- chol(A))
```

```{r}
t(res) %*% res
```

### 特征值分解 {#sec-spectral}

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

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

```{r}
(res <- eigen(A))
```

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

```{r}
res$vectors %*% diag(res$values) %*% solve(res$vectors)
```

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

$|A - \lambda I| = 0$

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

### SVD 分解 {#sec-svd}

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

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

```{r}
(res <- svd(A))
```

```{r}
# A = U D V'
res$u %*% diag(res$d) %*% t(res$v)
# D = U'AV
t(res$u) %*% A %*% res$v
# I = VV'
res$v %*% t(res$v)
# I = UU'
res$u %*% t(res$u)
```

## Eigen 库 {#sec-eigen-library}

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

``` {#rcpp-eigen .cpp}
#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()` 。

```{r}
#| message: false
# 编译代码
Rcpp::sourceCpp(file = "code/rcpp_eigen.cpp")
```

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

```{r}
# 计算特征值
getEigenValues(A)
```

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

```{r}
# 计算特征向量
getEigenVectors(A)
```

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

```{r}
t(getEigenVectors(A)) %*% diag(getEigenValues(A)) %*% getEigenVectors(A)
```

## 应用 {#sec-matrix-linear-regression}

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

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

```{r}
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}}$ ，从慢到快地优化。教科书版就是从左至右依次计算。

```{r}
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}
$$

```{r}
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}
$$

```{r}
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}
$$

```{r}
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}
$$

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

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