Самый быстрый способ умножения матрицы на вектор

У меня есть матрица mat и вектор v. Я хотел бы умножить первый столбец матрицы mat на первый элемент вектора v и умножить второй столбец матрицы mat на второй элемент вектора v. Я могу сделать это, как показано. Как я могу сделать это быстрее в R, так как мы получаем большую матрицу?

    mat = matrix(rnorm(1500000), ncol= 100)
    v= rnorm(100)
    > system.time( mat %*% diag(v))
      user  system elapsed 
      0.02    0.00    0.02

Ответы

Ответ 1

sweep кажется, работает немного быстрее на моей машине

sweep(mat, 2, v, FUN = "*")

Некоторые ориентиры:

> microbenchmark(mat %*% diag(v),sweep(mat, 2, v, FUN = "*"))

Unit: milliseconds
  expr       min        lq   median        uq      max neval
   %*% 214.66700 226.95551 231.2366 255.78493 349.1911   100
 sweep  42.42987  44.72254  62.9990  70.87403 127.2869   100

Ответ 2

Утилизация может ускоряться, но вы перерабатываете в столбцах, а не поперек, поэтому просто переносите и переносите назад.

t( t(mat) * v )

Это должно быть быстрее, чем sweep или %*%.

microbenchmark(mat %*% diag(v),sweep(mat, 2, v, FUN = "*"), t(t(mat)*v))
Unit: milliseconds
            expr       min        lq    median        uq      max neval
             %*% 150.47301 152.16306 153.17379 161.75416 281.3315   100
           sweep  35.94029  42.67210  45.53666  48.07468 168.3728   100
   t(t(mat) * v)  16.50813  23.41549  26.31602  29.44008 160.1651   100

Ответ 3

Немного поздно в игре, но кто-то сказал быстрее всего? Это может быть другим хорошим использованием Rcpp. Эта функция (называемая mmult) будет по умолчанию умножать каждый столбец матрицы на каждый последующий элемент вектора, но имеет возможность сделать это по столбцу, установив byrow = FALSE. Он также проверяет, что m и v имеют соответствующий размер с опцией byrow. Во всяком случае, он быстрый (примерно в 10-12 раз быстрее, чем лучший ответ на родной R)...

ИЗМЕНИТЬ

@chris предоставил этот отличный ответ к другому вопросу, который я задал, пытаясь заставить это работать с RcppArmadillo. Однако кажется, что функция Rcpp -одно, которую я размещал здесь, по-прежнему примерно в 8 раз быстрее, чем примерно, и примерно в 70 раз быстрее, чем метод OP. Нажмите ссылку для кода для функции @chris '- это красиво просто.

Я поставлю бенчмаркинг наверху.

require( microbenchmark )
m <- microbenchmark( mat %*% diag(v) , mmult( mat , v ) , sweep(mat, 2, v, FUN = "*") , chris( mat , v ) , t( t(mat) * v ) , times = 100L )
print( m , "relative" , order = "median" , digits = 3 )
Unit: relative
                        expr   min    lq median    uq   max neval
               mmult(mat, v)  1.00  1.00   1.00  1.00  1.00   100
               chris(mat, v) 10.74  9.31   8.15  7.27 10.44   100
               t(t(mat) * v)  9.65  8.75   8.30 15.33  9.52   100
 sweep(mat, 2, v, FUN = "*") 20.51 18.35  22.18 21.39 16.94   100
             mat %*% diag(v) 80.44 70.11  73.12 70.68 54.96   100

Просмотрите, как работает mmult и возвращает тот же результат, что и OP...

require( Rcpp )

#  Source code for our function
func <- 'NumericMatrix mmult( NumericMatrix m , NumericVector v , bool byrow = true ){
  if( byrow );
    if( ! m.nrow() == v.size() ) stop("Non-conformable arrays") ;
  if( ! byrow );
    if( ! m.ncol() == v.size() ) stop("Non-conformable arrays") ;

  NumericMatrix out(m) ;

  if( byrow ){
    for (int j = 0; j < m.ncol(); j++) {
      for (int i = 0; i < m.nrow(); i++) {
        out(i,j) = m(i,j) * v[j];
      }
    }
  }
  if( ! byrow ){
    for (int i = 0; i < m.nrow(); i++) {
      for (int j = 0; j < m.ncol(); j++) {
        out(i,j) = m(i,j) * v[i];
      }
    }
  }
  return out ;
}'

#  Make it available
cppFunction( func )

#  Use it
res1 <- mmult( m , v )

#  OP function
res2 <- mat %*% diag(v)

#  Same result?
identical( res1 , res2 ) # Yes!!
[1] TRUE