R 中使用 bandSparse 和矩阵乘法的线性卷积

Linear convolution in R using bandSparse & matrix multiplication

在 R 中,我可以使用

计算向量 y 与卷积核 K 的线性卷积
linconv = function (y,K) {   K=K/sum(K)
                             out=convolve(y, K, type="open")
                             out=head(out,-(length(out)-length(y)) )
                             return(out)
}

例如如果我有一个向量 y

t = 1:1000
y = 10000*exp(sin(t*0.02))

和卷积核K

K = dexp(1:41, rate=1/5)

y与卷积核K的线性卷积将是linconv(y, K)

现在我想计算我的线性卷积,而不是使用矩阵乘法,方法是设置一个带状卷积矩阵,其中包含我的卷积核 K 的时移副本作为列。理想情况下,我想使用 Matrix::bandSparse 来设置这个带状卷积矩阵。有没有人碰巧知道 bandSparse 的正确语法来完成这个,并获得与使用上面的线性卷积函数相同的结果,但使用矩阵乘法代替?现在我做了

linconv2 = function (y,K) {  K=K/sum(K)
                             require(Hmisc)
                             lags = 0:(N-1)
                             X = as.matrix(data.frame(do.call(cbind, lapply(lags, function(lag) Lag(c(K, rep(0, length(y))), shift=lag)))))
                             X[is.na(X)] = 0
                             ntoclip = (nrow(X)-length(y))/2
                             X = X[head(1:nrow(X), -(nrow(X)-length(y))),]
                             out = as.vector(y %*% X)                                 
                             return(out)
}

但我想在上面设置矩阵 X 而不是使用 Matrix::bandSparse。有人知道这样做的正确语法吗?

终于找到正确的语法了:

linconv2 = function (y,K) {  K=K/sum(K)
                             require(Matrix)
                             X = as.matrix(bandSparse(length(y), 
                                              k = seq(-(length(K)-1),0,1), 
                                              diag = t(replicate(length(y), rev(K))), symm=FALSE))

                             out = X %*% as.matrix(y, ncol=1)
                             return(out)
}