R语言中矩阵乘积替换求和为求积的高效实现方法

我懂你这种痛点——当矩阵规模上去之后，嵌套循环实现「行和列对应元素相乘再求积」的操作慢得让人抓狂，而默认的%*%只支持求和，没法直接拿来用。下面给你几个高效的替代方案，速度能赶上%*%的级别：

1. 对数转换法（最快，推荐正元素矩阵） #

这个方法利用了对数的数学性质：prod(x * y) = exp(sum(log(x) + log(y)))，完美借助%*%的底层优化实现加速，几乎和矩阵叉积一样快。

代码示例 #

set.seed(1)
a <- matrix(sample((1:100)/100, 15), ncol=3)
b <- matrix(sample((1:100)/100, 15), ncol=5)

# 对数转换+叉积+指数还原
log_a <- log(a)
log_b <- log(b)
result <- exp(log_a %*% log_b)

验证正确性 #

和你的嵌套循环结果对比：

# 原循环实现的基准结果
loop_result <- do.call('cbind', lapply(1:nrow(a), function(i) { 
  sapply(1:ncol(b), function(j) { 
    prod(a[i,] * b[,j]) 
  }) 
}))

all.equal(result, loop_result) # 返回TRUE，说明结果一致

注意：这个方法只适用于矩阵元素全为正数的情况，如果有0或负数，对数转换会报错，需要用其他方案。

2. Rcpp自定义函数（通用型，超大矩阵首选） #

如果你的矩阵包含非正数，或者追求极致速度，可以用Rcpp写C++层面的循环，避免R语言的循环开销，速度比R嵌套循环快几个数量级。

代码示例 #

library(Rcpp)
# 编译C++函数
cppFunction('NumericMatrix matrix_prod(NumericMatrix a, NumericMatrix b) {
  int n = a.nrow();    // a的行数
  int p = a.ncol();    // 公共维度
  int d = b.ncol();    // b的列数
  NumericMatrix result(n, d);
  
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < d; ++j) {
      double prod_val = 1.0;
      for (int k = 0; k < p; ++k) {
        prod_val *= a(i, k) * b(k, j);
      }
      result(i, j) = prod_val;
    }
  }
  return result;
}')

# 使用函数
result_cpp <- matrix_prod(a, b)
all.equal(result_cpp, loop_result) # 返回TRUE

3. 向量化apply优化（比循环快，但不如前两种） #

如果不想引入Rcpp，也可以用apply的向量化写法替代嵌套循环，虽然速度不如前两个方案，但比原生嵌套循环快很多。

代码示例 #

result_apply <- t(apply(a, 1, function(row) {
  apply(b, 2, function(col) {
    prod(row * col)
  })
}))

all.equal(result_apply, loop_result) # 返回TRUE

总结 #

• 若矩阵元素全为正数：对数转换法是最优选择，速度和%*%持平；
• 若有非正数或需要极致性能：Rcpp自定义函数是首选；
• 快速实现且不想依赖其他包：用向量化apply方案。

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