Matlab bsxfun(@times，.，.)等效R

Question

R是否与Matlab bsxfun(@times,a,b)相当？假设一个人想在矩阵element wise multiplication上执行a,b

Matlab:

a=[1 0 3 -4];
b=[0 1 5 7; 2 9 -3 4];

bsxfun(@times,a,b) = [0 0 15 -28; 2 0 -9 -16]

R:

a<-c(1,0,3,-4)
b<-matrix(c(0,2,1,9,5,-3,7,4),nrow = 2,ncol = 4)

a * b = matrix(c(0,0,3,-36,5,0,21,-16),nrow = 2,ncol = 4)

对于R获得上述a*b结果的方式有什么想法，因为我原以为它与Matlab bsxfun(@times,a,b)是完全相同的

编辑：

bsxfun("*",repmat(a,2,1),b) # using R {pracma}

最好的

Answer 1

使用列主要矩阵来实现它，因为这是R-约定：

> b<-matrix(c(0,2,1,9,5,-3,7,4),nrow = 4,ncol = 2)
> a*b
     [,1] [,2]
[1,]    0    5
[2,]    0    0
[3,]    3   21
[4,]  -36  -16

如果您使用最初的b结构，那么当您尝试使用sweep时，会有一些不愉快的惊喜。

> b2<-matrix(c(0,2,1,9,5,-3,7,4),nrow = 2,ncol = 4)
> sweep(b2, 2, a, '*')
     [,1] [,2] [,3] [,4]
[1,]    0    0   15  -28
[2,]    2    0   -9  -16

由于matrix函数使用列主要的位置填充，而且您没有在其调用中指定byrow=TRUE，所以b-matrix与Matlab矩阵不同。

> b3<-matrix(c(0,2,1,9,5,-3,7,4),nrow = 2,ncol = 4, byrow=TRUE)
> sweep(b3, 2, a, '*')
     [,1] [,2] [,3] [,4]
[1,]    0    0    3  -36
[2,]    5    0   21  -16

Answer 2

从执行速度的角度来看，这两种方法都不是最佳的。"bsxfun“的表现实际上是如此令人沮丧，以至于对于非玩具应用程序来说似乎是毫无用处的。

在我的基准测试中，最快的方法是

matrix(a, ncol = n_col, nrow = nsample, byrow = TRUE) * b

当矩阵大小增加时，sweep()变得更有竞争力，但即使对于中等大的矩阵(1e6乘4)，也要花费双倍的时间。

在下面找到完整的基准

# test of recycling efficiency

rm(list=ls())

library(microbenchmark)
library(pracma)

a = c(1,0,3,-4)
n_col = 4

# make example more realistic by expanding number of rows of b
nsample = 1e3
b = repmat(matrix(c(0,2,1,9,5,-3,7,4),nrow = 2,ncol = n_col), nsample / 2, 1)

print(microbenchmark(
  erg_1 = matrix(a, ncol = n_col,  nrow = nsample, byrow = TRUE) * b,
  erg_2 = matrix(rep.int(a, nsample), nrow = nsample, ncol = n_col, byrow = TRUE) * b,
  erg_3 = matrix(a * c(t(b)), nrow = nsample, ncol = n_col, byrow = TRUE),
  erg_4 = sweep(b, 2, a, '*'),
  erg_5 = bsxfun('*', repmat(a, nsample, 1), b) 
  ))


#same as above but now larger matrices
nsample = 1e6
b = repmat(matrix(c(0,2,1,9,5,-3,7,4),nrow = 2,ncol = n_col), nsample / 2, 1)

print(microbenchmark(
  erg_1 = matrix(a, ncol = n_col,  nrow = nsample, byrow = TRUE) * b,
  erg_2 = matrix(rep.int(a, nsample), nrow = nsample, ncol = n_col, byrow = TRUE) * b,
  erg_3 = matrix(a + c(t(b)), nrow = nsample, ncol = n_col, byrow = TRUE),
  erg_4 = sweep(b, 2, a, '*')
  #erg_5 = bsxfun('*', repmat(a, nsample, 1), b) #bsxfun is non-competitive
))

>Unit: microseconds
  expr      min        lq       mean    median        uq      max neval
 erg_1    9.057   10.1135   11.93394   11.0195   12.6790   36.226   100
 erg_2   14.189   15.3970   18.75324   16.9060   19.9250   41.358   100
 erg_3   26.263   28.8295   35.04538   30.9430   34.8675   86.941   100
 erg_4   40.452   44.0750   56.88289   51.4705   66.4130  109.279   100
 erg_5 2694.827 2968.4755 3243.76025 3208.9185 3417.5125 5575.306   100
>Unit: milliseconds
  expr      min       lq     mean   median       uq      max neval
 erg_1 10.85538 11.30668 20.58625 12.93408 13.28290 69.17918   100
 erg_2 16.07206 18.00058 29.17394 18.24751 20.09845 75.30993   100
 erg_3 22.41231 24.58957 30.83620 24.99544 26.49047 79.71910   100
 erg_4 20.74838 21.53673 29.52071 22.88867 23.30420 81.07150   100

Answer 3

更新:确实有一些错误，像(1,10,2)这样的形状会让扫描窒息。因此，我根据扫描的编写方式编写了另一个实现。速度在某种程度上由于重铸而受到损害，但速度仍然快了4倍.实现保存在github bsxfun.R (从@g获取基准标记)上。

我想缩略一下要扫描的边距，但是pracma::bsxfun（）没有使用sweep()。所以我写了以下文章(参见github链接)，希望没有bug。

标签: R，pracma，数组