發現好久沒有更新微信文了, 所謂才思枯竭, 黔驢技窮就是我現在的狀態. 記得看過這樣一句話: "如果你不知道寫什么東西, 那就寫不知道寫什么事情這件事吧". 深得我心.
分享一篇我CSND博客里面的R語言矩陣操作, 可以通過編程理解很多線性代數的概念. 這篇文章閱讀量2萬+, 而我的CSND博客閱讀量才10萬+, 可以看出博客的閱讀量分布不是正態的, 符合馬太效應.1.1 矩陣的生成生成一個4行4列的矩陣，這里用1~16數字。 mat <- matrix(1:16,4,4)
mat1.2 提取主對角線diag(mat)
1
6
11
16
1.3 生成對角線為1的對角矩陣m1 <- diag(4)
m11.4 提取矩陣的下三角mat[lower.tri(mat)]
2
3
4
7
8
12
1.5 提取矩陣上三角mat[upper.tri(mat)]
5
9
10
13
14
15
1.6 以矩陣下三角構建對角矩陣mat1 <- mat
mat1[upper.tri(mat1)] <- t(mat1)[upper.tri(mat1)]原矩陣mat： mat變換后的對角矩陣 mat11.7 將矩陣轉化為行列形式原矩陣，生成三列：行，列，值 mat相關代碼 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.9 將矩陣轉置t(mat)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 - B2.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)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假設B為3X3矩陣 B <- matrix(c(1,0,2,0,1,4,2,4,1),3,3)
B則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)
AB <- matrix(c(1,0,6,1),2,2)
Br1 <- 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
|
