将数组中矩阵的每个切片除以其自己的向量？答案

【问题标题】：Divide every slice of a matrix in an array by its own vector?将数组中矩阵的每个切片除以其自己的向量？
【发布时间】：2015-03-10 19:40:53
【问题描述】：

假设我有两个数组（或张量，如果需要张量包）

dim(Xbeta)
  products      draws Households
        13      20        10

dim(denom)
        1       20        10


set.seed(1)
Xbeta=array(rnorm(13*20*10,0,1),dim=c(13,20,10))
denom=array(rnorm(1*20*10,0,1),dim=c(1,20,10))

没有循环，我想做以下事情：

   for(i in 1:10){

    Xbeta[,,i]=t(t(Xbeta[,,i]) / denom[,,i])

       }

我想将Xbeta[,,i] 切片中的每一列除以denom[,,i]. 中的每个对应数字

例如...Xbeta[,1,i]/denom[,1,i]...等

【问题讨论】：

我也一直想知道如何矢量化这样的操作。一个干净、无循环的解决方案将非常受欢迎。

标签： arrays r matrix linear-algebra

【解决方案1】：

您可以通过 (1) 3 维转置分子数组和 (2) 将分母数组展平为向量来避免循环和复制，这样除法运算自然会在整个转置后循环不完整的分母向量分子数组，使数据按照您想要的方式排列。然后，您必须对结果进行 3 维“反转置”以恢复原始转置。

aperm(aperm(Xbeta,c(2,3,1))/c(denom),c(3,1,2));

第一次调用aperm() 将列转置为行，将z 切片转置为列，将行转置为z 切片。 c() 对 denom 的调用将分母数组展平为向量，因为在循环时，我们不关心维度。对aperm() 的最后调用会反转原来的转置。

要详细了解此解决方案的逻辑，您的输入基本上是分子数组的每个 z 切片的除数向量，并且您希望将每个除数应用于相应 z 的每一行- 切片和列。这意味着除数向量必须首先应用于列，然后，随着每个分母 z-slice 用尽，应用于分子 z-slice。在分子数组的一整行（覆盖该行中的所有 z 切片）用尽后，分母向量的整体已经用尽，导致它循环回到分子数组的下一行的开头。

见https://stat.ethz.ch/R-manual/R-devel/library/base/html/aperm.html。

关于性能的粗略想法：

r> set.seed(1);
r> Xbeta <- array(rnorm(13*20*10,0,1), dim=c(13,20,10) );
r> denom <- array(rnorm(1*20*10,0,1), dim=c(1,20,10) );
r> robert <- function() { result <- array(NA, dim=c(13,20,10) ); for (i in 1:10) { result[,,i] <- t(t(Xbeta[,,i]) / denom[,,i]); }; };
r> andre <- function() { denom_myVersion <- array(rep(c(denom), each=13 ), c(13,20,10) ); result <- Xbeta / denom_myVersion; };
r> bgoldst <- function() { result <- aperm(aperm(Xbeta,c(2,3,1))/c(denom),c(3,1,2)); };
r> N <- 99999;
r> system.time({ replicate(N, robert() ); });
   user  system elapsed
 25.421   0.000  25.440
r> system.time({ replicate(N, andre() ); });
   user  system elapsed
 12.578   0.594  13.283
r> system.time({ replicate(N, bgoldst() ); });
   user  system elapsed
  8.484   0.594   9.142

此外，作为一般建议，使用最少的样本输入来呈现这类问题（对于提问者和回答者）是有帮助的，例如：

r> n <- array(1:12,dim=c(2,3,2)); n;
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

r> d <- array(1:6,dim=c(1,3,2)); d;
, , 1

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

, , 2

     [,1] [,2] [,3]
[1,]    4    5    6

r> aperm(aperm(n,c(2,3,1))/c(d),c(3,1,2));
, , 1

     [,1] [,2]     [,3]
[1,]    1  1.5 1.666667
[2,]    2  2.0 2.000000

, , 2

     [,1] [,2]     [,3]
[1,] 1.75  1.8 1.833333
[2,] 2.00  2.0 2.000000

【讨论】：

【解决方案2】：

# Is this what you're looking for?
Xbeta <- array(rnorm(13*20*10,0,1),dim=c(13,20,10))
denom <- array(rnorm(1*20*10,0,1),dim=c(1,20,10))
div.list <- sapply(1:10, FUN = function(x) t(Xbeta[,,x]) / denom[,,x], simplify = FALSE) 
result <-   array(do.call(c, div.list), dim = dim(Xbeta)[c(2,1,3)])

【讨论】：

【解决方案3】：

我不确定您为什么为 denon 选择 3 维数组。无论如何，这可以通过密切关注这些数字在内存中的存储方式来完成。在数组中，第一个维度“移动最快”。通过将 denom 值“每个”复制 13 次，您可以创建一个与分子具有完全相同维度的数组。

所以，让我们测试一下：让我们保存随机值，以便我们可以将它们用于两种方法：

set.seed(1)
Num_2600 <- rnorm(13*20*10,0,1)
Denom_200 <- rnorm(20*10,0,1)
Xbeta=array(Num_2600,dim=c(13,20,10))
denom=array(Denom_200,dim=c(1,20,10))
Your_result <- array(NA, dim=c(13,20,10))

您的代码给出：

for(i in 1:10){
Your_result[,,i] <- t(t(Xbeta[,,i]) / denom[,,i])
   }

我的代码：

denom_myVersion <- array(rep( Denom_200 , each=13), c(13,20,10))

>  all(Your_result == Xbeta / denom_myVersion)
[1] TRUE
>

所以我们得到了相同的结果。困难的部分是如何决定如何复制，以便数字落在正确的位置。注意：

denom_myVersion <- array(rep( Denom_200 , times=13), c(13,20,10))

>  all(Your_result == Xbeta / denom_myVersion)
[1] FALSE
>

在rep 中使用'each' 作为参数，每个元素在转到下一个元素之前重复13 次。随着时间的推移，整个向量重复 13 次。比较：

> rep(1:3, each =3)
[1] 1 1 1 2 2 2 3 3 3
> rep(1:3, times=3)
[1] 1 2 3 1 2 3 1 2 3

【讨论】：