ループを使ってdata.table列を順次更新する

Question

私は、ベクトル化されたコードと同じくらい速くなることを願って、hereと記述された単純な動的プログラミングの例をdata.tableを使って実装しました。 s（s_next）のsとa、将来の値を条件a:(a+10)の一つ、確率p=1/(B + 1)と各として実現されています

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 

は私の質問にはほとんどコンテンツを与えるために。 u列は、各組み合わせについて u(s, a)を与える。 (s, a)。

• 各固有状態sの初期値V（常にn by 1ベクトル）が与えられると、V[s]=max(u(s, a)) + beta* sum(p*v(s_next))（ベルマン式）に従ってVアップデート。
• 最大化は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 

Answer 1

が、これは（あまりよくない）、約15秒かかり、私のマシンで

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] 
} 

...私はこれを行うだろうかだ...多分これはあなたにいくつかのアイデアを与えるだろう。たとえば、は実際にはVのnという一意の値しかないため、実際には大きすぎます。