У меня есть матрица m и вектор v. Я хотел бы умножить первый столбец матрицы m на первый элемент вектора v и умножить второй столбец матрицы m на второй элемент вектора v и так далее. Я могу сделать это со следующим кодом, но я ищу способ, который не требует двух транспозиционных вызовов. Как я могу сделать это быстрее в R?
m <- matrix(rnorm(120000), ncol=6)
v <- c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5)
system.time(t(t(m) * v))
# user system elapsed
# 0.02 0.00 0.02
Используйте некоторую линейную алгебру и выполняйте матричное умножение, которое довольно быстро выполняется в R.
eg
m %*% diag(v)
некоторый бенчмаркинг
m = matrix(rnorm(1200000), ncol=6)
v=c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5)
library(microbenchmark)
microbenchmark(m %*% diag(v), t(t(m) * v))
## Unit: milliseconds
## expr min lq median uq max neval
## m %*% diag(v) 16.57174 16.78104 16.86427 23.13121 109.9006 100
## t(t(m) * v) 26.21470 26.59049 32.40829 35.38097 122.9351 100
Если у вас больше столбцов, ваше решение t (t (m) * v) превосходит решение по умножению матриц с большим отрывом. Тем не менее, есть более быстрое решение, но оно требует больших затрат на использование памяти. Вы создаете матрицу размером m с помощью rep() и умножаете поэлементно. Вот сравнение, модифицирующее пример mnel:
m = matrix(rnorm(1200000), ncol=600)
v = rep(c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5), length = ncol(m))
library(microbenchmark)
microbenchmark(t(t(m) * v),
m %*% diag(v),
m * rep(v, rep.int(nrow(m),length(v))),
m * rep(v, rep(nrow(m),length(v))),
m * rep(v, each = nrow(m)))
# Unit: milliseconds
# expr min lq mean median uq max neval
# t(t(m) * v) 17.682257 18.807218 20.574513 19.239350 19.818331 62.63947 100
# m %*% diag(v) 415.573110 417.835574 421.226179 419.061019 420.601778 465.43276 100
# m * rep(v, rep.int(nrow(m), ncol(m))) 2.597411 2.794915 5.947318 3.276216 3.873842 48.95579 100
# m * rep(v, rep(nrow(m), ncol(m))) 2.601701 2.785839 3.707153 2.918994 3.855361 47.48697 100
# m * rep(v, each = nrow(m)) 21.766636 21.901935 23.791504 22.351227 23.049006 66.68491 100
Как видите, использование "each" в rep() жертвует скоростью для ясности. Разница между rep.int и rep кажется незначительной, обе реализации меняются местами при повторных запусках microbenchmark(). Имейте в виду, что ncol (m) == length (v).
Как отмечает @Arun, я не знаю, что вы будете бить свое решение с точки зрения эффективности времени. Что касается понятности кода, есть и другие варианты:
Один вариант:
> mapply("*",as.data.frame(m),v)
V1 V2 V3
[1,] 0.0 0.0 0.0
[2,] 1.5 0.0 0.0
[3,] 1.5 3.5 0.0
[4,] 1.5 3.5 4.5
И еще:
sapply(1:ncol(m),function(x) m[,x] * v[x] )
Ради полноты я добавил sweep в тест. Несмотря на несколько вводящие в заблуждение имена атрибутов, я думаю, что он может быть более читабельным, чем другие альтернативы, а также довольно быстрым:
n = 1000
M = matrix(rnorm(2 * n * n), nrow = n)
v = rnorm(2 * n)
microbenchmark::microbenchmark(
M * rep(v, rep.int(nrow(M), length(v))),
sweep(M, MARGIN = 2, STATS = v, FUN = '*'),
t(t(M) * v),
M * rep(v, each = nrow(M)),
M %*% diag(v)
)
Unit: milliseconds
expr min lq mean
M * rep(v, rep.int(nrow(M), length(v))) 5.259957 5.535376 9.994405
sweep(M, MARGIN = 2, STATS = v, FUN = '*') 16.083039 17.260790 22.724433
t(t(M) * v) 19.547392 20.748929 29.868819
M * rep(v, each = nrow(M)) 34.803229 37.088510 41.518962
M %*% diag(v) 1827.301864 1876.806506 2004.140725
median uq max neval
6.158703 7.606777 66.21271 100
20.479928 23.830074 85.24550 100
24.722213 29.222172 92.25538 100
39.920664 42.659752 106.70252 100
1986.152972 2096.172601 2432.88704 100
Как сделано bluegrue, достаточно простого представления, чтобы выполнить умножение по элементам.
Количество умножений и суммирования уменьшается на широкий диапазон, как если бы выполнялось простое матричное умножение с diag(), где для этого случая можно было бы избежать большого количества нулевых умножений.
m = matrix(rnorm(1200000), ncol=6)
v=c(1.5, 3.5, 4.5, 5.5, 6.5, 7.5)
v2 <- rep(v,each=dim(m)[1])
library(microbenchmark)
microbenchmark(m %*% diag(v), t(t(m) * v), m*v2)
Unit: milliseconds
expr min lq mean median uq max neval cld
m %*% diag(v) 11.269890 13.073995 16.424366 16.470435 17.700803 95.78635 100 b
t(t(m) * v) 9.794000 11.226271 14.018568 12.995839 15.010730 88.90111 100 b
m * v2 2.322188 2.559024 3.777874 3.011185 3.410848 67.26368 100 a
Если вы работаете с разреженными матрицами, вероятно, что поэлементные методы (т.е. расширение v до плотной матрицы) могут превысить ограничения памяти, поэтому я собираюсь исключить эти предлагаемые методы из этого анализа.
Пакет Matrix предоставляет метод Diagonal, который значительно улучшает метод диагонального умножения.
library(Matrix)
library(microbenchmark)
N_ROW = 5000
N_COL = 10000
M <- rsparsematrix(N_ROW, N_COL, density=0.1)
v <- rnorm(N_COL)
microbenchmark(
M %*% Diagonal(length(v), v),
t(t(M) * v),
sweep(M, MARGIN=2, STATS=v, FUN='*'))
# Unit: milliseconds
# expr min lq mean median uq max neval
# M %*% Diagonal(length(v), v) 36.46755 39.03379 47.75535 40.41116 43.30241 248.1036 100
# t(t(M) * v) 207.70560 223.35126 269.08112 250.25379 284.49382 511.0403 100
# sweep(M, MARGIN = 2, STATS = v, FUN = "*") 1682.45263 1869.87220 1941.54691 1924.80218 1999.62484 3104.8305 100