专栏首页育种数据分析之放飞自我R语言中矩阵常用的操作(笔记)

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 条评论
登录 后参与评论

相关文章

  • 线性混合模型系列四:矩阵求解

    这篇文章通过R语言代码的形式,介绍给定方差组分的情况下,如何根据两种矩阵求解的方法分别计算BLUE值和BLUP值。

    邓飞
  • 一步法中混合线性模型方程组构建和控制--blupf90

    http://nce.ads.uga.edu/wiki/lib/exe/fetch.php?media=singlestepblupf90.pdf

    邓飞
  • 软件介绍: BLUPF90的无敌和寂寞

    前同桌说, 要带粉丝了, 一个换一个, 我就把题目改了一下, 原来的题目是《遗传育种软件三剑客之一:BLUPF90》,但我还是喜欢现在这个题目《BLUPF90的...

    邓飞
  • 彻底理解Java IO

    谈到IO,我们会想到从磁盘读取的文件IO,网络请求的Socket IO,还有可能我们不怎么常用的跨进程通信的管道IO...... 这些在Java中都被抽象为“...

    三好码农
  • MYSQL 小内存, 大问题

    每种数据库都有自己的管理内存的方法,MYSQL 管理内存(仅仅讨论 INNODB 数据库引擎)的方法大部分都关注在 innodb_buffer_pool_siz...

    AustinDatabases
  • 为什么MySQL内存占用这么大? for InnoDB

    这是 Innodb 引擎最重要的缓存,也是提升查询性能的重要手段。一般是global共享内存中占用最大的部分。在进行 SQL 读和写的操作时,首先并不是对物理数...

    elontian田凌翔
  • 如何实现一个会讲笑话的Python程序?

    笑话从哪里来?自己写肯定是不现实的。在这个“云”的时代,各种云都有,自然是不缺开放API的(大部分都是免费的)。随意一搜,果然被我找到一个接口:易源_笑话大全h...

    小小科
  • 【调试】939- 5个Chrome调试混合应用的技巧

    对前端开发人员来说,Chrome 真是一个必备的开发工具,大到页面展示,小到 BUG 调试/HTTP 抓包等,本文我将和大家分享自己做混合应用开发过程中经常用到...

    pingan8787
  • 前端优化

    Pack Components into a Multipart Document

    只喝牛奶的杀手
  • 这个功能下班之前必须做好

    1. 老板:“这个功能下班之前必须做好!”程序猿:“好的。” 第二天早上老板:“这个功能怎么还没做好?”程序猿:“我还没下班呢。” 2. 世界上最遥远的距离不是...

    程序员互动联盟

扫码关注云+社区

领取腾讯云代金券