我有一个矩阵，其中包含几个国家的行业经济信息。我想做一系列取决于矩阵中元素的位置和相应的列/行名称的计算。
行/列名称包含三个字母的国家/地区代码，后跟行业编号。行名称和列名称相同。
计算如下：
我将举一个简单的例子，尽管我的数据比这个例子大得多。
假设我有矩阵a
set.seed(10)
a <- matrix(sample(36) , nrow = 6)
colnames(a) <- rownames(a) <- paste( rep(c("aaa" , "bbb" , "ccc") , each = 2), rep(c(1:2) , times = 3))
给予：
aaa 1 aaa 2 bbb 1 bbb 2 ccc 1 ccc 2
aaa 1 9 19 26 6 5 21
aaa 2 10 24 2 30 36 27
bbb 1 12 15 13 11 20 16
bbb 2 8 28 33 18 34 17
ccc 1 22 35 14 23 29 4
ccc 2 7 31 25 32 1 3
在这种情况下，有三个国家“aaa”、“bbb”和“ccc”，只有两个产业，1和2。我想计算在自己国家生产的产业的价值，并相互相乘。
例如，在a[1, 4]此将“aaa 1”与“bbb 2”匹配。我想将每个国家的行业组合 1 - 2 相乘，即“aaa 1” - “aaa 2”与“bbb 1” - “bbb 2”（19 * 11 = 209）。沿对角线块的值将被平方（同一国家/地区的两倍）。
最终矩阵如下所示：
aaa 1 aaa 2 bbb 1 bbb 2 ccc 1 ccc 2
aaa 1 81 361 117 209 261 76
aaa 2 100 576 330 432 10 72
bbb 1 117 209 169 121 377 44
bbb 2 330 432 1089 324 33 54
ccc 1 261 76 377 44 841 16
ccc 2 10 72 33 54 1 9
我使用以下低效的代码手动计算：
b <- kronecker(diag(1, 3), matrix(1, 2, 2))
b <- (a * b)^2
c <- b
c[1,3] <- 9 * 13
c[1,4] <- 19 * 11
c[1,5] <- 9 * 29
c[1,6] <- 19 * 4
c[2,3] <- 10 * 33
c[2,4] <- 24 * 18
c[2,5] <- 10 * 1
c[2,6] <- 24 * 3
.
.
.
c[6,1] <- 1 * 10
c[6,2] <- 3 * 24
c[6,3] <- 1 * 33
c[6,4] <- 3 * 18
是否有灵活的代码来计算此数据，无论数据中的国家和行业数量如何？我很感激任何帮助。
编辑：我正在使用的矩阵很大（尺寸超过 2000 * 2000）。我需要一个可以在不冻结 R 的情况下处理此类数据的代码。不幸的是， outer 可以解决此处针对较小的较小矩阵提出的问题，但在处理较大的矩阵时会陷入困境。
假设每个国家都有相同数量的行业，并且dimnames(a)排序良好，如 A1、A2...Am、B1、B2、...Bm。这种方法利用了向量化计算的优势，避免了用暗名索引矩阵，从而提高了效率并显着节省了内存。
m: 行业数量m <- 2
n <- nrow(a)
block <- matrix(a[kronecker(diag(n/m), matrix(1, m, m)) == 1], m)
block[rep(1:m, len = n), ] * c(block[, c(t(matrix(1:n, m)))])
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 81 361 117 209 261 76
# [2,] 100 576 330 432 10 72
# [3,] 117 209 169 121 377 44
# [4,] 330 432 1089 324 33 54
# [5,] 261 76 377 44 841 16
# [6,] 10 72 33 54 1 9
a如果你的实际矩阵小于10000，这样的东西应该表现得足够好。你可以进一步优化，特别是如果你知道它们是n排序的（例如，参见@DarrenTsai 的答案），但这并不是明显必要的n。ndimnames(a)
n <- nrow(a)
nms <- rownames(a)
u <- unlist(strsplit(nms, " "))
country <- u[c(TRUE, FALSE)]
industry <- u[c(FALSE, TRUE)]
i1 <- seq_len(n)
j2 <- rep(i1, each = n)
j1 <- match(paste(country, industry[j2]), nms)
i2 <- match(paste(country[j2], industry), nms)
matrix(a[cbind(i1, j1)] * a[cbind(i2, j2)], n, n)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 81 361 117 209 261 76
[2,] 100 576 330 432 10 72
[3,] 117 209 169 121 377 44
[4,] 330 432 1089 324 33 54
[5,] 261 76 377 44 841 16
[6,] 10 72 33 54 1 9
当然，如果dimnames(a)未排序，但您想获得排序的好处，那么您可以排列 的行和列a，对生成的块矩阵进行操作，并对答案的行和列进行逆排列。
n <- nrow(a)
nms <- rownames(a)
## Number of industries, countries:
ni <- sum(startsWith(nms, substr(nms[1L], 1L, 3L)))
nc <- n %/% ni
p <- order(nms)
a.pp <- a[p, p] # permute
## Essentially @DarrenTsai's answer {uglier but highly optimized}:
k <- sequence(rep.int(ni, n), rep(seq.int(1L, n, ni), each = ni) + seq.int(0L, n * n - 1L, n))
b <- matrix(a.pp[k], ni, n)
r.pp <- b[rep.int(seq_len(ni), nc), ] * as.vector(b[, sequence(rep.int(nc, ni), seq_len(ni), ni)])
q <- p
q[q] <- seq_along(q)
r <- r.pp[q, q] # inverse permute
实际上，这里有一些简洁的矩阵代数。Block [i,j]of是blocks和ofr.pp的元素乘积。因此可以看作 的一种块叉积，因为是 的对角块的水平串联。[i,i][j,j]a.ppr.ppbba.pp
另一种方法是将原始矩阵子集为所需的子矩阵并将它们相乘（GKi1）。
x <- do.call(rbind, strsplit(rownames(a), " "))
x <- apply(outer(unique(x[,1]), unique(x[,2]), paste), 1,
\(i) a[i, i], simplify = FALSE)
x <- do.call(rbind, lapply(x, \(y) do.call(cbind, lapply(x, `*`, y))))
rownames(x) <- colnames(x)
x
# aaa 1 aaa 2 bbb 1 bbb 2 ccc 1 ccc 2
#aaa 1 81 361 117 209 261 76
#aaa 2 100 576 330 432 10 72
#bbb 1 117 209 169 121 377 44
#bbb 2 330 432 1089 324 33 54
#ccc 1 261 76 377 44 841 16
#ccc 2 10 72 33 54 1 9
x[rownames(a), colnames(a)] #In case the order of the original matrix is needed
一种变体，以不同的方式提取原始矩阵的相关信息 - 类似于@Mikael Jagan（谢谢！）的评论，可能看起来像（GKi1B）。
n <- nrow(a)
x <- do.call(rbind, strsplit(rownames(a), " ", TRUE))
x1 <- unique(x[,1])
x2 <- unique(x[,2])
n1 <- length(x1)
n2 <- length(x2)
A <- a[rep(match(t(outer(unique(x[,1]), unique(x[,2]), paste)), rownames(a)), each=n2) +
sequence(rep.int(n2, n), sequence(rep.int(n2, n1), seq.int(0L, by=n*n2, length.out=n1), n-1L))]
do.call(rbind, lapply(seq.int(0, by=n2*n2, length.out=n1),
\(i) matrix(A * A[seq_len(n2*n2)+i], n2)) )
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] 81 361 117 209 261 76
#[2,] 100 576 330 432 10 72
#[3,] 117 209 169 121 377 44
#[4,] 330 432 1089 324 33 54
#[5,] 261 76 377 44 841 16
#[6,] 10 72 33 54 1 9
另一种使用outer但也可能适用于更大矩阵的方法可能是（GKi2）。
n <- rownames(a)
i <- do.call(rbind, strsplit(n, " "))
j <- outer(i[,1], i[,2], paste)
array(a[cbind(n, c(j))] * a[cbind(c(t(j)), n[col(a)])], dim(a))
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] 81 361 117 209 261 76
#[2,] 100 576 330 432 10 72
#[3,] 117 209 169 121 377 44
#[4,] 330 432 1089 324 33 54
#[5,] 261 76 377 44 841 16
#[6,] 10 72 33 54 1 9
基准
set.seed(10)
a <- matrix(sample(1:100, 2002*2002, TRUE) , nrow = 2002)
colnames(a) <- rownames(a) <- paste( rep(strrep(letters, 3) , each = 77), 1:77)
bench::mark(check=FALSE,
ThomasIsCoding = {p <- expand.grid(rep(list(row.names(a)), 2))
p1 <- as.matrix(transform(p, Var2 = paste(sub("\\d+", "", Var1), sub("\\D+", "", Var2), sep = "")))
p2 <- as.matrix(transform(p, Var1 = paste(sub("\\d+", "", Var2), sub("\\D+", "", Var1), sep = "")))
`dim<-`(a[p1] * a[p2], dim(a))},
"Mikael Jagan" = {n <- nrow(a)
nms <- rownames(a)
u <- unlist(strsplit(nms, " "))
country <- u[c(TRUE, FALSE)]
industry <- u[c(FALSE, TRUE)]
i1 <- seq_len(n)
j2 <- rep(i1, each = n)
j1 <- match(paste(country, industry[j2]), nms)
i2 <- match(paste(country[j2], industry), nms)
matrix(a[cbind(i1, j1)] * a[cbind(i2, j2)], n, n)},
GKi2 = {n <- rownames(a)
i <- do.call(rbind, strsplit(n, " "))
j <- outer(i[,1], i[,2], paste)
array(a[cbind(n, c(j))] * a[cbind(c(t(j)), n[col(a)])], dim(a))},
"Mikael & Darren" = {n <- nrow(a)
nms <- rownames(a)
ni <- sum(startsWith(nms, substr(nms[1L], 1L, 3L)))
nc <- n %/% ni
p <- order(nms)
a.pp <- a[p, p] # permute
b <- matrix(a.pp[as.logical(kronecker(diag(nc), matrix(1, ni, ni)))], ni, n)
r.pp <- b[rep.int(seq_len(ni), nc), ] * as.vector(b[, sequence(rep.int(nc, ni), seq_len(ni), ni)])
q <- p
q[q] <- seq_along(q)
r.pp[q, q]},
GKi1 = {x <- do.call(rbind, strsplit(rownames(a), " "))
x <- apply(outer(sort(unique(x[,1])), sort(unique(x[,2])), paste), 1,
\(i) a[i, i], simplify = FALSE)
x <- do.call(rbind, lapply(x, \(y) do.call(cbind, lapply(x, `*`, y))))
rownames(x) <- colnames(x)
x},
"Mikael & DarrenB" = {n <- nrow(a)
nms <- rownames(a)
ni <- sum(startsWith(nms, substr(nms[1L], 1L, 3L)))
nc <- n %/% ni
p <- order(nms)
a.pp <- a[p, p] # permute
b <- {k <- sequence(rep.int(ni, n), rep(seq.int(1L, n, ni), each = ni) + seq(0L, n * n - 1L, n)); matrix(a.pp[k], ni, n)}
r.pp <- b[rep.int(seq_len(ni), nc), ] * as.vector(b[, sequence(rep.int(nc, ni), seq_len(ni), ni)])
q <- p
q[q] <- seq_along(q)
r.pp[q, q]},
GKi1B = {n <- nrow(a)
x <- do.call(rbind, strsplit(rownames(a), " ", TRUE))
x1 <- unique(x[,1])
x2 <- unique(x[,2])
n1 <- length(x1)
n2 <- length(x2)
A <- a[rep(match(t(outer(unique(x[,1]), unique(x[,2]), paste)), rownames(a)), each=n2) +
sequence(rep.int(n2, n), sequence(rep.int(n2, n1), seq.int(0L, by=n*n2, length.out=n1), n-1L))]
do.call(rbind, lapply(seq.int(0, by=n2*n2, length.out=n1),
\(i) matrix(A * A[seq_len(n2*n2)+i], n2))) }
)
结果
expression min median `itr/sec` mem_alloc `gc/sec` n_itr n_gc
<bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl> <int> <dbl>
1 ThomasIsCoding 1.24s 1.24s 0.810 872.1MB 4.05 1 5
2 Mikael Jagan 743.23ms 743.23ms 1.35 367.1MB 1.35 1 1
3 GKi2 812.95ms 812.95ms 1.23 734.1MB 4.92 1 4
4 Mikael & Darren 62.56ms 64.66ms 15.5 140.7MB 9.69 8 5
5 GKi1 28.72ms 29.57ms 33.3 48.8MB 7.83 17 4
6 Mikael & DarrenB 28.21ms 32.37ms 30.9 48.9MB 5.79 16 3
7 GKi1B 26.7ms 28.3ms 35.0 49.8MB 9.71 18 5
在本例中，Mikael & DarrenB、GKi1 和 GKiB 速度最快且使用的内存量最少。