学习 R 语言：矩阵和数组
目录
本文内容来自《R 语言编程艺术》(The Art of R Programming)，有部分修改
矩阵 (Matrix) 是一种特殊的向量，包含两个附加属性：行数和列数。
数组 (Array) 是更一般的对象，可以有多个维度。矩阵是二维数组。
创建矩阵
R 中下标从 1 开始，矩阵按列存储。
使用 matrix() 函数创建矩阵
y <- matrix(
c(1, 2, 3, 4),
nrow=2,
ncol=2,
)
print(y) [,1] [,2]
[1,] 1 3
[2,] 2 4
在给定所有元素的情况下，nrow 和 ncol 两个参数只需要给出一个
y <- matrix(
1:4,
nrow=2
)
print(y) [,1] [,2]
[1,] 1 3
[2,] 2 4
另一种方法是创建空的矩阵，再对每一位赋值
y <- matrix(nrow=2, ncol=2)
print(y) [,1] [,2]
[1,] NA NA
[2,] NA NA
y[1, 1] <- 1
y[2, 1] <- 2
y[1, 2] <- 3
y[2, 2] <- 4
print(y) [,1] [,2]
[1,] 1 3
[2,] 2 4
设置参数 byrow=TRUE，可以按行提供的矩阵数据。
但矩阵实际还是按列存储
m <- matrix(
1:6,
nrow=2,
byrow=TRUE,
)
print(m) [,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
一般矩阵运算
线性代数运算
矩阵乘法
print(y %*% y) [,1] [,2]
[1,] 7 15
[2,] 10 22
矩阵元素乘法
print(3*y) [,1] [,2]
[1,] 3 9
[2,] 6 12
矩阵加法
print(y + y) [,1] [,2]
[1,] 2 6
[2,] 4 8
矩阵索引
向量索引同样适用于矩阵索引
提取矩阵列
z <- matrix(
c(
1, 2, 3, 4,
1, 1, 0, 0,
1, 0, 1, 0
),
nrow=4,
)
print(z) [,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 1 0
[3,] 3 0 1
[4,] 4 0 0
print(z[, 2:3]) [,1] [,2]
[1,] 1 1
[2,] 1 0
[3,] 0 1
[4,] 0 0
提取矩阵行
y <- matrix(
c(
11, 21, 31,
12, 22, 32
),
ncol=2,
)
print(y) [,1] [,2]
[1,] 11 12
[2,] 21 22
[3,] 31 32
print(y[2:3,]) [,1] [,2]
[1,] 21 22
[2,] 31 32
print(y[2:3, 2])[1] 22 32
对子矩阵赋值
y <- matrix(
1:6,
nrow=3
)
print(y) [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
y[c(1,3),] <- matrix(
c(1, 1, 8, 12),
nrow=2
)
print(y) [,1] [,2]
[1,] 1 8
[2,] 2 5
[3,] 1 12
另一个示例
x <- matrix(nrow=3, ncol=3)
y <- matrix(
c(4, 5, 2, 3),
nrow=2
)
print(y) [,1] [,2]
[1,] 4 2
[2,] 5 3
y <- matrix(
1:6,
nrow=3
)
print(y) [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
print(y[-2,]) [,1] [,2]
[1,] 1 4
[2,] 3 6
扩展案例：图像操作
library(pixmap)
mtrush1 <- read.pnm("mtrush1.pgm")
mtrush1Pixmap image
Type : pixmapGrey
Size : 194x259
Resolution : 1x1
Bounding box : 0 0 259 194
plot(mtrush1)
str(mtrush1)Formal class 'pixmapGrey' [package "pixmap"] with 6 slots
..@ grey : num [1:194, 1:259] 0.278 0.263 0.239 0.212 0.192 ...
..@ channels: chr "grey"
..@ size : int [1:2] 194 259
..@ cellres : num [1:2] 1 1
..@ bbox : num [1:4] 0 0 259 194
..@ bbcent : logi FALSE
mtrush1@grey[28, 88]0.796078431372549
mtrush2 <- mtrush1
mtrush2@grey[84:163, 135:177] <- 1
plot(mtrush2)
blurpart <- function(img, rows, cols, q) {
lrows <- length(rows)
lcols <- length(cols)
newimg <- img
randomnoise <- matrix(nrow=lrows, ncol=lcols, runif(lrows*lcols))
newimg@grey[rows, cols] <- (1-q) * img@grey[rows, cols] + q * randomnoise
return(newimg)
}mtrush3 <- blurpart(
mtrush1,
84:163,
135:177,
0.2
)
plot(mtrush3)
矩阵元素筛选
x <- matrix(
c(1, 2, 3, 2, 3, 4),
nrow=3
)
print(x) [,1] [,2]
[1,] 1 2
[2,] 2 3
[3,] 3 4
print(x[x[,2] >= 3, ]) [,1] [,2]
[1,] 2 3
[2,] 3 4
分析上面的命令
j <- x[,2] >= 3
print(j)[1] FALSE TRUE TRUE
print(x[j,]) [,1] [,2]
[1,] 2 3
[2,] 3 4
注意：求 j 是向量化计算
另一个示例
m <- matrix(
1:6,
nrow=3
)
print(m) [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
z <- c(5, 12, 13)
print(m[z %% 2 == 1,]) [,1] [,2]
[1,] 1 4
[2,] 3 6
print(m[m[,1] > 1 & m[,2] > 5, ])[1] 3 6
分析
a <- m[,1] > 1
print(a)[1] FALSE TRUE TRUE
b <- m[,2] > 5
print(b)
[1] FALSE FALSE TRUE
c <- a & b
print(c)[1] FALSE FALSE TRUE
print(m[c, ])[1] 3 6
注意：这里使用 & 表示向量的逻辑与运算。而 && 是用于 if 语句的标量逻辑与运算
向量运算也适用于矩阵
m <- matrix(
c(5, 2, 9, -1, 10, 11),
nrow=3
)
print(m) [,1] [,2]
[1,] 5 -1
[2,] 2 10
[3,] 9 11
print(which(m > 2))[1] 1 3 5 6
扩展案例：生成协方差矩阵
row() 函数和 col() 函数
协方差矩阵是对称的，对角线元素均为 1
makecov <- function(rho, n) {
m <- matrix(nrow=n, ncol=n)
m <- ifelse(row(m) == col(m), 1, rho)
return(m)
}示例
print(makecov(0.2, 3)) [,1] [,2] [,3]
[1,] 1.0 0.2 0.2
[2,] 0.2 1.0 0.2
[3,] 0.2 0.2 1.0
解释
z <- matrix(
3:8,
nrow=3
)
print(z) [,1] [,2]
[1,] 3 6
[2,] 4 7
[3,] 5 8
print(row(z)) [,1] [,2]
[1,] 1 1
[2,] 2 2
[3,] 3 3
print(col(z)) [,1] [,2]
[1,] 1 2
[2,] 1 2
[3,] 1 2
print(row(z) == col(z)) [,1] [,2]
[1,] TRUE FALSE
[2,] FALSE TRUE
[3,] FALSE FALSE
对矩阵的行和列调用函数
*apply 系列函数：apply, tapply, lapply
apply() 函数在各行或各列上调用指定函数
使用 apply() 函数
一般形式
apply(m, dimcode, f, fargs)
其中
m：矩阵dimcode：维度编号，1 表示对每行应用函数，2 表示对每列应用函数f：应用在行或列上的函数，第一个参数为一行或一列fargs：可选参数，用于传递给f函数
z <- matrix(
1:6,
nrow=3
)
print(z) [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
对每一列应用函数 mean
注：可以直接使用 colmean() 函数
print(apply(z, 2, mean))[1] 2 5
对每一行应用自定义函数
f <- function(x) x/c(2, 8)
y <- apply(z, 1, f)
print(y) [,1] [,2] [,3]
[1,] 0.5 1.000 1.50
[2,] 0.5 0.625 0.75
调用的函数返回 k 个元素的向量，则 apply 函数返回的结果就有 k 行。
执行几次函数，则返回结果就有多少列。
如果函数返回一个标量，则 apply 函数返回向量。
t() 函数可以实现转置
print(t(apply(z, 1, f))) [,1] [,2]
[1,] 0.5 0.500
[2,] 1.0 0.625
[3,] 1.5 0.750
附加参数
0 和 1 构成的矩阵，对于每行，前 d 个元素中 1 较多则取 1，否则取 0
copymaj <- function(rw, d) {
maj <- sum(rw[1:d]) / d
return(if(maj > 0.5) 1 else 0)
}x <- matrix(
c(
1, 1, 1, 0,
0, 1, 0, 1,
1, 1, 0, 1,
1, 1, 1, 1,
0, 0, 1, 0
),
nrow=4
)
print(x) [,1] [,2] [,3] [,4] [,5]
[1,] 1 0 1 1 0
[2,] 1 1 1 1 0
[3,] 1 0 0 1 1
[4,] 0 1 1 1 0
print(apply(x, 1, copymaj, 3))[1] 1 1 0 1
print(apply(x, 1, copymaj, 2))[1] 0 1 0 0
扩展案例：寻找异常值
识别偏离最远的观测值
findols <- function(x) {
findol <- function(xrow) {
mdn <- median(xrow)
devs <- abs(xrow-mdn)
return(which.max(devs))
}
return(apply(x, 1, findol))
}增加或删除矩阵的行或列
矩阵无法改变大小，增加或删除实际上是创建新的对象
改变矩阵大小
重新赋值可以改变向量的大小
x <- c(12, 5, 13, 16, 8)
print(x)[1] 12 5 13 16 8
x <- c(x, 20)
print(x)[1] 12 5 13 16 8 20
x <- c(x[1:3], 20, x[4:6])
print(x)[1] 12 5 13 20 16 8 20
x <- x[-2:-4]
print(x)[1] 12 16 8 20
使用 rbind() 和 cbind() 为矩阵增加行或列
one <- c(1, 1, 1, 1)
z <- matrix(
c(
1, 2, 3 ,4,
1, 1, 0, 0,
1, 0, 1, 0
),
nrow=4
)
print(z) [,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 1 0
[3,] 3 0 1
[4,] 4 0 0
print(cbind(one, z)) one
[1,] 1 1 1 1
[2,] 1 2 1 0
[3,] 1 3 0 1
[4,] 1 4 0 0
可能会用到循环补齐 (recycling)
print(cbind(1, z)) [,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 2 1 0
[3,] 1 3 0 1
[4,] 1 4 0 0
cbind() 和 rbind() 可以用来快速生成小矩阵
q <- cbind(
c(1, 2),
c(3, 4)
)
print(q) [,1] [,2]
[1,] 1 3
[2,] 2 4
注意：创建向量和矩阵均耗费时间，请谨慎使用 cbind 和 rbind
循环添加列或行的方式不可取，应该首先生成一个大矩阵，然后再往里填值。
可以通过重新赋值删除矩阵的行或列
m <- matrix(1:6, nrow=3)
print(m) [,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
m <- m[c(1,3),]
print(m) [,1] [,2]
[1,] 1 4
[2,] 3 6
扩展案例：找到图中距离最近的一对端点
矩阵元素 (i, j) 表示城市 i 和 城市 j 间的距离
imin 函数寻找当前行的最小值，因为矩阵是对称的，所以仅搜索从行号 i + 1 开始到结尾。
返回最小值的索引（即矩阵中的列数）和最小值
imin <- function(x) {
lx <- length(x)
i <- x[lx]
j <- which.min(x[(i+1):(lx-1)])
k <- i + j
return(c(k, x[k]))
}mind 函数为每行分别求最小值 (wmins)，再找出全局最小值，返回最小值，行号和列号
mind <- function(d) {
n <- nrow(d)
dd <- cbind(d, 1:n)
wmins <- apply(dd[-n,], 1, imin)
i <- which.min(wmins[2,])
j <- wmins[1, i]
return(c(d[i, j], i, j))
}示例
q <- matrix(
c(
0, 12, 13, 8, 20,
12, 0, 15, 28, 88,
13, 15, 0, 6, 9,
8, 28, 6, 0, 33,
20, 88, 9, 33, 0
),
nrow=5,
byrow=T,
)
print(q) [,1] [,2] [,3] [,4] [,5]
[1,] 0 12 13 8 20
[2,] 12 0 15 28 88
[3,] 13 15 0 6 9
[4,] 8 28 6 0 33
[5,] 20 88 9 33 0
print(mind(q))[1] 6 3 4
分析
增加记录当前行号的一列
n <- nrow(q)
qq <- cbind(q, 1:n)
print(qq) [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 12 13 8 20 1
[2,] 12 0 15 28 88 2
[3,] 13 15 0 6 9 3
[4,] 8 28 6 0 33 4
[5,] 20 88 9 33 0 5
为前 n-1 行，即前 4 行，寻找最小值，并记录最小值（第 2 行）及其对应的列数（第 1 行）。
返回结果中每列代表矩阵中的行，即第 1 列表示矩阵的第 1 行，最后一列表示矩阵的第 4 行
wmins <- apply(qq[-n,], 1, imin)
print(wmins) [,1] [,2] [,3] [,4]
[1,] 4 3 4 5
[2,] 8 15 6 33
求最小值，返回列号，即矩阵中的行号 i
i <- which.min(wmins[2,])
print(i)[1] 3
从 wmins 中取出矩阵的列号 j
j <- wmins[1, i]
print(j)[1] 4
print(c(qq[i, j], i, j))[1] 6 3 4
如果最小值唯一，则有更简单的方法
注：有两次矩阵计算，效率不是最优的
minda <- function(d) {
smallest <- min(d)
ij <- which(d == smallest, arr.ind=TRUE)
return(c(smallest, ij))
}示例，需要将 0 值和重复值都替换为比较大的数
p = q
p[row(p)>=col(p)] = 99999
print(p) [,1] [,2] [,3] [,4] [,5]
[1,] 99999 12 13 8 20
[2,] 99999 99999 15 28 88
[3,] 99999 99999 99999 6 9
[4,] 99999 99999 99999 99999 33
[5,] 99999 99999 99999 99999 99999
print(minda(p))[1] 6 3 4
向量与矩阵的差异
矩阵实际上就是向量，只不过多了两个属性：行数和列数
z <- matrix(1:8, nrow=4)
print(z) [,1] [,2]
[1,] 1 5
[2,] 2 6
[3,] 3 7
[4,] 4 8
求长度
print(length(z))[1] 8
但 z 不仅仅是向量
print(class(z))[1] "matrix" "array"
print(attributes(z))$dim
[1] 4 2
用 dim() 函数访问 dim 属性
print(dim(z))[1] 4 3
nrow() 和 ncol() 是 dim() 的简单封装
print(nrow(z))[1] 4
print(ncol(z))
[1] 2
nrowfunction (x)
dim(x)[1L]
避免意外降维
从矩阵中提取一行，会减少一维，即返回向量
z <- matrix(1:8, nrow=4)
print(z) [,1] [,2]
[1,] 1 5
[2,] 2 6
[3,] 3 7
[4,] 4 8
r <- z[2,]
print(r)[1] 2 6
print(attributes(z))$dim
[1] 3 2
print(attributes(r))NULL
str(z) int [1:3, 1:2] 1 2 3 4 5 6
str(r) int [1:2] 2 6
使用 drop 参数禁止自动减少维度
r <- z[2,, drop=FALSE]
print(r) [,1] [,2]
[1,] 2 5
print(dim(r))[1] 1 2
可见中括号 [ 实际上也是一个函数
print(z[3, 2])[1] 6
print("["(z, 3, 2))[1] 6
as.matrix() 函数可以将向量转为矩阵
u <- 1:3
print(u)[1] 1 2 3
v <- as.matrix(u)
print(v) [,1]
[1,] 1
[2,] 2
[3,] 3
print(attributes(u))NULL
print(attributes(v))$dim
[1] 3 1
矩阵的行和列的命名问题
矩阵的行和列也可以命名
z <- matrix(1:4, nrow=2)
print(z) [,1] [,2]
[1,] 1 3
[2,] 2 4
print(colnames(z))NULL
colnames(z) <- c("a", "b")
print(z) a b
[1,] 1 3
[2,] 2 4
print(colnames(z))[1] "a" "b"
print(z[, "a"])[1] 1 2
高维数组
三个学生，两次考试，每次考试记录两个成绩
第一次考试的数据，行表示学生，列表示单次考试的两个成绩
firsttest <- matrix(
c(
46, 30,
21, 25,
50, 50
),
nrow=3,
byrow=T
)
print(firsttest) [,1] [,2]
[1,] 46 30
[2,] 21 25
[3,] 50 50
第二次考试的成绩
secondtest <- matrix(
c(
46, 43,
41, 35,
50, 50
),
nrow=3,
byrow=T
)
print(secondtest) [,1] [,2]
[1,] 46 43
[2,] 41 35
[3,] 50 50
合并为一个数组，两次考试形成第三维，叫做层 (layer)
tests <- array(
data=c(firsttest, secondtest),
dim=c(3, 2, 2)
)print(attributes(tests))$dim
[1] 3 2 2
学生 3 (dim 1) 在第一次考试 (dim 3) 的第二个部分 (dim 2) 中的得分
print(tests[3, 2, 1])[1] 50
print(tests), , 1
[,1] [,2]
[1,] 46 30
[2,] 21 25
[3,] 50 50
, , 2
[,1] [,2]
[1,] 46 43
[2,] 41 35
[3,] 50 50
参考
本文的 Jupyter Notebook 代码请访问如下项目：
