6.231 Fall 2015 Lecture 4 Notes

6.231 DYNAMIC PROGRAMMING LECTURE 4 LECTURE OUTLINE • Review of approximation in value space • Approximate VI and PI • Projected Bellman equations • Matrix form of the projected equation • Simulation-based implementation • LSTD and LSPE methods • Optimistic versions • Multistep projected Bellman equations • Bias-variance tradeoff 1 DISCOUNTED MDP • System: Controlled Markov chain with states i = 1, . . . , n and finite set of controls u ∈ U (i) • Transition probabilities: pij (u)                • Cost of a policy π = {µ0 , µ1 , . . .} starting at state i: # !N " % $ k Jπ (i) = lim E α g ik , µk (ik ), ik+1 | i = i0 N →∞ k=0 with α ∈ [0, 1) • Shorthand notation for DP mappings (T J)(i) = min u∈U (i) (Tµ J)(i) = n " j=1 n " % pij (u) g(i, u, j)+αJ(j) , i = 1, . . . , n, j=1 $ $ pij µ(i) %$ $ % % g i, µ(i), j +αJ(j) , i = 1, . . . , n 2 “SHORTHAND” THEORY – A SUMMARY • Bellman’s equation: J ∗ = T J ∗ , Jµ = Tµ Jµ or J ∗ (i) Jµ (i) = j =1 g(i, u, j) + αJ ∗ (j) % ∀i $ %$ $ % % pij µ(i) g i, µ(i), j + αJµ (j) , ∀i u∈U (i) pij (u) $ , = min n " n " j=1 • Optimality condition: µ: optimal <==> Tµ J ∗ = T J ∗ i.e., µ(i) ∈ arg min u∈U (i) n " j=1 % $ ∗ pij (u) g(i, u, j)+αJ (j) , ∀ i 3 THE TWO MAIN ALGORITHMS: VI AND PI • Value iteration: For any J ∈ ,n J ∗ (i) = lim (T k J)(i), ∀ i = 1, . . . , n k→∞ • Policy iteration: Given µk − Policy evaluation: Find Jµk by solving Jµk (i) = n " j=1 $ k pij µ (i) %$ $ k % % g i, µ (i), j +αJµk (j) , i = 1, . . . , n or Jµk = Tµk Jµk − Policy improvement: Let µk+1 be such that µ k+1 (i) ∈ arg min u∈U (i) n " $ % pij (u) g(i, u, j)+αJµk (j) , ∀ i j =1 or Tµk+1 Jµk = T Jµk • Policy evaluation is equivalent to solving an n × n linear system of equations • For large n, exact PI is out of the question (even though it terminates finitely) 4 APPROXIMATION IN VALUE SPACE • Approximate J ∗ or Jµ from a parametric class J˜(i, r), where i is the current state and r = (r1 , . . . , rm ) is a vector of “tunable” scalars weights. • By adjusting r we can change the “shape” of J˜ so that it is close to the true optimal J ∗ . • Any r ∈ ,s defines a (suboptimal) one-step lookahead policy µ̃(i) = arg min u∈U (i) n " j=1 $ % ˜ pij (u) g(i, u, j)+αJ (j, r) , ∀ i • We will focus mostly on linear architectures J˜(r) = Φr where Φ is an n × s matrix whose columns are viewed as basis functions • Think n: HUGE, s: (Relatively) SMALL • For J˜(r) = Φr, approximation in value space means approximation of J ∗ or Jµ within the subspace S = {Φr | r ∈ ,s } 5 APPROXIMATE VI • Approximates sequentially Jk (i) = (T k J0 )(i), k = 1, 2, . . ., with J˜k (i, rk ) • The starting function J0 is given (e.g., J0 ≡ 0) • After a large enough number N of steps, J˜N (i, rN ) is used as approximation J˜(i, r) to J ∗ (i) • Fitted Value Iteration: A sequential “fit” to produce J˜k+1 from J˜k , i.e., J˜k+1 ≈ T J˜k or (for a single policy µ) J˜k+1 ≈ Tµ J˜k − For a “small” subset Sk of states i, compute n " (T J˜k )(i) = min u∈U (i) pij (u) g(i, u, j) + αJ˜k (j, r) $ j=1 % − “Fit” the function J˜k+1 (i, rk+1 ) to the “small” set of values (T J˜k )(i), i ∈ Sk − Simulation can be used for “model-free” implementation • Error Bound: If the fit is uniformly accurate within δ > 0 (i.e., maxi |J˜k+1 (i) − T J˜k (i)| ≤ δ), lim sup max $ k→∞ i=1,...,n % ∗ ˜ Jk (i, rk ) − J (i) ≤ 6 2αδ (1 − α)2 AN EXAMPLE OF FAILURE • Consider two-state discounted MDP with states 1 and 2, and a single policy. − Deterministic transitions: 1 → 2 and 2 → 2 − Transition costs ≡ 0, so J ∗ (1) = J ∗ (2) = 0. • Consider approximate VI &scheme that approxi' mates cost functions in S = (r, 2r) | r2∈ , 3 with 1 a weighted least squares fit; here Φ = 2 • Given Jk = (rk , 2rk ), we find Jk+1 = (rk+1 , 2rk+1 ), where for weights ξ1 , ξ2 > 0, rk+1 is obtained as rk+1 4 $ %2 $ %2 5 = arg min ξ1 r−(T Jk )(1) +ξ2 2r−(T Jk )(2) r • With straightforward calculation rk+1 = αβr k , where β = 2(ξ1 +2ξ2 )/(ξ1 +4ξ2 ) > 1 • So if α > 1/β, the sequence {rk } diverges and so does {Jk }. • Difficulty is that T is a contraction, butΠ T (= least squares fit composed with T ) is not • Norm mismatch problem 7 APPROXIMATE PI Guess Initial Policy Evaluate Approximate Cost Approximate Policy J˜µ (r) = Φr Using Simulation Generate “Improved” Policy µ Evaluation Policy Improvement • Evaluation of typical policy µ: Linear cost function approximation J˜µ (r) = Φr, where Φ is full rank n × s matrix with columns the basis functions, and ith row denoted φ(i)& . • Policy “improvement” to generate µ: n " $ % & µ(i) = arg min pij (u) g(i, u, j) + αφ(j) r u∈U (i) j=1 • Error Bound: If max |J˜µk (i, rk ) − Jµk (i)| ≤ δ, i k = 0, 1, . . . The sequence {µk } satisfies $ % ∗ lim sup max Jµk (i) − J (i) ≤ k→∞ i 8 2αδ (1 − α)2 POLICY EVALUATION • Let’s consider approximate evaluation of the cost of the current policy by using simulation. − Direct policy evaluation - Cost samples generated by simulation, and optimization by least squares − Indirect policy evaluation - solving the projected equation Φr = ΠTµ (Φr) where Π is projection w/ respect to a suitable weighted Euclidean norm Jµ Tµ(Φr) Projection on S Projection on S Φr = ΠTµ(Φr) ΠJµ 0 S: Subspace spanned by basis functions 0 S: Subspace spanned by basis functions Indirect method: Solving a projected form of Bellman’s equation Direct Mehod: Projection of cost vector Jµ • Recall that projection can be implemented by simulation and least squares 9 WEIGHTED EUCLIDEAN PROJECTIONS • Consider a weighted Euclidean norm 6 7 n 7" $ %2 8 ξi J(i) , )J)ξ = i=1 where ξ is a vector of positive weights ξ1 , . . . , ξn . • Let Π denote the projection operation onto S = {Φr | r ∈ ,s } with respect to this norm, i.e., for any J ∈ ,n , ΠJ = Φr∗ where r∗ = arg mins )J − Φr)2ξ r∈' 10 PI WITH INDIRECT POLICY EVALUATION Guess Initial Policy Evaluate Approximate Cost J˜µ (r) = Φr Using Simulation Generate “Improved” Policy µ Approximate Policy Evaluation Policy Improvement • Given the current policy µ: − We solve the projected Bellman’s equation Φr = ΠTµ (Φr) − We approximate the solution Jµ of Bellman’s equation J = Tµ J with the projected equation solution J˜µ (r) 11 KEY QUESTIONS AND RESULTS • Does the projected equation have a solution? • Under what conditions is the mapping ΠTµ a contraction, so ΠTµ has unique fixed point? • Assuming ΠTµ has unique fixed point Φr∗ , how close is Φr∗ to Jµ ? • Assumption: The Markov chain corresponding to µ has a single recurrent class and no transient states, i.e., it has steady-state probabilities that are positive N 1 " P (ik = j | i0 = i) > 0 ξj = lim N →∞ N k=1 • Proposition: (Norm Matching Property) (a) ΠTµ is contraction of modulus α with respect to the weighted Euclidean norm ) · )ξ , where ξ = (ξ1 , . . . , ξn ) is the steady-state probability vector. (b) The unique fixed point Φr∗ of ΠTµ satisfies )Jµ − Φr∗ )ξ 1 ≤√ )Jµ − ΠJµ )ξ 2 1−α 12 PRELIMINARIES: PROJECTION PROPERTIES • Important property of the projection Π on S with weighted Euclidean norm ) · )ξ . For all J ∈ ,n , J ∈ S, the Pythagorean Theorem holds: )J − J)2ξ = )J − ΠJ)2ξ + )ΠJ − J)2ξ Proof: Geometrically, (J − ΠJ) and ( ΠJ− J) are orthogonal in the scaled geometry of the norm ) · )ξ ,-where two vectors x, y ∈ ,n are orthogonal n if i=1 ξi xi yi = 0. Expand the quadratic in the RHS below: )J − J)2ξ = )(J − ΠJ) + ( ΠJ − J))2ξ • The Pythagorean Theorem implies that the projection is nonexpansive, i.e., )ΠJ − ΠJ¯)ξ ≤ )J − J¯)ξ , for all J, J¯ ∈ ,n . To see this, note that 92 9 92 9 92 9 9Π(J − J)9 ≤ 9Π(J − J)9 + 9(I − Π)(J − J)9 ξ ξ ξ = )J − J )2ξ 13 PROOF OF CONTRACTION PROPERTY • Lemma: If P is the transition matrix of µ, z ∈ ,n )P z)ξ ≤ )z)ξ , Proof: Let pij be the components of P . For all z ∈ ,n , we have )P z)2ξ = = n " i=1  ξi  n " n " n " j =1 2 pij zj  ≤ ξi pij zj2 = j=1 i=1 n " j=1 n " i=1 ξi n " pij zj2 j=1 ξj zj2 = )z)2ξ , where the inequality follows from the convexity of the quadratic function, and the next to equal-last n ity follows from the defining property i=1 ξi pij = ξj of the steady-state probabilities. • Using the lemma, the nonexpansiveness of Π, and the definition Tµ J = g + αP J, we have #ΠTµ J−ΠTµ J̄#ξ ≤ #Tµ J−Tµ J̄#ξ = α#P (J−J̄ )#ξ ≤ α#J−J̄ #ξ for all J, J¯ ∈ ,n . Hence ΠTµ is a contraction of modulus α. 14 PROOF OF ERROR BOUND • LetΦ r∗ be the fixed point of ΠT . We have )Jµ − Φr∗ )ξ 1 ≤√ )Jµ − ΠJµ )ξ . 2 1−α Proof: We have )Jµ − Φr∗ )ξ2 9 92 = )Jµ − + 9ΠJµ − Φr∗ 9ξ 92 9 2 = )Jµ − ΠJµ )ξ + 9ΠT Jµ − ΠT (Φr∗ )9ξ ΠJµ )ξ2 ≤ )Jµ − ΠJµ )2ξ + α2 )Jµ − Φr∗ )2ξ , where − The first equality uses the Pythagorean Theorem − The second equality holds because Jµ is the fixed point of T and Φr∗ is the fixed point of ΠT − The inequality uses the contraction property of ΠT . Q.E.D. 15 MATRIX FORM OF PROJECTED EQUATION • Its solution is the vector J = Φr∗ , where r∗ solves the problem 92 9 ∗ 9 mins Φr − (g + αP Φr )9ξ . r∈' • Setting to 0 the gradient with respect to r of this quadratic, we obtain Φ& Ξ $ Φr∗ − (g + αP Φr∗ ) % = 0, where Ξ is the diagonal matrix with the steadystate probabilities ξ1 , . . . , ξn along the diagonal. • This is just the orthogonality condition: The error Φr∗ − (g + αP Φr∗ ) is “orthogonal” to the subspace spanned by the columns of Φ. • Equivalently, Cr∗ = d, where C = Φ& Ξ(I − αP )Φ, 16 d = Φ& Ξg. PROJECTED EQUATION: SOLUTION METHODS • • Matrix inversion: r∗ = C −1 d Projected Value Iteration (PVI) method: Φrk+1 = ΠT (Φrk ) = Π(g + αP Φrk ) Converges to r∗ because ΠT is a contraction. Value Iterate T(Φrk) = g + αPΦrk Projection on S Φrk+1 Φrk 0 S: Subspace spanned by basis functions • PVI can be written as: rk+1 9 92 = arg mins 9Φr − (g + αP Φrk )9ξ r∈' By setting to 0 the gradient with respect to r, Φ& Ξ which yields $ % Φrk+1 − (g + αP Φrk ) = 0, rk+1 = rk − (Φ& ΞΦ)−1 (Crk − d) 17 SIMULATION-BASED IMPLEMENTATIONS • Key idea: Calculate simulation-based approximations based on k samples Ck ≈ C, dk ≈ d • Matrix inversion r∗ = C −1 d is approximated by r̂k = Ck−1 dk This is the LSTD (Least Squares Temporal Differences) Method. • PVI method rk+1 = rk − (Φ& ΞΦ)−1 (Crk − d) is approximated by rk+1 = rk − Gk (Ck rk − dk ) where Gk ≈ (Φ& ΞΦ)−1 This is the LSPE (Least Squares Policy Evaluation) Method. • Key fact: Ck , dk , and Gk can be computed with low-dimensional linear algebra (of order s; the number of basis functions). 18 SIMULATION MECHANICS • We generate an infinitely long trajectory (i0 , i1 , . . .) of the Markov chain, so states i and transitions (i, j) appear with long-term frequencies ξi and pij . • After generating the transition (it , it+1 ), we compute the row φ(it )& of Φ and the cost component g(it , it+1 ). • We form 1 Ck = k+1 k " t=0 $ φ(it ) φ(it )−αφ(it+1 ) %$ ≈ Φ$ Ξ(I−αP )Φ k 1 " dk = φ(it )g(it , it+1 ) ≈ Φ& Ξg k + 1 t=0 Also in the case of LSPE k 1 " Gk = φ(it )φ(it )& ≈ Φ& ΞΦ k + 1 t=0 • Convergence based on law of large numbers. • Ck , dk , and Gk can be formed incrementally. Also can be written using the formalism of temporal differences (this is just a matter of style) 19 OPTIMISTIC VERSIONS • Instead of calculating nearly exact approximations Ck ≈ C and dk ≈ d, we do a less accurate approximation, based on few simulation samples • Evaluate (coarsely) current policy µ, then do a policy improvement • This often leads to faster computation (as optimistic methods often do) • Very complex behavior (see the subsequent discussion on oscillations) • The matrix inversion/LSTD method has serious problems due to large simulation noise (because of limited sampling) • LSPE tends to cope better because of its iterative nature • A stepsize γ ∈ (0, 1] in LSPE may be useful to damp the effect of simulation noise rk+1 = rk − γGk (Ck rk − dk ) 20 MULTISTEP METHODS • Introduce a multistep version of Bellman’s equation J = T (λ) J, where for λ ∈ [0, 1), ∞ " λ" T "+1 T (λ) = (1 − λ) "=0 Geometrically weighted sum of powers of T . • Note that T " is a contraction with modulus α" , with respect to the weighted Euclidean norm )·)ξ , where ξ is the steady-state probability vector of the Markov chain. • Hence T (λ) is a contraction with modulus αλ = (1 − λ) ∞ " α"+1 λ" "=0 α(1 − λ) = 1 − αλ Note that αλ → 0 as λ → 1 • T t and T (λ) have the same fixed point Jµ and )Jµ − Φrλ∗ )ξ ≤> 1 1− αλ2 )Jµ − ΠJµ )ξ where Φrλ∗ is the fixed point of ΠT (λ) . • The fixed point Φrλ∗ depends on λ. 21 BIAS-VARIANCE TRADEOFF . Solution of projected equation Φ ∗ Φr = ΠT (λ) (Φr) Slope Jµ )λ=0 Simulation error ΠJµ λ=10 Simulation error Simulation error Bias Simulation error Solution of Subspace S = {Φr | r ∈ "s } Set • Error bound )Jµ −Φrλ∗ )ξ ≤ √ 1 1−α2λ )Jµ −ΠJµ )ξ • As λ ↑ 1, we have αλ ↓ 0, so error bound (and the quality of approximation) improves as λ ↑ 1. In fact lim Φrλ∗ = ΠJµ λ↑1 • But the simulation noise in approximating ∞ " T (λ) = (1 − λ) λ" T "+1 "=0 increases • Choice of λ is usually based on trial and error 22 MULTISTEP PROJECTED EQ. METHODS • The projected Bellman equation is Φr = ΠT (λ) (Φr) • In matrix form: C (λ) r = d(λ) , where C (λ) = Φ& Ξ with $ I− P (λ) = (1 − λ) ∞ " αP (λ) % Φ, d(λ) = Φ& Ξg (λ) , α" λ" P "+1 , g (λ) = "=0 α " λ" P " g "=0 • The LSTD(λ) method is $ (λ) %−1 (λ) Ck dk , (λ ) ∞ " (λ ) where Ck and dk are simulation-based approximations of C (λ) and d(λ) . • The LSPE(λ) method is rk+1 = rk − (λ) γGk Ck rk $ − (λ) % dk where Gk is a simulation-based approx. to (Φ& ΞΦ)−1 • TD(λ): An important simpler/slower iteration [similar to LSPE(λ) with Gk = I - see the text]. 23 MORE ON MULTISTEP METHODS (λ ) (λ) • The simulation process to obtain Ck and dk is similar to the case λ = 0 (single simulation trajectory i0 , i1 , . . . more complex formulas) (λ ) Ck k k " $ %& 1 " m−t m−t φ(im )−αφ(im+1 ) , = φ(it ) α λ k + 1 t=0 m=t (λ) dk k k " 1 " αm−t λm−t gim = φ(it ) k + 1 t=0 m=t • In the context of approximate policy iteration, we can use optimistic versions (few samples between policy updates). • Many different versions (see the text). • Note the λ-tradeoffs: (λ) (λ) − As λ ↑ 1, Ck and dk contain more “simulation noise”, so more samples are needed for a close approximation of rλ (the solution of the projected equation) − The error bound )Jµ −Φrλ )ξ becomes smaller − As λ ↑ 1,Π T (λ) becomes a contraction for arbitrary projection norm 24