将一个向量加到矩阵的所有行

Question

我正在最大化一个似然函数并尝试减少循环次数。 我想将向量（要估计的参数）添加到矩阵（数据）的所有行中。向量的长度等于矩阵的列数。 a+b会产生错误的结果，因为R的循环利用规则是按列而不是按行。

a<-c(1,2,0,0,0)  # parameters to be optimized
b<-matrix(1,ncol=5,nrow=6) # data
t(a+t(b)) # my code would work, anything more intuitive?

期望的输出结果

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

错误的输出

a+b
    [,1] [,2] [,3] [,4] [,5]
[1,]    2    3    1    1    1
[2,]    3    1    1    1    2
[3,]    1    1    1    2    3
[4,]    1    1    2    3    1
[5,]    1    2    3    1    1
[6,]    2    3    1    1    1

- MLE

---

t(apply(b, 1 , function(x) x+a)) - Ronak Shah

Answer 1

我们可以使用 col 来复制 'a' 元素。

b + a[col(b)]
#     [,1] [,2] [,3] [,4] [,5]
#[1,]    2    3    1    1    1
#[2,]    2    3    1    1    1
#[3,]    2    3    1    1    1
#[4,]    2    3    1    1    1
#[5,]    2    3    1    1    1
#[6,]    2    3    1    1    1

或者更快的选择是使用rep

b + rep(a, each = nrow(b))

---

或者使用 sweep

sweep(b, 2, a, "+")

基准测试

set.seed(24)
b <- matrix(sample(0:9, 8000*5000, replace=TRUE), ncol=5000)
a <- sample(0:3, 5000, replace=TRUE)
system.time(b + a[col(b)])
#  user  system elapsed 
#  1.08    0.06    1.14 
system.time(b + rep(a, each = nrow(b)))
#   user  system elapsed 
#   0.83    0.03    0.86 

system.time(t(a+t(b)))
#   user  system elapsed 
#   1.14    0.03    1.17 

system.time(sweep(b, 2, a, "+"))
#  user  system elapsed 
#  0.62    0.06    0.69 

- akrun

---

我尝试过使用b+a[col(b)], 当数据较大时，运行速度似乎很慢。我将仔细检查我的代码并获得反馈。 - MLE

@Phdaml 在这种情况下，你可以尝试使用 rep，即 b + rep(a, each = nrow(b)) 应该更快。我更新了帖子。 - akrun

@Phdaml 我更新了一些使用 8000*5000 的基准测试。看起来 sweep 和 rep 更快。 - akrun

Answer 2

这些使用outer()或collapse::TRA()的解决方案比使用rep或sweep明显更快。

- Sebastian

Answer 3

另一个基准测试：

es=2:7

r=sapply(es,function(e){
  m=matrix(rnorm(10*10^e),ncol=10)
  v=rnorm(10)
  b=microbenchmark(times=10,
    t(t(m)+v),
    v[col(m)]+m,
    m+rep(v,each=nrow(m)),
    sweep(m,2,v,"+"),
    m+outer(rep(1,nrow(m)),v),
    collapse::TRA(m,v,"+"),
    Rfast::eachrow(m,v,"+")
  )
  a=aggregate(b$time,list(b$expr),median)
  setNames(a[,2],gsub(" ","",a[,1]))/1e6
})

r2=apply(r,2,function(x)formatC(x,max(0,2-ceiling(log10(min(x,na.rm=T)))),format="f"))
r3=apply(rbind(paste0("1e",es),r2),2,function(x)formatC(x,max(nchar(x)),format="s"))
writeLines(apply(cbind(r3,c("",rownames(r))),1,paste,collapse=" "))

输出（中位数时间，以毫秒为单位，1e7表示1千万行）：

   1e2   1e3  1e4  1e5 1e6  1e7 
0.0102 0.116 1.10 13.4  96 1478 t(t(m)+v)
0.0053 0.105 0.98  9.3  65 1225 v[col(m)]+m
0.0182 0.197 1.93 19.7 173 2044 m+rep(v,each=nrow(m))
0.0397 0.151 1.16 12.9  93 1431 sweep(m,2,v,"+")
0.0088 0.053 0.46  3.1  41  610 m+outer(rep(1,nrow(m)),v)
0.0049 0.037 0.32  3.1  16  402 collapse::TRA(m,v,"+")
0.0043 0.036 0.31  3.3  13  382 Rfast::eachrow(m,v,"+")

- nisetama