两个 3D 数组和一个 2D 矩阵的乘积答案

【问题标题】：Product of two 3D array and a 2D matrix两个 3D 数组和一个 2D 矩阵的乘积
【发布时间】：2019-08-04 10:23:48
【问题描述】：

我正在尝试找到一种更有效的方法来在 R 中对以下矩阵进行编码：设 A 和 C 是两个维数为 (n, n, m) 的 3D 数组， B 是维数 (m, m) 的矩阵，则 M 是一个 (n, n) 矩阵，这样：

M_ij = SUM_kl A_ijk * B_kl * C_ijl

for (i in seq(n)) {
  for (j in seq(n)) {
    M[i, j] <- A[i,j,] %*% B %*% C[i,j,]
    }
  }

可以使用 TensorA 包使用 i 和 j 作为平行维度来编写此代码，但我更愿意使用基本 R 对象。

einstein.tensor(A %e% log(B), C, by = c("i", "j"))

谢谢！

【问题讨论】：

A[I,j,] 不是向量吗？您不能通过省略 I 或 j 然后进行矩阵乘法将其编写为单个循环吗？如果你想“放弃”循环，应该有一种方法 apply

标签： arrays r

【解决方案1】：

我不知道这是否会更快，但它会避免一级循环：

for (i in seq(n))
  M[i,] <- diag(A[i,,] %*% B %*% t(C[i,,]))

在这个例子中它给出了和你一样的答案：

n <- 2
m <- 3
A <- array(1:(n^2*m), c(n, n, m))
C <- A + 1
B <- matrix(1:(m^2), m, m)

M <- matrix(NA, n, n)
for (i in seq(n))
  M[i,] <- diag(A[i,,] %*% B %*% t(C[i,,]))
M
#      [,1] [,2]
# [1,] 1854 3216
# [2,] 2490 4032

编辑添加：基于https://stackoverflow.com/a/42569902/2554330，这是一个稍快的版本：

for (i in seq(n))
  M[i,] <- rowSums((A[i,,] %*% B) * C[i,,])

我对@987654325@ 和m <- 300 做了一些计时，最快的时间是 3.1 秒，而我原来的解决方案是 4.7 秒，而问题中的解决方案是 17.4 秒。

【讨论】：