如何用矩阵运算或内置函数实现两矩阵对应行的叉积计算

原始数据与需求 #

首先定义两个输入矩阵：

A <- matrix(nrow = 5, ncol = 3, byrow = TRUE, c(
  5, 1, 2, 
  3, 1, 5, 
  5, 1, 3, 
  7, 2, 8, 
  2, 4, 2
))

B <- matrix(nrow = 5, ncol = 2, byrow = TRUE, c(
  2, 1,
  3, 5,
  1, 4,
  3, 3,
  3, 4
))

需求是创建一个三维数组arr，维度为c(ncol(A), ncol(B), nrow(A))，其中第i个切片对应A的第i行与B的第i行的外积（即A[i, ] %*% t(B[i, ])）。例如第2切片的预期结果为：

A[2, ] %*% t(B[2, ])
#      [,1] [,2]
# [1,]    9   15
# [2,]    3    5
# [3,]   15   25

现有实现方式 #

你已经通过两种方式实现需求：

For循环实现 #

arr <- array(NA, dim = c(ncol(A), ncol(B), nrow(A))) 
for (i in 1:nrow(A)) {
  arr[, , i] <- A[i, ] %*% t(B[i, ])
}

sapply函数实现 #

temp <- sapply(1:nrow(A), function(i) A[i, ] %*% t(B[i, ]), simplify = FALSE)
arr <- array(unlist(temp), dim = c(ncol(A), ncol(B), nrow(A))) 

更高效/简洁的实现方法 #

方法1：向量化矩阵广播操作 #

利用R的矩阵广播特性，通过元素级乘法直接生成所有外积元素，再调整维度得到目标数组，避免显式循环：

# 生成元素级乘积矩阵，再重塑为三维数组并调整维度顺序
arr <- aperm(
  array(
    A[, rep(1:ncol(A), each = ncol(B))] * B[, rep(1:ncol(B), ncol(A))],
    dim = c(nrow(A), ncol(B), ncol(A))
  ),
  c(3, 2, 1)
)

方法2：结合Map与tcrossprod #

Map可以对两组行数据逐对执行外积操作，再通过unlist和array重塑结果，代码更简洁：

arr <- array(
  unlist(Map(tcrossprod, split(A, 1:nrow(A)), split(B, 1:nrow(B)))),
  dim = c(ncol(A), ncol(B), nrow(A))
)

方法3：apply+tcrossprod #

通过apply按行处理合并后的A和B，逐行计算外积，最后调整维度：

arr <- apply(
  cbind(A, B), 1,
  function(x) tcrossprod(x[1:ncol(A)], x[(ncol(A)+1):(ncol(A)+ncol(B))])
)
dim(arr) <- c(ncol(A), ncol(B), nrow(A))

这些方法中，向量化广播操作在大数据量下的性能最优，因为它完全利用了R的底层向量化计算能力，避免了循环的开销。

内容的提问来源于stack exchange，提问作者little_sky