沿任意维重复的m维数组与(m-1)-dimensional阵列的比较

Question

我已经实现了一个基于多维数组的计算，它取代了一些循环代码。在这个过程中，我做了一些我认为可以做得更好的事情--但我不知道怎么做。

其中之一是将生成的3d数组与沿第三维空间重复的2d数组进行比较。

items12 = c(1,2,3,4,5,6)
items3 = c(1,2,3)

m2d = outer(items12, items12, "-")
m3d = outer(items3, m2d, "*")

经过一些操作之后，我想比较一下m2d和m3d，m2d沿着第三个dim重复。我知道有两种选择，两种方法都没有那么优雅，我很好奇有没有更好的方法。

实例化重复的3d数组。记忆沉重但速度快。

m2d.z.3d = outer(
  m2d, 
  rep(1, length(items3)), "*"
)

m3d - m2d.z.3d

循环播放。轻但慢。

apply(m3d, 3, function(x) {
    x - m2d
})

有什么建议吗？你会选哪一个？

更新示例澄清了任意索引要求。

items12 = c(1,2,3)
items3 = c(1,2)

m2d = outer(items12, items12, "-")
m3d = outer(m2d,items3, "*")

m3d - (m3d - items.3)

# items.3 wrapped along rows
, , 1

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    3
[3,]    1    2    3

, , 2

     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    2    3
[3,]    1    2    3

m3d.yx = aperm(m3d, c(2,1,3))
aperm(m3d.yx - (m3d.yx - c(items.3)), c(2,1,3)) 

#items.3 wrapped around columns
, , 1

     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3

, , 2

     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3

更新

在这种情况下的一些基准。

items.3 = rep(c(1,2,3), n)
items.2 = rep(c(1,2), n)

m2d = outer(items.3, items.3, "-")
m3d = outer(m2d, items.2, "*")

funRecycle = function() # items.3 wraps around the columns (index 1, then 2, then 3 etc.)
  m3d - (m3d - c(items.3)) 
funAperm = function() { # temporarily interchange index 1 and 2 to apply along desired index
  m3d.yx = aperm(m3d, c(2,1,3))
  aperm(m3d.yx - (m3d.yx - c(items.3)), c(2,1,3)) 
}
funOuter = function() { # assign the 3d matrix
  m2d.z.3d = outer(
    m2d, 
    rep(1, length(items.2)), "*"
  )
  m3d - m2d.z.3d
}
funArray = function() { # assign the 3d matrix with array
  m2d.z.3d = array(m2d, dim=c(dim(m2d)[1:2], length(items.2)))
  m3d - m2d.z.3d
}
funSweep <- function() sweep(m3d, c(1, 2), m2d, "-")

N=1

Unit: microseconds
         expr    min      lq     mean  median      uq    max neval   cld
 funRecycle()  1.110  1.3875  1.65388  1.6650  1.9420  2.775   100 a    
   funAperm() 17.200 19.1420 21.23113 20.2520 20.9455 69.077   100    d 
   funOuter() 14.426 15.8130 17.58316 17.2005 18.1710 35.232   100   c  
   funArray()  2.774  3.3300  3.95079  3.8840  4.1610 14.148   100  b   
   funSweep() 31.903 32.7360 34.84129 33.5680 34.4000 62.141   100     e

n=100

Unit: milliseconds
         expr       min        lq      mean    median        uq       max
 funRecycle()  28.51351  32.35671  37.13257  33.98931  39.94408  85.94085
   funAperm() 232.69297 276.07494 344.70083 352.40273 395.50492 569.54978
   funOuter()  35.25947  43.98674  53.06895  49.72790  55.93677  95.38608
   funArray()  96.78482 110.10501 119.68267 116.50378 120.70943 172.53973
   funSweep() 150.88675 168.90293 193.06270 178.11013 216.79349 291.23719

我对这样的结果感到惊讶，从某种程度上说，与数组()复制数组相比，将所有内容乘以1与外层相比要快得多。(总的来说，n outer()看起来可能比循环方法更快)。

如果我们必须对不同的索引(funAperm)进行比较，那么在所有情况下，使用外层构建数组的速度都要快得多。

除了aperm之外，还有什么建议可以在任意的索引中进行这种比较？

Answer 1

假设您是指(我假设这是因为否则m3d - m2d.z.3d不工作)：

m3d = outer(m2d, items3, "*") # note how I switched the arguments

然后这就起作用了：

m3d - c(m2d)

为了证明这一点：

all.equal(m3d - c(m2d), m3d - m2d.z.3d)
# [1] TRUE

在这里，我们只是利用向量循环，因为我们想重复沿最后一个维度。我们需要做c()来消除维度，否则R会抱怨数组是不可变形的(尽管它们实际上在我们这里想要的特定意义上)。

基于对R源代码(src/main/arithmetic.c:real_binary())的敷衍了事的审查，看来向量回收不会重复回收的向量，因此这应该是快速和内存高效的。

如果我们想在任意维上这样做，我们必须用aperm对数组进行所有维的重新洗牌，以使相关维度保持最后的位置，然后将结果重新洗牌回原来的维度顺序。这会增加一些开销。

至于选择哪种方法，如果没有耗尽内存，请使用快速方法(即避免循环，以利于完全向量化操作)。

另外，使用items12 <- seq(100)和items3 <- seq(50)的一些基准测试

funOuter <- function() {
  m2d.z.3d = outer(
    m2d, 
    rep(1, length(items3)), "*"
  )
  m3d - m2d.z.3d
}
funRecycle <- function() m3d - c(m2d)
funLoop <- function() apply(m3d, 3, "-", m2d)    # this does not appear correct because `apply` doesn't reconstruct dimensions like `sapply`
funSweep <- function() sweep(m3d, c(1, 2), m2d)  # this is the same type of thing but works properly

library(microbenchmark)
microbenchmark(funOuter(), funRecycle(), funLoop(), funSweep())

生产：

Unit: milliseconds
         expr       min        lq      mean    median
   funOuter()  2.297287  2.673768  3.232277  2.835404
 funRecycle()  1.327101  1.485082  2.093252  1.599543
    funLoop() 22.579010 24.586667 27.211804 26.840069
   funSweep() 11.251656 12.012664 13.516147 13.736908       

并核对结果：

all.equal(funOuter(), funRecycle())
# [1] TRUE
all.equal(funOuter(), funSweep())
# [1] TRUE
all.equal(funOuter(), funLoop())
# Nope, not equal