发现好久没有更新微信文了, 所谓才思枯竭, 黔驴技穷就是我现在的状态. 记得看过这样一句话: "如果你不知道写什么东西, 那就写不知道写什么事情这件事吧". 深得我心.
分享一篇我CSND博客里面的R语言矩阵操作, 可以通过编程理解很多线性代数的概念. 这篇文章阅读量2万+, 而我的CSND博客阅读量才10万+, 可以看出博客的阅读量分布不是正态的, 符合马太效应.
生成一个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 |
diag(mat)
m1 <- diag(4)
m1
1 | 0 | 0 | 0 |
---|---|---|---|
0 | 1 | 0 | 0 |
0 | 0 | 1 | 0 |
0 | 0 | 0 | 1 |
mat[lower.tri(mat)]
mat[upper.tri(mat)]
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 |
原矩阵,生成三列:行,列,值
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 |
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 |
t(mat)
1 | 2 | 3 | 4 |
---|---|---|---|
5 | 6 | 7 | 8 |
9 | 10 | 11 | 12 |
13 | 14 | 15 | 16 |
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 |
c <- 2c*A
2 | 10 | 18 | 26 |
---|---|---|---|
4 | 12 | 20 | 28 |
6 | 14 | 22 | 30 |
8 | 16 | 24 | 32 |
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 |
第一种,直接计算
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 |
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 |
对于奇异阵,并不存在逆矩阵,但是可以计算其广义逆矩阵
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 |
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 |
公式:$ 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