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)
  1. 1
  2. 6
  3. 11
  4. 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)]
  1. 2
  2. 3
  3. 4
  4. 7
  5. 8
  6. 12

1.5 提取矩阵上三角

mat[upper.tri(mat)]
  1. 5
  2. 9
  3. 10
  4. 13
  5. 14
  6. 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

欢迎关注我的公众号:R-breeding

本文分享自微信公众号 - 育种数据分析之放飞自我(R-breeding)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-06-22

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券

年度创作总结 领取年终奖励