# R 语言中的矩阵计算

## 前言

R 是作为统计语言，生来就对数学有良好的支持。矩阵计算作为底层的数学工具，有非常广泛的使用场景。用R语言很好地封装了，矩阵的各种计算方法，一个函数一行代码，就能完成复杂的矩阵分解等操作。让建模人员可以更专注于模型推理和业务逻辑实现，把复杂的矩阵计算交给R语言来完成。 本文总结了 R 语言用于矩阵的各种计算操作。

## 1. 基本操作

```# 生成矩阵
> m<-matrix(1:20,4,5);m
[,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

# 矩阵行
> nrow(m)
[1] 4
# 矩阵列
> ncol(m)
[1] 5
```

```# 对角线元素
> diag(m)
[1]  1  6 11 16

# 以对角线元素，生成对角矩阵
> diag(diag(m))
[,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    6    0    0
[3,]    0    0   11    0
```

```# 上三角
> m<-matrix(1:20,4,5)
> upper.tri(m)
[,1]  [,2]  [,3]  [,4] [,5]
[1,] FALSE  TRUE  TRUE  TRUE TRUE
[2,] FALSE FALSE  TRUE  TRUE TRUE
[3,] FALSE FALSE FALSE  TRUE TRUE
[4,] FALSE FALSE FALSE FALSE TRUE

# 上三角值
> m[which(upper.tri(m))]
[1]  5  9 10 13 14 15 17 18 19 20

# 上三角矩阵
> m[!upper.tri(m)]<-0;m
[,1] [,2] [,3] [,4] [,5]
[1,]    0    5    9   13   17
[2,]    0    0   10   14   18
[3,]    0    0    0   15   19
[4,]    0    0    0    0   20

# 下三角矩阵
> m[!lower.tri(m)]<-0;m
[,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    2    0    0    0    0
[3,]    3    7    0    0    0
[4,]    4    8   12    0    0
```

```> m<-matrix(1:20,4,5);m
[,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20

# 转置，行列互转
> t(m)
[,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
[4,]   13   14   15   16
[5,]   17   18   19   20
```

```# 创建方阵
> m<-matrix(1:16,4,4);m
[,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 用上三角填充下三角
> m[lower.tri(m)]<-m[upper.tri(m)]
> m
[,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    5    6   10   14
[3,]    9   13   11   15
[4,]   10   14   15   16
```

```> m<-matrix(1:16,4,4);m
[,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 上三角取值
> m[upper.tri(m)]
[1]  5  9 10 13 14 15

# 下三角取值
> m[lower.tri(m)]
[1]  2  3  4  7  8 12
```

```> m<-matrix(1:20,4,5)

# 转置后，取下三角，填充下三角
> m[lower.tri(m)]<-t(m)[lower.tri(m)]
> m
[,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    5    6   10   14
[3,]    9   10   11   15
[4,]   13   14   15   16
```

```> m<-matrix(1:12,4,3);m
[,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

# 矩阵转行列data.frame
> row<-rep(1:nrow(m),ncol(m))              # 行
> col<-rep(1:ncol(m),each=nrow(m))         # 列
> df<-data.frame(row,col,v=as.numeric(m))  # 行列索引数据框
> df
row col  v
1    1   1  1
2    2   1  2
3    3   1  3
4    4   1  4
5    1   2  5
6    2   2  6
7    3   2  7
8    4   2  8
9    1   3  9
10   2   3 10
11   3   3 11
12   4   3 12

# 行列索引数据框转矩阵
> m<-matrix(0,length(unique(df\$row)),length(unique(df\$col)) )
> apply(df,1,function(dat){
+     m[dat[1],dat[2]]<<-dat[3]
+     invisible()
+ })
# 打印矩阵
> m
[,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12
```

## 2. 矩阵计算

```# 加载矩阵计算工具包
> library(matrixcalc)

# 新建2个矩阵，行列长度相同
> m0<-matrix(1:20,4,5);m0
[,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
> m1<-matrix(sample(100,20),4,5);m1
[,1] [,2] [,3] [,4] [,5]
[1,]   40   79   97   57   78
[2,]   93   32   48    8   95
[3,]   63    6   56   12    9
[4,]   28   31   72   27   26

# 矩阵加法
> m0+m1
[,1] [,2] [,3] [,4] [,5]
[1,]   41   84  106   70   95
[2,]   95   38   58   22  113
[3,]   66   13   67   27   28
[4,]   32   39   84   43   46

# 矩阵减法
> m0-m1
[,1] [,2] [,3] [,4] [,5]
[1,]  -39  -74  -88  -44  -61
[2,]  -91  -26  -38    6  -77
[3,]  -60    1  -45    3   10
[4,]  -24  -23  -60  -11   -6
```

```> m0*m1
[,1] [,2] [,3] [,4] [,5]
[1,]   40  395  873  741 1326
[2,]  186  192  480  112 1710
[3,]  189   42  616  180  171
[4,]  112  248  864  432  520
```

```> t(m0)%*%m1
[,1] [,2] [,3] [,4] [,5]
[1,]  527  285  649  217  399
[2,] 1423  877 1741  633 1231
[3,] 2319 1469 2833 1049 2063
[4,] 3215 2061 3925 1465 2895
[5,] 4111 2653 5017 1881 3727

# 通过函数实现矩阵相乘
> crossprod(m0,m1)
[,1] [,2] [,3] [,4] [,5]
[1,]  527  285  649  217  399
[2,] 1423  877 1741  633 1231
[3,] 2319 1469 2833 1049 2063
[4,] 3215 2061 3925 1465 2895
[5,] 4111 2653 5017 1881 3727
```

```> m<-matrix(1:6,2);m
[,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> t(m)
[,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6

> m %o% t(m)  # 外积，同outer(m,t(m))
, , 1, 1

[,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2, 1

[,1] [,2] [,3]
[1,]    3    9   15
[2,]    6   12   18

, , 3, 1

[,1] [,2] [,3]
[1,]    5   15   25
[2,]   10   20   30

, , 1, 2

[,1] [,2] [,3]
[1,]    2    6   10
[2,]    4    8   12

, , 2, 2

[,1] [,2] [,3]
[1,]    4   12   20
[2,]    8   16   24

, , 3, 2

[,1] [,2] [,3]
[1,]    6   18   30
[2,]   12   24   36
```

```> m0<-matrix(1:4,2,2);m0
[,1] [,2]
[1,]    1    3
[2,]    2    4
> m1<-matrix(1:9,3,3);m1
[,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> m0 %s% m1       #矩阵直和，同 direct.sum(m0,m1)
[,1] [,2] [,3] [,4] [,5]
[1,]    1    3    0    0    0
[2,]    2    4    0    0    0
[3,]    0    0    1    4    7
[4,]    0    0    2    5    8
[5,]    0    0    3    6    9
```

```> m0<-matrix(1:4,2,2);m0
[,1] [,2]
[1,]    1    3
[2,]    2    4
> m1<-matrix(1:9,3,3);m1
[,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> m0 %x% m1     #矩阵直积，同 kronecker(m0,m1)，direct.prod(m0,m1)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    4    7    3   12   21
[2,]    2    5    8    6   15   24
[3,]    3    6    9    9   18   27
[4,]    2    8   14    4   16   28
[5,]    4   10   16    8   20   32
[6,]    6   12   18   12   24   36
```

## 3. 矩阵性质

### 3.1 奇异矩阵

```# 创建一个非奇异矩阵
> m1<-matrix(sample(1:16),4)

# 行列式不等于0
> det(m1)     # !=0
[1] 14193

# 有逆矩阵
> solve(m1)   # 可逆
[,1]        [,2]        [,3]         [,4]
[1,] -0.026210104  0.09666737  0.02473050 -0.119988727
[2,] -0.002325090  0.08384415 -0.07038681  0.008173043
[3,] -0.007186641 -0.13478475  0.05516804  0.176777285
[4,]  0.074614246  0.03663778 -0.01395054 -0.080462200

# 满秩
> library(matrixcalc)
> matrix.rank(m1)    # 长度=n,满秩
[1] 4
```

```# 奇异矩阵
> m0<-matrix(1:16,4)

# 计算行列式
> det(m0)
[1] 0

# 不可逆
> solve(m0)
Error in solve.default(m0) :
Lapack routine dgesv: system is exactly singular: U[3,3] = 0

# 不是满秩
> matrix.rank(m0)     # 长度<4
[1] 2
```

### 3.2 逆矩阵

```# 创建方阵，非奇异矩阵
> m0<-matrix(sample(100,16),4,4);m0
[,1] [,2] [,3] [,4]
[1,]   24   31   80   37
[2,]   84   13   42   71
[3,]   95   62   93   86
[4,]   69   16   94   35

# 计算矩阵的逆矩阵
> solve(m0)
[,1]          [,2]        [,3]         [,4]
[1,] -0.03083680 -0.0076561475  0.01258023  0.017218514
[2,] -0.01710957 -0.0270246488  0.03152548 -0.004553923
[3,]  0.01384721 -0.0003070371 -0.00886117  0.007757524
[4,]  0.03142440  0.0282722871 -0.01541411 -0.024126340

# 逆矩阵的性质，逆矩阵与原矩阵进行矩阵乘法，得到对角矩阵。
> round(solve(m0) %*% m0)  # 对角矩阵
[,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1```

```# 创建奇异矩阵
> m<-matrix(1:16,4,4);m
[,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 逆矩阵计算失败
> solve(m)
Error in solve.default(m) :
Lapack routine dgesv: system is exactly singular: U[3,3] = 0

# 广义逆矩阵
> library(MASS) # 加载MASS包
> ginv(m)
[,1]    [,2]  [,3]    [,4]
[1,] -0.285 -0.1075  0.07  0.2475
[2,] -0.145 -0.0525  0.04  0.1325
[3,] -0.005  0.0025  0.01  0.0175
[4,]  0.135  0.0575 -0.02 -0.0975
```

`ginv` 函数计算非奇异矩阵，和 `solve()` 函数比较。

```# 非奇异矩阵
> m0<-matrix(sample(100,16),4,4);m0
[,1] [,2] [,3] [,4]
[1,]    5   61   75   74
[2,]   10   67   11   21
[3,]   29   32   92   17
[4,]   72   23   87   36

# 计算广义矩阵的逆矩阵，与solve()结果相同
> ginv(m0)
[,1]         [,2]         [,3]         [,4]
[1,] -0.0080590817  0.006534814 -0.011333076  0.018105645
[2,] -0.0041284695  0.017615769  0.005333304 -0.004308072
[3,]  0.0009226026 -0.006671759  0.016312821 -0.005707878
[4,]  0.0165261737 -0.008200731 -0.020163889  0.008112906

# 计算矩阵的逆矩阵
> solve(m0)
[,1]         [,2]         [,3]         [,4]
[1,] -0.0080590817  0.006534814 -0.011333076  0.018105645
[2,] -0.0041284695  0.017615769  0.005333304 -0.004308072
[3,]  0.0009226026 -0.006671759  0.016312821 -0.005707878
[4,]  0.0165261737 -0.008200731 -0.020163889  0.008112906
```

```> m<-matrix(1:16,4,4);m
[,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 原矩阵%*%广义逆矩阵%*%原矩阵=原矩阵
> m %*% ginv(m) %*% m
[,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 广义逆矩阵%*%原矩阵%*%广义逆矩阵=广义逆矩阵
> ginv(m) %*% m %*% ginv(m)
[,1]    [,2]  [,3]    [,4]
[1,] -0.285 -0.1075  0.07  0.2475
[2,] -0.145 -0.0525  0.04  0.1325
[3,] -0.005  0.0025  0.01  0.0175
[4,]  0.135  0.0575 -0.02 -0.0975

# 原矩阵%*%广义逆矩阵=转置后的，原矩阵%*%广义逆矩阵
> m %*% ginv(m)
[,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7
> t(m %*% ginv(m))
[,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7

# 广义逆矩阵%*%原矩阵=转置后的，广义逆矩阵%*%原矩阵
> ginv(m) %*% m
[,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7
> t(ginv(m) %*% m)
[,1] [,2] [,3] [,4]
[1,]  0.7  0.4  0.1 -0.2
[2,]  0.4  0.3  0.2  0.1
[3,]  0.1  0.2  0.3  0.4
[4,] -0.2  0.1  0.4  0.7
```

### 3.3 特征值和特征向量

```# 创建一个方阵
> m0<-matrix(sample(100,16),4,4);m0
[,1] [,2] [,3] [,4]
[1,]   97   93   11   70
[2,]   19   58   23   90
[3,]   17   94    6   35
[4,]   79   71   43   88

# 计算特征值和特征向量
> eigen(m0)
eigen() decomposition
\$values
[1] 236.14449+ 0.00000i  40.51501+ 0.00000i -13.82975+32.81166i -13.82975-32.81166i

\$vectors
[,1]          [,2]                  [,3]                  [,4]
[1,] 0.6055193+0i -0.6428709+0i  0.2114869-0.0613776i  0.2114869+0.0613776i
[2,] 0.4115920+0i  0.3939259+0i -0.0781518+0.3993580i -0.0781518-0.3993580i
[3,] 0.3054253+0i  0.6482306+0i  0.7557355+0.0000000i  0.7557355+0.0000000i
[4,] 0.6088134+0i -0.1064728+0i -0.3210016-0.3342655i -0.3210016+0.3342655i
```

`symmetric=TRUE` 时，计算对称矩阵的特征值和特征向量，当 m0 不是对称矩阵时，则取下三角对称结构进行计算。

```> eigen(m0,symmetric = TRUE)
eigen() decomposition
\$values
[1] 231.824838  87.176816  -3.422899 -66.578756

\$vectors
[,1]       [,2]       [,3]        [,4]
[1,] -0.4801643  0.7155875  0.5040624 -0.05742733
[2,] -0.5024315 -0.5470752  0.2262475 -0.63014551
[3,] -0.3634220 -0.4138943  0.3287896  0.76714627
[4,] -0.6204267  0.1316618 -0.7659181  0.10538188

# 生成下三角矩阵的对称矩阵
> m0[upper.tri(m0)]<-t(m0)[upper.tri(m0)];m0
[,1] [,2] [,3] [,4]
[1,]   97   19   17   79
[2,]   19   58   94   71
[3,]   17   94    6   43
[4,]   79   71   43   88

# 计算特征值，与eigen(m0,symmetric = TRUE) 一致
> eigen(m0)
eigen() decomposition
\$values
[1] 231.824838  87.176816  -3.422899 -66.578756

\$vectors
[,1]       [,2]       [,3]        [,4]
[1,] -0.4801643  0.7155875  0.5040624 -0.05742733
[2,] -0.5024315 -0.5470752  0.2262475 -0.63014551
[3,] -0.3634220 -0.4138943  0.3287896  0.76714627
[4,] -0.6204267  0.1316618 -0.7659181  0.10538188
```

## 4. 矩阵分解

### 4.1 三角分解 LU

```# 创建一个矩阵
> m0<-matrix(sample(100,16),4);m0
[,1] [,2] [,3] [,4]
[1,]   25   88   35   87
[2,]   26   75   73   15
[3,]   36   62   80   53
[4,]   70   44   21   47

# 三角分解
> lu<-lu.decomposition(m0)
> lu
\$L                                        # 下三角矩阵
[,1]      [,2]     [,3] [,4]
[1,] 1.00  0.000000 0.000000    0
[2,] 1.04  1.000000 0.000000    0
[3,] 1.44  3.917676 1.000000    0
[4,] 2.80 12.251816 4.617547    1

\$U                                        # 上三角矩阵
[,1]   [,2]      [,3]      [,4]
[1,]   25  88.00   35.0000   87.0000
[2,]    0 -16.52   36.6000  -75.4800
[3,]    0   0.00 -113.7869  223.4262
[4,]    0   0.00    0.0000 -303.5137

# 三角分解的下三角矩阵L %*% 三角分解的上三角矩阵U = 原矩阵
> lu\$L %*% lu\$U
[,1] [,2] [,3] [,4]
[1,]   25   88   35   87
[2,]   26   75   73   15
[3,]   36   62   80   53
[4,]   70   44   21   47
```

### 4.2 choleskey 分解

Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的。Cholesky 分解法又称平方根法，是当A为实对称正定矩阵时，LU 三角分解法的变形。

```# 创建对称方阵
> m1<-diag(5)+2;m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# choleskey分解
> chol(m1)
[,1]     [,2]      [,3]      [,4]      [,5]
[1,] 1.732051 1.154701 1.1547005 1.1547005 1.1547005
[2,] 0.000000 1.290994 0.5163978 0.5163978 0.5163978
[3,] 0.000000 0.000000 1.1832160 0.3380617 0.3380617
[4,] 0.000000 0.000000 0.0000000 1.1338934 0.2519763
[5,] 0.000000 0.000000 0.0000000 0.0000000 1.1055416

# choleskey分解的性质
# chol分解的逆矩阵 %*% chol分解矩阵 = 原矩阵
> t(chol(m1))%*%chol(m1)
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# chol分解矩阵的对角线值的平方的积 = 行列式的模数
> prod(diag(chol(m1))^2)
[1] 11
> det(m1)
[1] 11

# chol分解的逆
> chol2inv(m1)
[,1]         [,2]        [,3]         [,4]         [,5]
[1,]  0.166658199 -0.055580958 -0.01859473 -0.006401463 -0.002743484
[2,] -0.055580958  0.166590459 -0.05578418 -0.019204390 -0.008230453
[3,] -0.018594726 -0.055784179  0.16598080 -0.057613169 -0.024691358
[4,] -0.006401463 -0.019204390 -0.05761317  0.160493827 -0.074074074
[5,] -0.002743484 -0.008230453 -0.02469136 -0.074074074  0.111111111

> chol2inv(chol(m1))
[,1]       [,2]       [,3]       [,4]       [,5]
[1,]  0.8181818 -0.1818182 -0.1818182 -0.1818182 -0.1818182
[2,] -0.1818182  0.8181818 -0.1818182 -0.1818182 -0.1818182
[3,] -0.1818182 -0.1818182  0.8181818 -0.1818182 -0.1818182
[4,] -0.1818182 -0.1818182 -0.1818182  0.8181818 -0.1818182
[5,] -0.1818182 -0.1818182 -0.1818182 -0.1818182  0.8181818
```

### 4.3 QR 分解

QR 分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵，所以称为 QR 分解法，与此正规正交矩阵的通用符号 Q 有关。

```# 创建对称方阵
> m1<-diag(5)+2;m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# QR分解
> q<-qr(m1);q
\$qr
[,1]       [,2]       [,3]       [,4]       [,5]
[1,] -5.0 -4.8000000 -4.8000000 -4.8000000 -4.8000000
[2,]  0.4 -1.4000000 -0.6857143 -0.6857143 -0.6857143
[3,]  0.4  0.2142857 -1.2205720 -0.4012839 -0.4012839
[4,]  0.4  0.2142857  0.1560549 -1.1527216 -0.2852095
[5,]  0.4  0.2142857  0.1560549  0.1246843  1.1168808

\$rank
[1] 5

\$qraux
[1] 1.600000 1.928571 1.975343 1.992196 1.116881

\$pivot
[1] 1 2 3 4 5

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

# 计算QR分解矩阵的R值
> qr.R(q)
[,1] [,2]       [,3]       [,4]       [,5]
[1,]   -5 -4.8 -4.8000000 -4.8000000 -4.8000000
[2,]    0 -1.4 -0.6857143 -0.6857143 -0.6857143
[3,]    0  0.0 -1.2205720 -0.4012839 -0.4012839
[4,]    0  0.0  0.0000000 -1.1527216 -0.2852095
[5,]    0  0.0  0.0000000  0.0000000  1.1168808

# 计算QR分解矩阵的Q值
> qr.Q(q)
[,1]        [,2]        [,3]        [,4]       [,5]
[1,] -0.6  0.62857143  0.36784361  0.26144202 -0.2030692
[2,] -0.4 -0.77142857  0.36784361  0.26144202 -0.2030692
[3,] -0.4 -0.05714286 -0.85272836  0.26144202 -0.2030692
[4,] -0.4 -0.05714286 -0.03344033 -0.89127960 -0.2030692
[5,] -0.4 -0.05714286 -0.03344033 -0.02376746  0.9138115

# QR分解的性质
# Q的值 %*% R值 = 原矩阵
> qr.Q(q) %*% qr.R(q) #=m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# X值 = 原矩阵，同Q的值 %*% R值
> qr.X(q) #=m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# Q值的转置矩阵 * Q值 =
> t(qr.Q(q)) %*% qr.Q(q)
[,1]          [,2]          [,3]         [,4]          [,5]
[1,]  1.000000e+00 -3.469447e-17 -2.081668e-17 1.214306e-17  5.551115e-17
[2,] -3.469447e-17  1.000000e+00  7.199102e-17 5.204170e-17 -7.632783e-17
[3,] -2.081668e-17  7.199102e-17  1.000000e+00 3.382711e-17 -6.938894e-18
[4,]  1.214306e-17  5.204170e-17  3.382711e-17 1.000000e+00  3.122502e-17
[5,]  5.551115e-17 -7.632783e-17 -6.938894e-18 3.122502e-17  1.000000e+00
```

### 4.4 奇异值分解 SVD

```# 创建对称方阵
> m1<-diag(5)+2;m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# 奇异值分解
> s<-svd(m1);s
\$d
[1] 11  1  1  1  1

\$u
[,1]          [,2]          [,3]          [,4]       [,5]
[1,] -0.4472136  4.440892e-17 -3.330669e-17 -3.108624e-16  0.8944272
[2,] -0.4472136 -3.219647e-16  2.414735e-16  8.660254e-01 -0.2236068
[3,] -0.4472136 -5.773503e-01  5.773503e-01 -2.886751e-01 -0.2236068
[4,] -0.4472136 -2.113249e-01 -7.886751e-01 -2.886751e-01 -0.2236068
[5,] -0.4472136  7.886751e-01  2.113249e-01 -2.886751e-01 -0.2236068

\$v
[,1]          [,2]          [,3]       [,4]       [,5]
[1,] -0.4472136  0.000000e+00  0.000000e+00  0.0000000  0.8944272
[2,] -0.4472136 -1.665335e-16  1.457168e-16  0.8660254 -0.2236068
[3,] -0.4472136 -5.773503e-01  5.773503e-01 -0.2886751 -0.2236068
[4,] -0.4472136 -2.113249e-01 -7.886751e-01 -0.2886751 -0.2236068
[5,] -0.4472136  7.886751e-01  2.113249e-01 -0.2886751 -0.2236068

# 奇异分解性质
# svd的u矩阵 %*% svd的d矩阵的对角性值 %*% svd的v的转置矩阵 = 原矩阵
> s\$u %*% diag(s\$d) %*% t(s\$v) #=m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# svd的u矩阵的转置矩阵 %*% svd的u矩阵 = svd的s矩阵的转置矩阵 %*% svd的v矩阵
> t(s\$u) %*% s\$u
[,1]          [,2]          [,3]          [,4]          [,5]
[1,]  1.000000e+00  0.000000e+00 -5.551115e-17  0.000000e+00 -2.775558e-17
[2,]  0.000000e+00  1.000000e+00 -2.775558e-17 -1.387779e-16  8.326673e-17
[3,] -5.551115e-17 -2.775558e-17  1.000000e+00  6.245005e-17 -6.938894e-17
[4,]  0.000000e+00 -1.387779e-16  6.245005e-17  1.000000e+00 -1.387779e-16
[5,] -2.775558e-17  8.326673e-17 -6.938894e-17 -1.387779e-16  1.000000e+00
> t(s\$v) %*% s\$v
[,1]          [,2]         [,3]          [,4]          [,5]
[1,]  1.000000e+00  0.000000e+00 5.551115e-17 -1.665335e-16  2.775558e-17
[2,]  0.000000e+00  1.000000e+00 8.326673e-17 -8.326673e-17  0.000000e+00
[3,]  5.551115e-17  8.326673e-17 1.000000e+00  1.595946e-16  2.775558e-17
[4,] -1.665335e-16 -8.326673e-17 1.595946e-16  1.000000e+00 -8.326673e-17
[5,]  2.775558e-17  0.000000e+00 2.775558e-17 -8.326673e-17  1.000000e+00```

## 5. 特殊矩阵

### 5.1 Hankel Matrix

```> hankel.matrix(6, 1:50)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    2    3    4    5    6    7
[3,]    3    4    5    6    7    8
[4,]    4    5    6    7    8    9
[5,]    5    6    7    8    9   10
[6,]    6    7    8    9   10   11
```

### 5.2 Hilbert Matrix

```> hilbert.matrix(4)
[,1]      [,2]      [,3]      [,4]
[1,] 1.0000000 0.5000000 0.3333333 0.2500000
[2,] 0.5000000 0.3333333 0.2500000 0.2000000
[3,] 0.3333333 0.2500000 0.2000000 0.1666667
[4,] 0.2500000 0.2000000 0.1666667 0.1428571
```

### 5.3 Creation Matrix

```> creation.matrix(5)
[,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    1    0    0    0    0
[3,]    0    2    0    0    0
[4,]    0    0    3    0    0
[5,]    0    0    0    4    0
```

### 5.4 Stirling Matrix

```> stirling.matrix(4)
[,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    1    1    0    0
[4,]    0    1    3    1    0
[5,]    0    1    7    6    1```

### 5.5 Pascal matrix

```> pascal.matrix(4)
[,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    1    0    0
[3,]    1    2    1    0
[4,]    1    3    3    1
```

### 5.6 Fibonacci matrix

```> fibonacci.matrix(4)
[,1] [,2] [,3] [,4] [,5]
[1,]    1    1    0    0    0
[2,]    2    1    1    0    0
[3,]    3    2    1    1    0
[4,]    5    3    2    1    1
[5,]    8    5    3    2    1
```

### 5.7 Frobenius Matrix

Frobenius 矩阵也称为伴随矩阵，它出现在线性一阶微分方程组的解中。此函数返回一个在数值数学中有用的 Fronenius 矩阵。 矩阵定义：

```> frobenius.matrix(4)
[,1] [,2] [,3] [,4]
[1,]    0    0    0   -1
[2,]    1    0    0    4
[3,]    0    1    0   -6
[4,]    0    0    1    4
```

### 5.8 Duplication matrix

```> D.matrix(3)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    0    0    0    0
[2,]    0    1    0    0    0    0
[3,]    0    0    1    0    0    0
[4,]    0    1    0    0    0    0
[5,]    0    0    0    1    0    0
[6,]    0    0    0    0    1    0
[7,]    0    0    1    0    0    0
[8,]    0    0    0    0    1    0
[9,]    0    0    0    0    0    1
```

### 5.9 K matrix

k 矩阵是由 `H.matrices()` 函数构造的，利用直积进行计算子列表的分量。`K.matrix(r, c=r)` ，返回阶数为 `p=r*c` 的方阵，对于 r 行 c 列的矩阵 A，计算 `A``t(A)` 的直积。 计算公式：

```> K.matrix(2,3)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    0    0    0    0
[2,]    0    0    1    0    0    0
[3,]    0    0    0    0    1    0
[4,]    0    1    0    0    0    0
[5,]    0    0    0    1    0    0
[6,]    0    0    0    0    0    1
```

### 5.10 E Matrices

E 矩阵列表， `E.matrices(n)` 使得每个子列表的分量，是从 n 阶单位矩阵计算向量的外积导出的方阵。 计算公式：

```> E.matrices(3)
[[1]]
[[1]][[1]]
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0
[3,]    0    0    0

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

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

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

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

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

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

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

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

### 5.11 H Matrices

H 矩阵列表， `H.matrices(r, c=r)` 使得 r 阶 c 阶的子列表的分量，计算从 r 行和 c 列的单位矩阵的列向量的外积导出的方阵。

```> H.matrices(2,3)
[[1]]
[[1]][[1]]
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0

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

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

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

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

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

### 5.12 T Matrices

T 矩阵列表， `T.matrices(n)` 高级别列表中的组件数为 n。n 个组件中的每一个也是列表。每个子列表具有 n 个分量，每个分量是 n 阶矩阵。 计算公式：

```> T.matrices(3)
[[1]]
[[1]][[1]]
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0
[3,]    0    0    0

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

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

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

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

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

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

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

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

—END—

0 条评论

• ### 「Workshop」第十七期 奇异值分解

奇异值分解(Singular Value Decomposition，以下简称SVD)是在机器学习领域广泛应用的算法，它不光可以用于降维算法中的特征分解，还可以...

• ### 「R」使用NMF包绘制热图

学习文档: https://cran.r-project.org/web/packages/NMF/vignettes/heatmaps.pdf

• ### 绘图工具？你应该有‘掰弯’的能力

写推文，我从来就不是纯粹地炫技，因为一旦有了TBtools，所有人都能轻易掌握并做得出来我推出来的所有图。

• ### 《Unity3D 实战核心技术详解》书中关于矩阵的错误

最近一直在学习实时渲染，不免要接触线性代数。而渲染中，一定会用到矩阵，当我再次去复习我之前看的书时，发现《Unity3D 实战核心技术详解》关于矩阵就有几处错误...

• ### 彻底理解矩阵乘法

今天的角度比较清奇，我们来讲讲矩阵的乘法。当然了，我告诉你的肯定不是大学教科书上那些填鸭式的云里雾里的计算规则，你可能将规则背下来了，但完全不理解为什么会这样。...

• ### 矩阵相乘法则和技巧

只有当第一个矩阵（左矩阵）的列数等于第二矩阵（右矩阵）时，两矩阵才能相乘。因为得到的结果矩阵的i一行的第j个元素（Cij）是左矩阵第i行所有元素分别与右矩阵第...

• ### 数据分析与数据挖掘 - 06线性代数

导数是高等数学中非常重要的知识点，也是人工智能的算法应用中比较常用的一个知识，这一章我们的重点就是讲解一下导数和其求导法则。首先我们来看一下导数的基本概念:函数...

• ### 数学实验（预习）

也可以用初等变换求逆矩阵,构造一个n行2n列的矩阵(A E),并进行初等变换,A编程单位矩阵的时候,E就变成了A的逆矩阵.

• ### 吹弹牛皮之Unity 引擎基础 - 矩阵（三）

上图中展示了p，q两个基向量（单位向量）绕原点旋转后得到的新基向量p'和q'。根据勾股定理有：