将它们组合成一个数组并使用apply：

A <- matrix(c(2,4,3,5), 2)
B <- matrix(c(6,8,7,9), 2)

X <- list(A, B)
Y <- do.call(cbind, X)
Y <- array(Y, dim=c(dim(X[[1]]), length(X)))

apply(Y, c(1, 2), mean, na.rm = TRUE)
#     [,1] [,2]
#[1,]    4    5
#[2,]    6    7

如果apply的效率不够高，可以将colMeans（提供NA处理）与aperm一起使用：

colMeans(aperm(Y, c(3, 1, 2)), na.rm = TRUE)
#     [,1] [,2]
#[1,]    4    5
#[2,]    6    7

• • 更新：**任何喜欢base-R的人都会想要Reduce（请在下面的邮件中投票）：

Reduce(`+`, x) / length(x)

• • 原件：**

我认为现在有一个更好的解决方案，多亏了purrr::reduce。

require(purrr)
x = list(volcano+10, volcano-15, volcano-5)
v = reduce(x, `+`) / length(x)
image(v)

如果希望结果是具有相同列名和行名称的矩阵，可以创建一个函数将矩阵列表传递到。还可以将应用的函数从mean更改为sd或任何其他函数。

reduceApplyListOfArrays<- function(x){
     y<-apply(array(unlist(x), c(dim(x[[1]]), dim(x[[2]]), length(x))), 
              c(1,2), mean)
     colnames(y)<-colnames(x[[1]])
     rownames(y)<-rownames(x[[1]])
  return(y)
}

简单的解决办法是

A <- matrix(c(2,4,3,5), 2)
B <- matrix(c(6,8,7,9), 2)
(A+B)/2

得到了与“罗兰解”相同的结果