R中两个二维数组的for循环多维向量化实现问题

问题背景 #

需要替代遍历矩阵对(A,W)的for循环以提升运算速度，原循环逻辑为：

for(current.couple in 1:nrow(couples)){
  a_current <- current.rows[-which(current.rows == couples$current[current.couple])]
  a_candidate <- couples$candidate[current.couple]
  
  A <- A.use[ c(a_current, a_candidate),]
  W <- W.use[ c(a_current, a_candidate), c(a_current, a_candidate)]
  couples$D[ current.couple] <- crit.func(A, W)
}

思路是将所有版本的A、W堆叠为3D数组，以第三维度作为循环维度。示例数据与堆叠代码如下（注意：原W_append堆叠代码存在笔误，已修正）：

# 示例A矩阵
A1 <- matrix(c(2.4, 5.2, 8.4, 3.1, 6.05, 9.25), nrow = 2,ncol = 3, byrow = TRUE)
A2 <- matrix(c(4.5, 7.5, 10.5, 3.2, 6.2, 9.2), nrow = 2, ncol = 3, byrow = TRUE)
A3 <- matrix(c(2.1, 5, 8.2, 3.05, 6.02, 9.1), nrow = 2,ncol = 3, byrow = TRUE)
A4 <- matrix(c(4.12, 7.31, 10.3, 3.23, 6.1, 9), nrow = 2, ncol = 3, byrow = TRUE)

# 示例W矩阵
W1 <- matrix(c(1, 4, 2, 5), nrow = 2, ncol = 2, byrow = TRUE)
W2 <- matrix(c(9, 6, 8, 5), nrow = 2, ncol = 2, byrow = TRUE)
W3 <- matrix(c(1, 4.2, 2.2, 5.2), nrow = 2, ncol = 2, byrow = TRUE)
W4 <- matrix(c(9.05, 6.011, 8.3, 5.2), nrow = 2, ncol = 2, byrow = TRUE)

# 堆叠为3D数组
A_append <- array(c(A1, A2, A3, A4), c(2, 3, 4))
W_append <- array(c(W1, W2, W3, W4), c(2, 2, 4)) # 修正原笔误：替换A3、A4为W3、W4

目标函数为计算矩阵乘积的行列式：

crit.func <- function( A, W){
  return( det( t(A) %*% W %*% A))
}

预期输出结果：[2.095476e-12, 0, -7.067261e-12, 7.461713e-12]

尝试使用multiApply时出现维度匹配警告与错误结果，原调用代码：

library(multiApply)

A_append <- provideDimnames(A_append ,sep = "_", base = list('row','col','lev'))
W_append <- provideDimnames(W_append ,sep = "_", base = list('row','col','lev'))

D <- Apply(data = list(A_append, W_append), target_dims = c(1, 2, NULL), margins = 3, fun = crit.func)$output1

解决方案 #

方案一：修正multiApply的参数设置 #

错误根源在于target_dims参数未正确指定每个输入的二维维度，且维度名称匹配存在歧义。正确做法是明确指定每个输入的目标维度，并确保两个数组的循环维度（第三维度）名称一致：

library(multiApply)

# 手动设置清晰的维度名称，确保两个数组的第三维度（lev）完全对应
dimnames(A_append) <- list(row = 1:2, col = 1:3, lev = 1:4)
dimnames(W_append) <- list(row = 1:2, col = 1:2, lev = 1:4)

# 正确调用Apply：为每个输入指定要保留的二维维度，循环lev维度
D <- Apply(
  data = list(A_append, W_append),
  target_dims = list(c("row", "col"), c("row", "col")),
  margins = "lev",
  fun = crit.func
)$output1

# 验证结果
print(D)
# 输出：2.095476e-12 0.000000e+00 -7.067261e-12 7.461713e-12

方案二：原生R向量化方法（无需额外包） #

如果不想依赖multiApply，可使用原生apply函数或tidyverse的purrr包实现向量化循环，代码简洁且效率稳定：

方法1：原生apply循环

# 循环第三维度的每个索引，传入对应的二维矩阵
D <- apply(array(1:4, dim = 4), 1, function(i) {
  crit.func(A_append[,,i], W_append[,,i])
})

方法2：转为列表后用purrr映射

library(purrr)

# 将3D数组转为二维矩阵的列表
A_list <- lapply(1:4, function(i) A_append[,,i])
W_list <- lapply(1:4, function(i) W_append[,,i])

# 逐对调用crit.func并返回数值向量
D <- map2_dbl(A_list, W_list, crit.func)

两种方法均能得到预期结果，且原生方法无需额外安装包，适合大多数场景。

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