迭代矩阵计算场景下的R语言运算加速方案寻求

我在优化例程中发现生成res的代码是最大性能瓶颈，当前代码如下：

res <- exp(X %*% log(pr.t) + mX %*% log(1 - pr.t)) %*% wts

其中矩阵X元素固定（迭代中不变化），mX = 1 - X可预先存储复用；仅pr.t会在每次迭代时更新。我尝试过Rcpp，但性能和当前R代码相当，希望找到更高效的优化方案。

可复现示例代码：

X <- matrix(sample(c(0,1), 5000, replace = TRUE), 1000, 5)
mX <- 1 - X
pr.t <- matrix(runif(75), 5, 15)
wts <- runif(15)
res <- exp(X %*% log(pr.t) + mX %*% log(1 - pr.t)) %*% wts

优化方案1：数学公式简化，减少矩阵乘法次数 #

原公式中两次矩阵乘法可合并为一次，利用代数变形减少计算量：
因为 X %*% log(pr.t) + (1-X) %*% log(1-pr.t) 等价于 X %*% (log(pr.t) - log(1-pr.t)) + matrix(colSums(log(1-pr.t)), nrow(X), ncol(pr.t), byrow=TRUE)，这样只需要一次矩阵乘法，替代原来的两次：

# 预计算X（固定）
X_fixed <- X

# 每次迭代时执行
log_pr <- log(pr.t)
log_1pr <- log(1 - pr.t)
delta <- log_pr - log_1pr
log_term <- X_fixed %*% delta + matrix(colSums(log_1pr), nrow(X_fixed), ncol(pr.t), byrow = TRUE)
res <- exp(log_term) %*% wts

优化方案2：利用稀疏矩阵加速（若X中0元素占比高） #

如果X是稀疏的（大量0元素），将其转换为稀疏矩阵可以跳过0元素的计算，大幅提升矩阵乘法效率：

library(Matrix)
# 预计算稀疏矩阵（固定）
X_sparse <- Matrix(X, sparse = TRUE)

# 每次迭代时执行
log_pr <- log(pr.t)
log_1pr <- log(1 - pr.t)
log_term <- X_sparse %*% log_pr + (1 - X_sparse) %*% log_1pr
res <- exp(log_term) %*% wts

优化方案3：合并重复行减少计算量 #

如果X存在大量重复的行，可以预计算唯一行及其索引，每次迭代仅对唯一行计算，再扩展结果：

# 预计算唯一行和索引（固定）
unique_X <- unique(X)
row_map <- match(apply(X, 1, paste, collapse = ","), apply(unique_X, 1, paste, collapse = ","))

# 每次迭代时执行
log_pr <- log(pr.t)
log_1pr <- log(1 - pr.t)
unique_log_term <- unique_X %*% log_pr + (1 - unique_X) %*% log_1pr
log_term <- unique_log_term[row_map, ]
res <- exp(log_term) %*% wts

优化方案4：优化Rcpp实现（使用Armadillo库） #

之前的Rcpp性能不佳可能是未利用高效线性代数库，改用Armadillo实现矩阵运算，其底层优化的BLAS/LAPACK调用通常比R base更快：

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
arma::vec compute_res_arma(const arma::mat& X, const arma::mat& pr_t, const arma::vec& wts) {
  arma::mat log_pr = arma::log(pr_t);
  arma::mat log_1pr = arma::log(1.0 - pr_t);
  arma::mat log_term = X * log_pr + (1.0 - X) * log_1pr;
  arma::mat exp_term = arma::exp(log_term);
  arma::vec res = exp_term * wts;
  return res;
}

在R中调用：

# 预编译后，每次迭代调用
res <- compute_res_arma(X, pr.t, wts)

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