回火分数差分计算的矢量化答案

【问题标题】：Vectorization of tempered fractional differencing calculation回火分数差分计算的矢量化
【发布时间】：2016-03-31 05:14:04
【问题描述】：

我正在尝试加速这种缓和分数差分的近似。这控制了时间序列的长/准长内存。鉴于第一个 for 循环是迭代的，我不知道如何对其进行矢量化。此外，尝试向量化的输出与未更改的原始代码略有不同。感谢您的帮助。

原始代码

tempfracdiff= function (x,d,eta) {

n=length(x);x=x-mean(x);PI=numeric(n)
PI[1]=-d;TPI=numeric(n);ydiff=x

for (k in 2:n) {PI[k]=PI[k-1]*(k-1-d)/k}
for (j in 1:n) {TPI[j]=exp(-eta*j)*PI[j]}
for (i in 2:n) {ydiff[i]=x[i]+sum(TPI[1:(i-1)]*x[(i-1):1])}
return(ydiff)                  }

尝试向量化

tempfracdiffFL=function (x,d,eta) {

n=length(x);x=x-mean(x);PI=numeric(n)
PI[1]=-d;TPI=numeric(n);ydiff=x

for (k in 2:n) {PI[k]=PI[k-1]*(k-1-d)/k}
TPI[1:n]=exp(-eta*1:n)*PI[1:n]
ydiff[2:n]=x[2:n]+sum(TPI[1:(2:n-1)]*x[(2:n-1):1])
return(ydiff)          }

【问题讨论】：

标签： r for-loop vectorization

【解决方案1】：

对于 PI，您可以使用cumprod：

k <- 1:n
PI <- cumprod((k-1-d)/k)

TPI 可以不用索引来表示：

TPI <- exp(-eta*k)*PI

而ydiff 是x 加上x 和TPI 的卷积：

ydiff <- x+c(0,convolve(x,rev(TPI),type="o")[1:n-1])

所以，把它们放在一起：

mytempfracdiff = function (x,d,eta) {
  n <- length(x)
  x <- x-mean(x)
  k <- 1:n
  PI <- cumprod((k-1-d)/k)
  TPI <- exp(-eta*k)*PI
  x+c(0,convolve(x,rev(TPI),type="o")[1:n-1])
}

测试用例示例

set.seed(1)
x <- rnorm(100)
d <- 0.1
eta <- 0.5

all.equal(mytempfracdiff(x,d,eta), tempfracdiff(x,d,eta))
# [1] TRUE

library(microbenchmark)
microbenchmark(mytempfracdiff(x,d,eta), tempfracdiff(x,d,eta))

单位：微秒 expr min lq 平均中位数 uq mytempfracdiff(x, d, eta) 186.220 198.0025 211.9254 207.473 219.944 tempfracdiff(x, d, eta) 961.617 978.5710 1117.8803 1011.257 1061.816 最大内瓦尔 302.548 100 3556.270 100

【讨论】：

【解决方案2】：

对于 PI[k]，Reduce 很有帮助

n <- 5; d <- .3
fun <- function( a,b ) a * (b-1-d)/b
Reduce( fun, c(1,1:n), accumulate = T )[-1] # Eliminates PI[0]

[1] -0.30000000 -0.10500000 -0.05950000 -0.04016250 -0.02972025

【讨论】：