R语言 匹配行和列名称中包含的数字的矩阵计算

xvw2m8pv  于 2023-03-05  发布在  其他
关注(0)|答案(4)|浏览(145)

我有一个矩阵，其中包含几个国家的工业经济信息。我想根据矩阵中元素的位置和相应的列/行名称进行一系列计算。
行/列名称具有三个字母的国家/地区代码，后跟行业编号。行和列名称相同。
计算如下：
1.对于行中的每个元素，都将存在匹配的国家/地区和行业
1.保持产业组合不变，在比赛中把纯粹同一个国家的产业与另一个国家的同一个产业相乘。
我将给予一个简单的例子，尽管我的数据比这个例子大得多。
假设有矩阵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可以解决这里提出的较小矩阵的问题，但在处理较大矩阵时会卡住。

r

4条答案

按热度按时间
fae0ux8s

fae0ux8s1#

如果实际的矩阵a是n-x-n，且n小于10000，那么类似这样的方法应该可以很好地执行。您可以进一步优化，特别是如果您知道dimnames(a)是排序的（例如，参见@DarrenTsai的answer），但这显然不是必要的。

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

实际上，这里有一些简洁的矩阵代数，r.pp的块[i,j]是[i,i]和a.pp的块[j,j]的元素积，因此r.pp可以看作是b的一种块叉积，因为b是a.pp的对角块的水平连接。

赞(0）分享  回复(0）举报  2023-03-05
qlckcl4x

qlckcl4x2#

假设每个国家的产业数量相等，dimnames(a)的排序为A1，A2...Am，B1，B2，... Bm，该方法利用了向量化计算的优势，避免了用dimname索引矩阵，从而大大提高了效率并节省了内存。

  • 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
赞(0）分享  回复(0）举报  2023-03-05
ejk8hzay

ejk8hzay3#

另一种方式是将原始矩阵子集化为所需的子矩阵并将它们相乘（GKi 1）。

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的评论（谢谢！），可能看起来像（GKi 1B）。

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但也可以使用更大矩阵的方法可以是（GKi 2）。

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、GKi 1和GKiB速度最快，使用的内存量最低。

赞(0）分享  回复(0）举报  2023-03-05
b09cbbtk

b09cbbtk4#

您可以尝试下面的代码

> 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))
     [,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
赞(0）分享  回复(0）举报  2023-03-05
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com