使用循环顺序更新data.table列

我已经使用data.table实现了here所描述的一个简单的动态编程示例，希望它能像矢量化代码一样快。使用循环顺序更新data.table列

library(data.table) 
B=100; M=50; alpha=0.5; beta=0.9; 
n = B + M + 1 
m = M + 1 
u <- function(c)c^alpha 
dt <- data.table(s = 0:(B+M))[, .(a = 0:min(s, M)), s] # State Space and corresponging Action Space 
dt[, u := (s-a)^alpha,]        # rewards r(s, a) 
dt <- dt[, .(s_next = a:(a+B), u = u), .(s, a)]  # all possible (s') for each (s, a) 
dt[, p := 1/(B+1), s]         # transition probs 

#   s a s_next u  p 
#  1: 0 0  0 0 0.009901 
#  2: 0 0  1 0 0.009901 
#  3: 0 0  2 0 0.009901 
#  4: 0 0  3 0 0.009901 
#  5: 0 0  4 0 0.009901 
# ---       
#649022: 150 50 146 10 0.009901 
#649023: 150 50 147 10 0.009901 
#649024: 150 50 148 10 0.009901 
#649025: 150 50 149 10 0.009901 
#649026: 150 50 150 10 0.009901

给一点内容，我的问题：在s和a的s（s_next）未来值的条件被实现为a:(a+10)一个，每个概率p=1/(B + 1)。 u列给出了每种组合(s, a)的u(s, a)。

给定初始值V（总是n by 1向量）为每个唯一的状态s，V根据V[s]=max(u(s, a)) + beta* sum(p*v(s_next))（贝尔曼方程）的更新。
最大化为a，因此，[, `:=`(v = max(v), i = s_next[which.max(v)]), by = .(s)]在下面的迭代中。

实际上有非常高效的vectorized solution。我认为data.table解决方案的性能与矢量化方法相当。

我知道主要罪魁祸首是dt[, v := V[s_next + 1]]。唉，我不知道如何解决它。

# Iteration starts here 
system.time({ 
    V <- rep(0, n) # initial guess for Value function 
    i <- 1 
    tol <- 1 
    while(tol > 0.0001){ 
    dt[, v := V[s_next + 1]] 
    dt[, v := u + beta * sum(p*v), by = .(s, a) 
     ][, `:=`(v = max(v), i = s_next[which.max(v)]), by = .(s)] # Iteration 
    dt1 <- dt[, .(v[1L], i[1L]), by = s] 
    Vnew <- dt1$V1   
    sig <- dt1$V2 
    tol <- max(abs(V - Vnew)) 
    V <- Vnew 
    i <- i + 1 
    }  
}) 
# user system elapsed 
# 5.81 0.40 6.25

令我沮丧的是，data.table解决方案甚至比以下非矢量化解决方案还要慢。作为一个马虎data.table用户，我必须错过data.table功能。有没有办法改善事情，或者，data.table不适合这种计算？

S <- 0:(n-1) # StateSpace 
VFI <- function(V){ 
    out <- rep(0, length(V)) 
    for(s in S){ 
    x <- -Inf 
    for(a in 0:min(s, M)){ 
     s_next <- a:(a+B)  # (s') 
     x <- max(x, u(s-a) + beta * sum(V[s_next + 1]/(B+1))) 
    } 
    out[s+1] <- x 
    } 
    out 
} 
system.time({ 
V <- rep(0, n) # initial guess for Value function 
i <- 1 
tol <- 1 
while(tol > 0.0001){ 
    Vnew <- VFI(V)   
    tol <- max(abs(V - Vnew)) 
    V <- Vnew 
    i <- i + 1 
}  
}) 
# user system elapsed 
# 3.81 0.00 3.81

来源

2017-01-01 Khashaa

请参阅https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example。有人可能会抽出时间来解决这个问题，但减少到最简单的问题演示（在你的情况下，使用data.table缓慢）将会得到更好的结果。 –

@JackWasey你真的有一些神经。你真的认为这个链接是需要的吗？我认为Khashaa知道r/data.table没有比你更糟，并知道如何创建MWE。如果你不能帮助，你可以继续前进 - 不需要自负的评论。 –

如果你的问题的主要目标是如何提高data.table方法的性能，那么也许别人可以提供帮助。但是，如果你只是在寻找改善这些动态模型性能的方法，那么我个人总是使用RCpp来处理这种事情。向量化动态模型通常很棘手，而且通常不可能。如果需要速度，RCpp通常是最好的选择。 – dww

这是我会怎么做这个...

DT = CJ(s = seq_len(n)-1L, a = seq_len(m)-1L, s_next = seq_len(n)-1L) 
DT[ , p := 0] 
#p is 0 unless this is true 
DT[between(s_next, a, a + B), p := 1/(B+1)] 
#may as well subset to eliminate irrelevant states 
DT = DT[p>0 & s>=a] 
DT[ , util := u(s - a)] 

#don't technically need by, but just to be careful 
DT[ , V0 := rep(0, n), by = .(a, s_next)] 

while(TRUE) { 
    #for each s, maximize given past value; 
    # within each s, have to sum over s_nexts, 
    # to do so, sum by a 
    DT[ , V1 := max(.SD[ , util[1L] + beta*sum(V0*p), by = a], 
       na.rm = TRUE), by = s] 
    if (DT[ , max(abs(V0 - V1))] < 1e-4) break 
    DT[ , V0 := V1] 
}

在我的机器大约需要15秒（所以不太好）......但也许这会给你一些想法。例如，此data.table太大，因为最终只有n的唯一值V。

来源

2017-01-03 22:27:18 MichaelChirico

使用循环顺序更新data.table列

回答

相关问题