按列滚动和乘积

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

您正在寻找 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 2

使用 colSums:

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

基于更大的示例数据（在 ZheyuanLi 的回答中提供），微基准测试：

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 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)

按列滚动和乘积

Rolling sum product by column

r

matrix

rcpp