r - 动态高阶逐行乘积

Question

我对以下问题有疑问:

假设我有一些矩阵

A1 <- matrix(runif(rowsA1*T), rowsA1, T)
…
AD <- matrix(runif(rowsAD*T), rowsAD, T)

矩阵的数量是可变的(但肯定不会太大)。 有没有一种方法可以更有效地执行以下操作(但在允许可变数量矩阵的设置中):

f1 <- function(A1, A2, ..., AD) {
  for(i in 1:nrow(A1)) {
    for(j in 1:nrow(A2)) {
      ...
        for(d in 1:nrow(AD)) {
          ret[i,j,...,d] <- \sum_{t=1}^T (A1[i,t]*A2[j,t]*...*AD[d,t])
        }
      ...
    }
  }
ret
}

非常感谢您的帮助!

罗曼

-------------------------------- 编辑示例-------- ------------------------

 A1 <- |a b c|  A2 <- |j k l|  A3 <- |s t u|
       |d e f|        |m n o|        |v w x|
       |g h i|        |p q r|        |y z ä|

我想获得以下信息:

ret[1,1,1] <- a*j*s + b*k*t + c*l*u
ret[2,1,3] <- d*j*y + e*k*z + f*l*ä

希望这能让我的观点更清楚。

-------------------------------- 2013 年 11 月 26 日编辑 ------ --------------------------

嗨@flodel。我尝试实现您的代码，但一旦矩阵超过三个，似乎就会出现问题。

假设，我有以下矩阵

A1 <- matrix(runif(4*3), nrow = 4, ncol = 3)
A2 <- matrix(runif(3*3), nrow = 3, ncol = 3)
A3 <- matrix(runif(2*3), nrow = 2, ncol = 3)
A4 <- matrix(runif(1*3), nrow = 1, ncol = 3)

并将它们插入到您的代码中

output.f1 <- f1(A1,A2,A3,A4)

提供正确的维数

dim(output)
# [1] 4 3 2 1

但是输出全是NA

output.f1
# , , 1, 1

           # [,1] [,2] [,3]
# [1,] 0.13534704   NA   NA
# [2,] 0.07360135   NA   NA
# [3,] 0.07360135   NA   NA
# [4,] 0.07360135   NA   NA

# , , 2, 1

     # [,1] [,2] [,3]
# [1,]   NA   NA   NA
# [2,]   NA   NA   NA
# [3,]   NA   NA   NA
# [4,]   NA   NA   NA

感谢您的帮助...

最好的， 罗曼

Answer 1

试一试。对于大的 apply 循环，对于大矩阵它可能会很慢，但它可以完成对任意数量的矩阵通用的工作，而不需要相同的行数:

f1 <- function(...) {
  args    <- list(...)
  nrows   <- sapply(args, nrow)
  idx     <- do.call(expand.grid, lapply(nrows, seq.int))
  get.row <- function(i, mat) mat[i, ]
  get.val <- function(i.vec) sum(Reduce(`*`, Map(get.row, i.vec, args)))
  idx$val <- apply(idx, 1, get.val)
  ret     <- array(NA, dim = nrows)
  ret[as.matrix(idx[, seq_along(args)])] <- idx$val
  ret
}

示例用法:

A1 <- matrix(1:12, nrow = 4, ncol = 3)
A2 <- matrix(1:9,  nrow = 3, ncol = 3)
A3 <- matrix(1:6,  nrow = 2, ncol = 3)

out <- f1(A1, A2, A3)

检查:

identical(out[3, 2, 1],
          sum(A1[3, ] * A2[2, ] * A3[1, ]))
# [1] TRUE