【原】R语言中矩阵常用的操作（笔记）

发现好久没有更新微信文了, 所谓才思枯竭, 黔驴技穷就是我现在的状态. 记得看过这样一句话: "如果你不知道写什么东西, 那就写不知道写什么事情这件事吧". 深得我心.

分享一篇我CSND博客里面的R语言矩阵操作, 可以通过编程理解很多线性代数的概念. 这篇文章阅读量2万+, 而我的CSND博客阅读量才10万+, 可以看出博客的阅读量分布不是正态的, 符合马太效应.

1.1 矩阵的生成

生成一个4行4列的矩阵，这里用1~16数字。

mat <- matrix(1:16,4,4)

mat

1	5	9	13
2	6	10	14
3	7	11	15
4	8	12	16

1.2 提取主对角线

diag(mat)

2. 1




4. 6




6. 11




8. 16

1.3 生成对角线为1的对角矩阵

m1 <- diag(4)

m1

1	0	0	0
0	1	0	0
0	0	1	0
0	0	0	1

1.4 提取矩阵的下三角

mat[lower.tri(mat)]

2. 2




4. 3




6. 4




8. 7




10. 8




12. 12

1.5 提取矩阵上三角

mat[upper.tri(mat)]

2. 5




4. 9




6. 10




8. 13




10. 14




12. 15

1.6 以矩阵下三角构建对角矩阵

mat1 <- mat

mat1[upper.tri(mat1)] <- t(mat1)[upper.tri(mat1)]

原矩阵mat：

mat

1	5	9	13
2	6	10	14
3	7	11	15
4	8	12	16

变换后的对角矩阵

mat1

1	2	3	4
2	6	7	8
3	7	11	12
4	8	12	16

1.7 将矩阵转化为行列形式

原矩阵，生成三列：行，列，值

mat

1	5	9	13
2	6	10	14
3	7	11	15
4	8	12	16

相关代码

nrow <- dim(mat)[1]

ncol <- dim(mat)[2]

row <- rep(1:nrow,ncol)

col <- rep(1:ncol, each=nrow)

frame <- data.frame(row,col,value =as.numeric(mat))

frame

row	col	value
1	1	1
2	1	2
3	1	3
4	1	4
1	2	5
2	2	6
3	2	7
4	2	8
1	3	9
2	3	10
3	3	11
4	3	12
1	4	13
2	4	14
3	4	15
4	4	16

1.8 将三列形式转化为矩阵

nrow <- max(frame[, 1])

ncol <- max(frame[, 2])

y <- rep(0, nrow * ncol)

y[(frame[, 2] - 1) * nrow + frame[, 1]] <- frame[, 3]

y[(frame[, 1] - 1) * nrow + frame[, 2]] <- frame[, 3]

matrix(y, nrow = nrow, ncol = ncol, byrow = T)

1	5	9	13
2	6	10	14
3	7	11	15
4	8	12	16

1.9 将矩阵转置

t(mat)

1	2	3	4
5	6	7	8
9	10	11	12
13	14	15	16

2.1 矩阵相加减

A=B=matrix(1:16,nrow=4,ncol=4)

A + B

2	10	18	26
4	12	20	28
6	14	22	30
8	16	24	32

A - B

0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0

2.2 数与矩阵相乘

c <- 2c*A

2	10	18	26
4	12	20	28
6	14	22	30
8	16	24	32

3.3 矩阵相乘

A 为m × n矩阵，B为n× k矩阵，用符合“%*%”

A <- matrix(1:12,3,4)

B <- matrix(1:20,4,5)

A%*%B

70	158	246	334	422
80	184	288	392	496
90	210	330	450	570

3.4 计算t(A)%*%B的方法

第一种，直接计算

A <- matrix(1:12,3,4)

B <- matrix(1:15,3,5)

t(A)%*%B

14	32	50	68	86
32	77	122	167	212
50	122	194	266	338
68	167	266	365	464

第二种方法，用crossprod函数，数据量大时效率更高

A <- matrix(1:12,3,4)

B <- matrix(1:15,3,5)

crossprod(A,B)

14	32	50	68	86
32	77	122	167	212
50	122	194	266	338
68	167	266	365	464

3.5 矩阵求逆

a <- matrix(rnorm(16),4,4)

solve(a)

-3.542393	5.8825038	-3.2421870	6.9619170
1.081745	-2.2446318	1.4850549	-2.0828270
-1.577580	2.4698567	-0.7070850	2.5241525
-0.830685	0.5105919	-0.3352182	0.5344842

矩阵与其逆矩阵的乘积为对角矩阵

round(solve(a)%*%a)

1	0	0	0
0	1	0	0
0	0	1	0
0	0	0	1

3.6 矩阵的广义逆矩阵

对于奇异阵，并不存在逆矩阵，但是可以计算其广义逆矩阵

a <- matrix(1:16,4,4)

solve(a)Error in solve.default(a): Lapack例行程序dgesv: 系统正好是奇异的: U[3,3] = 0
Traceback:


1. solve(a)

2. solve.default(a)

显示矩阵奇异，这里可以使用MASS包的ginv计算其广义逆矩阵

library(MASS)

a <- matrix(1:16,4,4)

ginv(a)

-0.285	-0.1075	0.07	0.2475
-0.145	-0.0525	0.04	0.1325
-0.005	0.0025	0.01	0.0175
0.135	0.0575	-0.02	-0.0975

3.7 矩阵的直积（Kronecker，克罗内克积），使用函数kronecker计算

A 与B的直积：LaTex写作 “A \bigotimes B”

假设A为2X2矩阵

A <- matrix(c(10,5,5,20),2,2)

A

10	5
5	20

假设B为3X3矩阵

B <- matrix(c(1,0,2,0,1,4,2,4,1),3,3)

B

1	0	2
0	1	4
2	4	1

则A和B的直积就是6X6的矩阵

kronecker(A,B)

10	0	20	5	0	10
0	10	40	0	5	20
20	40	10	10	20	5
5	0	10	20	0	40
0	5	20	0	20	80
10	20	5	40	80	20

3.8 矩阵的直和（direct sum）

公式：$ A\oplus B$，在LaTex中是 “A \oplus B “

A <- matrix(c(1,2,3,3,2,1),2,3)

A

1	3	2
2	3	1

B <- matrix(c(1,0,6,1),2,2)

B

1	6
0	1

r1 <- dim(A)[1];c1 <- dim(A)[2]

r2 <- dim(B)[1];c2 <- dim(B)[2]

direct_sum <- rbind(cbind(A,matrix(0,r2,c2)),cbind(matrix(0,r1,c1),B))

direct_sum

1	3	2	0	0
2	3	1	0	0
0	0	0	1	6
0	0	0	0	1