Прокачка суммы продукта по столбцу

Question

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

x<-matrix(seq(1:30), ncol=3) 
x 
     [,1] [,2] [,3] 
[1,] 1 11 21 
[2,] 2 12 22 
[3,] 3 13 23 
[4,] 4 14 24 
[5,] 5 15 25 
[6,] 6 16 26 
[7,] 7 17 27 
[8,] 8 18 28 
[9,] 9 19 29 
[10,] 10 20 30 
y<-matrix(rep(seq(1:3), 4), ncol=3)/10 
y 
    [,1] [,2] [,3] 
[1,] 0.1 0.2 0.3 
[2,] 0.2 0.3 0.1 
[3,] 0.3 0.1 0.2 
[4,] 0.1 0.2 0.3 

, так что результат будет выглядеть следующим образом:

1.8 9.9  20.3 
2.5 10.7 21.2 
3.2 11.5 22.1 
3.9 12.3 23 
4.6 13.1 23.9 
5.3 13.9 24.8 
6 14.7 25.7 

В примере вывода выше значение 10.7 рассчитывается как:

output[2, 2] = 12 * 0.2 + 13 * 0.3 + 14 * 0.1 + 15 * 0.2 

Кто-нибудь знает, как сделать это? Я играю с пакетом RcppRoll, но не могу получить правильный ответ. Чем быстрее решение, тем лучше, поскольку это часть оптимизации, которая займет много итераций.

Answer 1

Использование colSums:

t(
    sapply(1:(nrow(x) - nrow(y) + 1), function(i){ 
    colSums(x[i:((nrow(y)) + i - 1), ] * y) 
    }) 
) 

на основе больших примеров данных (при условии, в ответ ZheyuanLi в), microbenchmark:

Unit: milliseconds 
expr  min  lq  mean median  uq  max neval cld 
    zx 179.8928 186.8033 202.5204 192.3973 199.7500 299.5910 100 a 
    ZL 365.9814 368.3878 391.8303 370.0935 373.4502 489.5045 100 b 

Answer 2

Вы ищете convolution. В R функция convolve вычисляет свертку двух векторов через FFT (быстрое преобразование Фурье). Читайте ?convolve. Обратите внимание, что нам нужно type = "filter".

Например, свертка для x[,1] и y[,1] является:

convolve(x[,1], y[,1], type = "filter") 
# [1] 1.8 2.5 3.2 3.9 4.6 5.3 6.0 

Это просто обернуть вещи с sapply:

sapply(seq_len(ncol(x)), function (i) convolve(x[,i], y[,i], type = "filter")) 

#  [,1] [,2] [,3] 
#[1,] 1.8 9.9 20.3 
#[2,] 2.5 10.7 21.2 
#[3,] 3.2 11.5 22.1 
#[4,] 3.9 12.3 23.0 
#[5,] 4.6 13.1 23.9 
#[6,] 5.3 13.9 24.8 
#[7,] 6.0 14.7 25.7 

---

Я думаю, что в вашем контексте, ваша матрица x является тонкой, то есть она имеет гораздо больше строк, чем столбцы. Мой номер sapply находится вдоль столбца. Почему бы вам не пройти практический тест и не профилировать?

x <- matrix(rnorm(3000 * 100), 3000) ## `3000 * 100` matrix 
y <- matrix(rnorm(100 * 100), 100) ## `100 * 100` matrix 

Rprof("foo.out") 
sapply(seq_len(ncol(x)), function (i) convolve(x[,i], y[,i], type = "filter")) 
Rprof(NULL) 

summaryRprof("foo.out")$by.total 

       total.time total.pct self.time self.pct 
"sapply"    1.32 100.00  0.00  0.00 
"FUN"     1.30  98.48  0.02  1.52 
"lapply"    1.30  98.48  0.00  0.00 
"convolve"    1.28  96.97  0.08  6.06 
"fft"     1.12  84.85  1.12 84.85 
"rep.int"    0.04  3.03  0.04  3.03 
"array"    0.02  1.52  0.02  1.52 
"c"     0.02  1.52  0.02  1.52 
"Re"     0.02  1.52  0.02  1.52 
"simplify2array"  0.02  1.52  0.00  0.00 

96%+ времени тратится на convolve, таким образом, накладные расходы sapply ничтожна.

Answer 3

Это может быть сделано путем rollapply в одной строке, как это. Он использует весь объектный подход, то есть нет явного подписи.

library(zoo) 
rollapply(x, nrow(y), function(x) colSums(x*y), by.column = FALSE) 

дает:

 [,1] [,2] [,3] 
[1,] 1.8 9.9 20.3 
[2,] 2.5 10.7 21.2 
[3,] 3.2 11.5 22.1 
[4,] 3.9 12.3 23.0 
[5,] 4.6 13.1 23.9 
[6,] 5.3 13.9 24.8 
[7,] 6.0 14.7 25.7 

Примечание: Хотя не любой короче, используя magrittr это может попеременно быть записана в виде:

library(magrittr) 
library(zoo) 
x %>% rollapply(nrow(y), . %>% `*`(y) %>% colSums, by.column = FALSE)