更新！

假设每个国家的产业数量相等，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

旧版本

comb <- function(x, y, mat) {
  prod(mat[outer(c(x[1], y[1]), c(x[2], y[2]), paste)])
}

nms <- strsplit(row.names(a), ' ')
outer(nms, nms, Vectorize(comb, c('x', 'y')), mat = a)

如果实际的矩阵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 also working for n = 0}:
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 <- r.pp[q, q] # inverse permute

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

您可以尝试下面的代码

> 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