有效地在数组上“应用”并保留结构

我有一个矩阵数组。

dims <- c(10000,5,5)
mat_array <- array(rnorm(prod(dims)), dims)

我想solve对每个矩阵执行基于矩阵的操作（例如通过函数求逆），但保留数组的完整结构。

到目前为止，我提出了 3 个选项：

选项 1：一个循环，它完全符合我的要求，但笨重且效率低下。

mat_inv <- array(NA, dims)
for(i in 1:dims[1]) mat_inv[i,,] <- solve(mat_array[i,,])

选项 2：该apply函数更快、更清晰，但将每个矩阵压缩为一个向量。

mat_inv <- apply(mat_array, 1, solve)
dim(mat_inv)
[1]    25 10000

我知道我可以设置输出维度以匹配输入的维度，但是我很担心这样做会弄乱索引，特别是如果我必须在不相邻的维度上应用（例如，如果我想在维度 2 上反转））。

选项 3：包中的aaply函数plyr，它完全符合我的要求，但比其他函数慢得多（4-5 倍）。

mat_inv <- plyr::aaply(mat_array, 1, solve)

是否有任何选项可以将的速度base::apply与的多功能性相结合plyr::aaply？