matrix-operations.Rmd

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

<!-- 介绍矩阵代数内容：各类矩阵分解及其几何解释，比如矩阵行列式 -->

参考 [matlib](https://github.com/friendly/matlib) 和 Matrix 包，[SparseM](https://CRAN.R-project.org/package=SparseM) 更加强调稀疏矩阵的 Cholesky 分解和后退法，矩阵取子集和 Kronecker 积。[fastmatrix](https://github.com/faosorios/fastmatrix)、[abind](https://cran.r-project.org/package=abind) 各种矩阵操作。

矩阵计算一般介绍参考在线书籍 Stephen Boyd and Lieven Vandenberghe 最新著作 Introduction to Applied Linear Algebra -- Vectors, Matrices, and Least Squares [@Boyd_2018_matrix] 及其 [Julia 语言实现](https://web.stanford.edu/~boyd/vmls/vmls-julia-companion.pdf)，

矩阵分解部分参考 [Introduction to Linear Algebra, 5th Edition](http://math.mit.edu/~gs/linearalgebra/linearalgebra5_Matrix.pdf)、[Linear Algebra for Data Science with examples in R](https://shainarace.github.io/LinearAlgebra)

[Evan Chen](https://web.evanchen.cc/) 的书[An Infinitely Large Napkin](https://github.com/vEnhance/napkin) 介绍矩阵代数的内积空间、群、环、域等高级内容，作者提供免费的电子书。

<!-- 解释基本的矩阵概念，比如秩、迹、行列式、直积等 -->

分块矩阵操作，各类分解算法，及其 R 实现

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

以 attitude 数据集为例介绍各种矩阵操作

```{r}
head(attitude)
```

rating 总体评价 complaints 处理员工投诉 privileges 不允许特权 learning 学习机会 raises 根据表现晋升 critical 批评 advancel 进步

```{r}
fit <- lm(rating ~. , data = attitude)
summary(fit) # 模型是显著的，很多变量的系数不显著
anova(fit)
```

<!-- 矩阵变换对应的几何意义，和实际数据对应的意义 -->

```{r}
attitude_mat <- as.matrix.data.frame(attitude)
# 生成演示用的矩阵
demo_mat <- t(attitude_mat[, -1]) %*% attitude_mat[, -1]
```

## SVD 分解 {#subsec-Singular-Value-Decomposition}

[Fast truncated singular value decompositions](https://github.com/bwlewis/irlba) 奇异值分解是特征值分解的推广

邱怡轩将奇异值分解用于图像压缩 <https://cosx.org/2014/02/svd-and-image-compression> 并制作了 [Shiny App](https://yihui.shinyapps.io/imgsvd/) 交互式演示

## 特殊函数 {#sec-special-functions}

### 阶乘 {#factorial}

-   阶乘 $n! = 1\times 2\times 3\cdots n$
-   双阶乘 $(2n+1)!! = 1 \times 3\times 5 \times \cdots \times (2n+1), n = 0,1,2,\cdots$

```{r}
factorial(5) # 阶乘
seq(from = 1, to = 5, length.out = 3)
prod(seq(from = 1, to = 5, length.out = 3)) # 连乘 双阶乘

seq(5)
cumprod(seq(5)) # 累积
cumsum(seq(5)) # 累和
```

此外还有 `cummax` 和 `cummin`

-   组合数 $C_{n}^{k} = \frac{n(n-1)…(n-k+1)}{k!}$

$C_{5}^{3} = \frac{5 \times 4 \times 3}{3 \times 2 \times 1}$

```{r}
choose(5,3)
```

> 斯特林公式

### 伽马函数 {#gamma}

$\Gamma(x) = \int_{0}^{\infty} t^{x-1}\exp(-t)dt$ $\Gamma(n) = (n-1)!, n \in \mathbb{Z}^{+}$

```{r gamma-ex}
gamma(2)
gamma(10)
gamma2 <- function(t,x){
  t^(x-1)*exp(-t)
}
integrate(gamma2, lower = 0, upper = + Inf, x = 10)
```

-   `psigamma(x, deriv)` 表示 $\psi(x)$ 的 `deriv` 阶导数

$\mathrm{digamma}(x) \triangleq \psi(x) = \frac{d}{dx}{\ln \Gamma(x)} = \Gamma'(x) / \Gamma(x)$

```{r psigamma-ex}
# 例1
x <- 2
eval(deriv(~ gamma(x), "x"))/gamma(x)
# 与此等价
psigamma(2, 0)

digamma(x) # psi(x) 的一阶导数
trigamma(x) # psi(x) 的二阶导数

# 例2
eval(deriv(~ psigamma(x, 1), "x"))
# 与此等价
psigamma(2, 2)

# 注意与下面这个例子比较
dx2x <- deriv(~ x^3, "x")
eval(dx2x)
```

### 贝塔函数 {#beta}

$B(a,b) = \Gamma(a)\Gamma(b)/\Gamma(a+b) = \int_{0}^{1} t^{a-1} (1-t)^{b-1} dt$

```{r beta-ex}
beta(1, 1)
beta(2, 3)
beta2 <- function(t, a, b) {
  t^(a - 1) * (1 - t)^(b - 1)
}
integrate(beta2, lower = 0, upper = 1, a = 2, b = 3)
```

### 贝塞尔函数 {#bessel}

``` r
besselI(x, nu, expon.scaled = FALSE) # 修正的第一类
besselK(x, nu, expon.scaled = FALSE) # 修正的第二类
besselJ(x, nu) # 第一类
besselY(x, nu) # 第二类
```

-   $\nu$ 贝塞尔函数的阶，可以是分数
-   `expon.scaled` 是否使用指数表示

```{r modified-bessel}
nus <- c(0:5, 10, 20)
x <- seq(0, 4, length.out = 501)
plot(x, x,
  ylim = c(0, 6), ylab = "", type = "n",
  main = "Bessel Functions  I_nu(x)"
)
for (nu in nus) lines(x, besselI(x, nu = nu), col = nu + 2)
legend(0, 6, legend = paste("nu=", nus), col = nus + 2, lwd = 1)
```

介绍复数矩阵的计算，矩阵的指数计算，介绍一点分形

```{r,eval=FALSE}
# 考虑用 gganimate 实现，去掉 caTools 依赖
library(caTools)
jet.colors <- colorRampPalette(c(
  "green", "blue", "red", "cyan", "#7FFF7F",
  "yellow", "#FF7F00", "red", "#7F0000"
))
m <- 1000 # define size
C <- complex(
  real = rep(seq(-1.8, 0.6, length.out = m), each = m),
  imag = rep(seq(-1.2, 1.2, length.out = m), m)
)
C <- matrix(C, m, m) # reshape as square matrix of complex numbers
Z <- 0 # initialize Z to zero
X <- array(0, c(m, m, 20)) # initialize output 3D array
for (k in 1:20) { # loop with 20 iterations
  Z <- Z^2 + C    # the central difference equation
  X[, , k] <- exp(-abs(Z)) # capture results
}
write.gif(X, "Mandelbrot.gif", col = jet.colors, delay = 100)
```