Aktualisieren einer data.table Spalte sequentiell mit Schleife

Question

Ich habe ein einfaches dynamisches Programmierbeispiel beschrieben here, mit data.table, in der Hoffnung, dass es so schnell wie vektorisierter Code sein würde. Voraussetzung für s und a, zukünftige Werte von s (s_next) gilt als einer der a:(a+10), die jeweils mit einer Wahrscheinlichkeit von p=1/(B + 1) realisiert:

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 

Um einen wenig Inhalt auf meine Frage zu geben. u Spalte gibt die u(s, a) für jede Kombination (s, a).

• die Anfangswerte V (immer n by 1 Vektor) für jeden einzelnen Zustand s, Gegeben V Updates nach V[s]=max(u(s, a)) + beta* sum(p*v(s_next)) (Bellman Gleichung).
• Maximierung ist wrt a, daher [, `:=`(v = max(v), i = s_next[which.max(v)]), by = .(s)] in der Iteration unten.

Eigentlich gibt es sehr effizient vectorized solution. Ich dachte data.table Lösung wäre vergleichbar in der Leistung als vektorisierte Ansatz.

Ich weiß, dass der Hauptschuldige ist dt[, v := V[s_next + 1]]. Ach, ich habe keine Ahnung, wie ich es beheben kann.

# 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 

Zu meiner Bestürzung, die data.table Lösung ist sogar langsamer als die folgenden äußerst nicht vektorisiert Lösung. Als schlampiger data.table-user muss ich einige data.table Funktionen vermissen. Gibt es eine Möglichkeit, Dinge zu verbessern, oder data.table ist nicht für diese Art von Berechnungen geeignet?

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

Hier ist, wie ich dies tun würde ...

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

Auf meinem Rechner dies etwa 15 Sekunden in Anspruch nimmt (nicht so gut) ... aber vielleicht wird diese Ihnen ein paar Ideen. Zum Beispiel ist dies data.table viel zu groß, da es wirklich nur n eindeutige Werte von V letztlich gibt.