Matlab bsxfun(@times,...,...) 等价于 R

Matlab bsxfun(@times,...,...) equivalent in R

R 是否得到了 Matlab bsxfun(@times,a,b) 的等价物?假设有人想在矩阵 a,b:

上执行 element wise multiplication

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 的结果,因为我希望它与 Matlabbsxfun(@times,a,b)

编辑:

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

最佳

用列主矩阵来做,因为那是 R-convention:

> 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-矩阵与您的 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

"bsxfun" 和 "sweep" 这两种方法从执行速度的角度来看都不是最优的。 "bsxfun"的性能其实很惨淡,对非玩具类的应用好像没什么用。

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

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

sweep() 在矩阵大小增加时变得更具竞争力,但即使对于中等大小的矩阵 (1e6-by-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

更新:确实存在一些错误,扫描会在 (1,10,2) 等形状上阻塞。因此,我根据 sweep 的编写方式编写了另一个实现。由于重铸,速度以某种方式受到影响,但仍快 4 倍。该实现存放在 github 回购 bsxfun.R(基准来自@g g)

我想缩写要扫描的边距,但 pracma::bsxfun() 没有使用 sweep()。所以我写了以下内容(参见 github link),希望没有错误。

标签:R、pracma、数组