【R语言干货】概括总结矩阵操作计算
来源：http://blog.fens.me/r-matrix/
前言
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
填充后，发现矩阵并不是对称的，原因是上三角取值按列取值，所以先取 10 后取 13，导致上三角和下三角取值顺序不完全一致。
> 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
矩阵和 data.frame 转换，用行列形成索引结构。
> 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
矩阵乘法，满足第二个矩阵的列数和第一个矩阵的行数相等，所以把上面生成的 m0 矩阵( 4 行 5 列)转置为( 5 行 4 列)，再用 m1 矩阵( 4 行 5 列)，进行矩阵乘法，得到一个 5 行 5 列的结果矩阵。
> 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 奇异矩阵
首先，我们线创建一个非奇异矩阵，判断非奇异矩阵方法，行列式不等于 0，矩阵可逆，满秩。
# 创建一个非奇异矩阵
> 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
再创建一个奇异矩阵，判断奇异矩阵方法包括，行列式等于 0，矩阵不可逆，不是满秩。
# 奇异矩阵
> 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 种矩阵常用的分解的方法，包括三角分解 LU，choleskey 分解，QR 分解，奇异值分解 SVD。
4.1 三角分解 LU
三角分解法是将原方阵分解成一个上三角形矩阵和一个下三角形矩阵，这样的分解法又称为 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
奇异值分解 (singular value decomposition, SVD) 是一种正交矩阵分解法。SVD 是最可靠的分解法，但是它比 QR 分解法要花上近十倍的计算时间。[U,S,V]=svd(A)，其中 U 和 V 分别代表两个正交矩阵，而 S 代表一对角矩阵。和 QR 分解法相同， 原矩阵 A 不必为正方矩阵。使用 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+005. 特殊矩阵
下面介绍的多种特殊矩阵，都是在 matrixcalc 库中提供的。
5.1 Hankel Matrix
汉克尔矩阵 (Hankel Matrix) 是具有恒定倾斜对角线的方形矩阵。Hankel 矩阵的行列式称为 catalecticant。该函数根据 n 向量 x 的值构造 n 阶 Hankel 矩阵。矩阵的每一行是前一行中值的循环移位。
矩阵定义：
> 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
希尔伯特矩阵是一种数学变换矩阵，正定，且高度病态（即，任何一个元素发生一点变动，整个矩阵的行列式的值和逆矩阵都会发生巨大变化），病态程度和阶数相关。希尔伯特矩阵是一种特殊的汉克尔矩阵，该函数返回 n 乘 n 希尔伯特矩阵。
矩阵定义：
> 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
创造矩阵，n 阶创建矩阵也称为推导矩阵，该函数返回阶数 n 创建矩阵，在主对角线下方的子对角线上具有序列 1,2，…，n-1 的方阵。
矩阵定义：
> 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’s approximation）是一条用来取n的阶乘的近似值的数学公式。一般来说，当 n 很大的时候，n 阶乘的计算量十分大，所以斯特林公式十分好用，而且，即使在 n 很小的时候，斯特林公式的取值已经十分准确。
斯特林矩阵（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 15.5 Pascal matrix
帕斯卡矩阵：由杨辉三角形表组成的矩阵称为帕斯卡（Pascal）矩阵。此函数返回 n 乘以 Pascal 矩阵。在数学中，尤其是矩阵理论和组合学，Pascal 矩阵是一个下三角矩阵，行中有二项式系数。通过对相同顺序的对称 Pascal 矩阵执行 LU 分解并返回下三角矩阵，可以容易地获得它。
帕斯卡的三角形是由数字行组成的三角形。第一行具有条目1.每个后续行通过添加前一行的相邻条目而形成，替换为 0，其中不存在相邻条目。pascal 函数通过选择与指定矩阵维度相对应的 Pascal 三角形部分来形成 Pascal 矩阵。
矩阵定义：
> 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 序列导出的 n + 1 平方 Fibonacci 矩阵。
计算公式：
> 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
复制矩阵，当 A 是对称矩阵时，该函数构造将 vech（A）映射到 vec（A）的线性变换 D。
计算公式：
> 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
通过 R 语言，我们实现了对于矩阵的各种计算操作，非常方便！有了好的工具，不管是学习还是应用，都会事半功倍。本文只是列举了矩阵的操作方法，没有解释计算的推到过程，推到过程请参考线性代码的教科书。
往期精彩
统计分析，数据处理
长按二维码.关注Leosay卫统