\documentclass[12pt,amstags,fleqn]{article} % Uncomment only one of these lines: \usepackage[mathlf,textlf,minionint]{MinionPro} \usepackage[T1]{fontenc} \usepackage{textcomp} \usepackage{amsthm} %\usepackage{amssymb} \usepackage{url} \usepackage{xspace} \usepackage{fancyvrb} \usepackage{color} \include{pygments} \usepackage{amsmath} \usepackage{amsthm} \theoremstyle{plain} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{conjecture}[theorem]{Conjecture} \theoremstyle{definition} \newtheorem{definition}[theorem]{Definition} \newtheorem{construction}[theorem]{Construction} \newtheorem{convention}[theorem]{Convention} % CH: I like Remarks to be in normal font, not italics. \theoremstyle{definition} \newtheorem{example}[theorem]{Example} \newtheorem{remark}[theorem]{Remark} \newtheorem{question}[theorem]{Question} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This stuff is needed to correctly render the metapost % diagrams depending on whether latex or pdflatex was % called. Don't change any of it! \usepackage{ifpdf} \ifpdf \pdfcompresslevel=9 \usepackage[pdftex]{graphicx} %\usepackage{thumbpdf} \usepackage[pdftex,colorlinks,bookmarks,backref]{hyperref} %\usepackage[pdftex,backref]{hyperref} \DeclareGraphicsRule{*}{mps}{*}{} \else %\usepackage{hyperref} %\usepackage[dvips]{graphicx} %\usepackage[dvips]{color} \usepackage{epsfig} \usepackage{graphicx} \DeclareGraphicsRule{*}{eps}{*}{} \fi %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \title{Markov Decision Processes and \\the UCT algorithm} \author{Carlo H\"{a}m\"{a}l\"{a}inen\\ {\texttt{carlo.hamalainen@gmail.com}}} %\date{} \begin{document} \maketitle \begin{abstract} In this note we introduce the concept of a Markov Decision Process. We show how to solve for the optimal policy and value vectors for two examples (iPod shuffle and sailing shortest path) using value iteration. We then show how to use the UCT algorithm on the same two problems. Complete source code (in Python) is available. \end{abstract} \section{Introduction} This file and related source code is available at\\ \url{http://carlo-hamalainen.net/stuff/mdpnotes} \\ An MDP (Markov Decision Process) is fully described by the following items: \begin{itemize} \item A set of states $S$. Here we will only consider the case that $S$ is finite, but it may in general be infinite. \item A set of actions $A$. An action is a function that takes us from some state to a new state. \item For each action $a \in A$ and state $s \in S$, the transition probability $P_a(s,s')$ is the probability of moving from state $s$ to $s'$ under action $a$. \item For each state $s$, the (local) cost of taking action $a \in A$ is denoted $\textnormal{cost}(s,a)$. \end{itemize} Here we will assume that the cost is fixed for each state-action pair $(s,a)$. (It is possible to extend the definition to stochastic costs.) If we reach a {\em target state} $t \in S$ then we stop. The goal is to find a {\em policy} $\pi$ that gives the optimal action $\pi(s) \in A$ for each state. The Markov property means that the transition probabilities $P_a(s,s')$ only depend on the current state $s$ and the action $a$ (so no information about past states is used). Once we have an optimal policy $\pi$ we can use it in a simulation as shown in the pseudo-Python of Figure~\ref{mdpsim}. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[htb] \hrule ~ \begin{verbatim} s = initial_state this_cost = 0 while s != target_state: best_action = pi[s] this_cost += cost(s, best_action) s = next_state(s, best_action) print this_cost \end{verbatim} ~ \hrule \caption{Using an optimal policy $\pi$ to run a simulation of an MDP.}\label{mdpsim} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The \texttt{next\_state} function assumes that we have available a {\em generative model} of the MDP that allows us to sample from the possible next states given a state and action. In particular we do not need to know the explicit transition probabilities $P_a(s,s')$. \section{iPod example} A ``real world'' example~\cite{ipodnorvig} comes from the iPod shuffle, an MP3 player with two modes: sequential (in order) play and shuffle. There is no screen so finding a particular song can be difficult -- you can either jump around randomly or try to move forwards (or backwards) in a linear manner. We can assume that the $N$ songs on our iPod have been sorted in some sensible way, and we'll label the songs $0$, $1$, \dots, $N-1$. There are two actions that we can take at any given song (state) \begin{itemize} \item Sequential: set the iPod to non-random playback and use the forward or backward track buttons to get to the target song. If we are at song $s$ and the target song is $t$, then this strategy will take $\left| s - t \right|$ button presses to find the target. \item Shuffle: set the iPod to random playback and move to the next track if we are not at the target song. This incurs a cost $T$ because it will take us some time to recognise the new song. \end{itemize} In a typical execution of this algorithm we will start at some initial song, randomly skip around a bit, and then use the sequential action to get to the final song. Or, if the target song is very close, we would only use the sequential strategy. Suppose that the iPod has $N$ songs. The set of states for this MDP is \[ S = \{0, 1, \dots, N-1\}. \] There are two actions, so \[ A = \{ \textnormal{sequential}, \textnormal{shuffle} \}. \] For the costs, the sequential strategy takes $\left| s - t \right|$ button presses from state $s$, so \[ \textnormal{cost}(s, \textnormal{sequential}) = \left| s - t \right|. \] The shuffle action only takes us to a different track but we have to recognise if this is the target song, so \[ \textnormal{cost}(s, \textnormal{shuffle}) = T \] for some constant $T$. Finally we have the transition probabilities. The sequential action is guaranteed to take us to the target song, so $P_{\textnormal{sequential}}(s,s') = 1$ if $s' = t$ and $0$ otherwise. The shuffle action merely takes us to another track, with equal probability, so $P_{\textnormal{shuffle}}(s,s') = 1/N$ for any state $s'$. \section{Value iteration} To find the optimal policy $\pi$ we will compute the value vector $V$. For each state $s \in S$, $V(s)$ is the expected cost of reaching the target state, using the best possible sequence of actions starting at state $s$. We can start from an initial value vector $V_0(s) = 0$ for all $s \in S$. Then the update step is \[ V_{t+1}(s) = \min_{a \in A} \left\{ \textnormal{cost}(s,a) + \sum_{s' \in S} P_a(s,s') V_t(s')\right\}. \] In other words, we minimise the sum of the local cost of taking some action, along with the expected cost from the possible new states. The policy is updated at each step by \[ \pi_{t+1}(s) = \underset{a \in A}{\operatorname{argmin}} \left\{ \textnormal{cost}(s,a) + \sum_{s' \in S} P_a(s,s') V_t(s')\right\}. \] The value iteration algorithm is guaranteed to converge~\cite{viconvergence}. We can also apply a discounting factor $\gamma$ to reflect the fact that costs (or rewards) incurred in the future are worth less than ones incurred in the present. If $0 \leq \gamma \leq 1$ then the value iteration update is simply \[ V_{t+1}(s) = \min_{a \in A} \left\{ \textnormal{cost}(s,a) + \gamma \sum_{s' \in S} P_a(s,s') V_t(s')\right\}. \] \subsection{Value iteration for the iPod example} See Section~\ref{secipodpython} for the Python source code for value iteration on the iPod example. In this example we have $\gamma = 1$, so there is no discounting. With $N = 10$ songs, recognition cost of $T = 0.5$, and target song $N/2$, we get the following value and policy vectors: \begin{Verbatim} >>> from ipod_mdp import * >>> N = 10; value_iteration(N, 0.5, N/2) ([2.1984130546875003, 2.1984130546875003, 2.1984130546875003, 2.0, 1.0, 0.0, 1.0, 2.0, 2.1984130546875003, 2.1984130546875003], ['shuffle', 'shuffle', 'shuffle', 'sequential', 'sequential', 'sequential', 'sequential', 'sequential', 'shuffle', 'shuffle'], 2) \end{Verbatim} For example, if the initial state is $1$, an execution may proceed as follows: \begin{itemize} \item $s = 1$, $\pi(1) = \textnormal{shuffle}$. Jump to random song $s = 8$. \item $s = 8$, $\pi(8) = \textnormal{shuffle}$. Jump to random song $s = 3$. \item $s = 3$, $\pi(3) = \textnormal{sequential}$. Click $5-3 = 2$ times to get to the target state $5$. \end{itemize} The cost of this run is $T + T + 2 = 0.5 + 0.5 + 2 = 3$, while the value vector has $V(1) = 2.1992065273437502$, which is reasonably close. A general execution with $N = 250$ and $T = 0.5$: \begin{Verbatim} $ python ipod_mdp.py 250 0.5 mean(V) = 10.6647967086 mean (simulation): 10.679298 difference: -0.0145012913545 shuffle when: 12 or more away \end{Verbatim} \section{Sailing} Original reference: \cite{sail}. Imagine a sailing boat on a rectangular lake. Our goal is to get to the north-east corner of the lake from some starting position. To simplify things we assume that the lake is divided into a grid of {\em waypoints}, and that we sail from one waypoint to another (so the boat can only travel in one of eight directions north, north-east, east, etc). To further simplify the situation, we assume that the wind blows in only one direction for the duration of a journey from one waypoint to another waypoint. The wind changes just as we set off towards a new waypoint. Directions will be the integers $D = \{0, 1, \dots, 7 \}$ where $0$ is north, $1$ is north-east, $2$ is east, etc. A state $s = (x,y,d,w_1,w_2)$ consists of a location $(x,y) \in \mathbb{N} \times \mathbb{N}$, the previous boat direction $d \in D$, previous wind direction $w_1 \in D$, next wind direction $w_2 \in D$. An action is the direction $d \in D$ to sail the boat on the next leg. For each action $a \in A$ and state $s \in S$, the transition probability \[ P_{d'}((x,y,d,w_1,w_2),\,(x',y',d',w_1',w_2')) = 0 \] unless $w_2 = w_1'$, $(x',y') = (x,y) + \vec{d'}$, in which case the probability is $P(w_1', w_2')$, the probability that the wind changes from direction $w_1'$ to $w_2'$. If the boat is in direction $d$ and the wind is in direction $w$ then we have the following costs: \[ \textnormal{cost}(d,w) = \begin{cases} 1\textnormal{ minute}, &\textnormal{if $\alpha = 0$} \\ 2\textnormal{ minutes}, &\textnormal{if $\alpha = 1$} \\ 3\textnormal{ minutes}, &\textnormal{if $\alpha = 2$} \\ 4\textnormal{ minutes}, &\textnormal{if $\alpha = 3$} \end{cases} \] where $\alpha = \min(d-w, 8-(d-w))$. We do not allow the case $\alpha = 4$ which corresponds to sailing directly into the wind. The original description of the sailing problem also includes a cost for changing tack of the boat (i.e. when the wind goes from being on the left side of the sails to the right side or vice versa) but we ignore that cost to simplify the model. The wind transition probabilities are given in an $8 \times 8$ matrix $W$. The probability that the wind will be in direction $w_2$ given that it is currently in direction $w_1$ is $W_{w_1,w_2}$. In our simulations we fix $W$ to the following matrix: \[ W = \left[ {\begin{array}{cccccccc} 0.4 & 0.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.3 \\ 0.4 & 0.3 & 0.3 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.4 & 0.3 & 0.3 & 0.0 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.4 & 0.3 & 0.3 & 0.0 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.4 & 0.2 & 0.4 & 0.0 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.3 & 0.3 & 0.4 & 0.0 \\ 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.3 & 0.3 & 0.4 \\ 0.4 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.3 & 0.3 \\ \end{array} } \right] \] We can now use value iteration to solve for the value and policy vectors. \subsection{Example of sailing} See Section~\ref{pythonsailing} for Python code for value iteration on the sailing example. We run value iteration and 10,000 simulations with a lake of size $5 \times 5$, target location $(4,4)$, and all other parameters as described above. \begin{Verbatim} $ python sailing.py "value_iteration_example()" Top of value_iteration(), max difference: 1.0 Top of value_iteration(), max difference: 4.0 Top of value_iteration(), max difference: 2.89 Top of value_iteration(), max difference: 2.2653 Top of value_iteration(), max difference: 1.935648 Top of value_iteration(), max difference: 1.66572639 Top of value_iteration(), max difference: 1.4497415721 Top of value_iteration(), max difference: 1.25418853297 Top of value_iteration(), max difference: 0.998954591375 Top of value_iteration(), max difference: 0.769139503728 Top of value_iteration(), max difference: 0.510518276155 Top of value_iteration(), max difference: 0.337118955482 Top of value_iteration(), max difference: 0.193325220838 Top of value_iteration(), max difference: 0.0315788407471 Done with value iteration Running simulations... Lake size: 5 x 5 Value iteration: Mean cost: 7.3 Median cost: 7.1 Standard dev: 3.7 Simulations (run 1000 times): Mean cost: 7.3 Median cost: 6.8 Standard dev: 4.1 Mean cost to sail across lake from (1, 1) to (4, 4): 8.5 \end{Verbatim} \section{UCT: Upper Confidence bounds on Trees}\label{uctsec} Solving an MDP system using value iteration can become computationally intensive on large examples (e.g. sailing on a lake of size $40 \times 40$) because each update step necessarily reads and changes every element of the value and policy vectors. An alternative approach is to use Monte Carlo planning. See Figure~\ref{mcplanning} for pseudo-Python. The algorithm alternates randomly between between trying new actions at each state (in order to search for better policies) or using the current best policy to improve the estimate of the policy. The policy vector $\pi$ is not explicitly computed. Instead we have the state-action vector $Q(s,a)$, the average seen cost (reward) of taking action $a$ when in state $s$. The update line \begin{center} \begin{verbatim} new_average = old_average + (0.5)*(q - old_average) \end{verbatim} \end{center} has a bias of $0.5$ instead of $1/n$ and this may need to be tuned for each application of Monte Carlo planning. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[htb] \begin{verbatim} def select_action(state): if random() < 0.01: return random_action(state) else: return best_action(state) def search(state): if state == terminal: return 0 action = select_action(state) new_state, cost = simulate_action(state, action) state_visit_counts[(state)] += 1 state_action_counts[(state, action)] += 1 try: old_average = Q[(state, action)] n = state_action_counts[(state, action)] new_average = old_average + (0.5)*(q - old_average) except KeyError: new_average = q self.Q[(state, action)] = new_average return q def monte_carlo_planning(initial_state): while True: search(initial_state) \end{verbatim} \caption{Pseudo-Python for Monte Carlo planning.}\label{mcplanning} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% There are some subtleties to the implementation. See Section~\ref{pythonmc} for source code. The UCT algorithm~\cite{Kocsis06banditbased} just changes the way that actions are selected during a rollout. The algorithm used is called UCB, and described in the next section. \subsection{UCB: Upper Confidence Bounds} Imagine that we have $K$ gambling machines with arbitrary reward distributions $P_1$, \dots, $P_K$. At each time step we can play any machine $j$ that we choose, and receive a reward according to the distribution $P_j$. We would like a strategy that maximises our total reward over $n$ plays. Since the distributions $P_j$ are fixed but unknown, we want to avoid sampling too many times from machines with low reward. In other words, we have to balance exploration versus exploitation. The UCB1 strategy~\cite{Auer:2002} is: \begin{enumerate} \item Play each machine once. \item Play machine $j$ that maximises $\bar x_j + \sqrt{\frac{2 \ln n}{n_j}}$ where $\bar x_j$ is the average reward from machine $j$, machine $j$ has been played $n_j$ times so far, and $n$ is the total number of plays so far. \end{enumerate} The $\bar x_j$ term gives preference to machines that have performed well in past plays, while the $\sqrt{{2 \ln n}/{n_j}}$ term gives preference to machines that have not been played many times so far, relative to $\ln n$. Write $\mu_j$ for the mean of distribution $P_j$. Obviously the best possible strategy is to only play the machine with the maximum mean $\mu^* = \max_j{\mu_j}$. Figure~\ref{figucbreward} shows best and actual reward for some simulations using the UCB1 strategy. See Section~\ref{pythonucb} for Python source code. %%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[htb] \begin{center} \includegraphics[width=.75\textwidth]{best_and_total_reward.pdf} \end{center} \caption{Comparing best and actual reward using the UCB1 strategy with $10$ machines. See Section~\ref{pythonucb} for details on the machines' reward distributions.}\label{figucbreward} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%% The {\em regret} is the lost reward due to not having played the optimal machine at each step: \[ n \mu^* - \mu_j \sum_{j=1}^K \mathbb{E}(T_j(n)) \] where $T_j(n)$ is the number of times that machine $j$ has been played after $n$ plays. Importantly, the regret can be bounded to be logarithmic in the number of plays so far: \begin{theorem}[\cite{Auer:2002}]\label{thmUCB} For all $K > 1$, if policy \textnormal{UCB1} is run on $K$ machines having arbitrary reward distributions $P_1$, \dots, $P_K$ with support in $[0,1]$ then its expected regret after $n$ plays is at most \[ \left[ 8 \sum_{i:\mu_i < \mu^*} \left( \frac{\ln n}{\Delta_i} \right) \right] + \left( 1 + \frac{\pi^2}{3} \right) \left( \sum_{j=1}^K \Delta_j \right) \] where $\Delta_i = \mu^* - \mu_i$. \end{theorem} Figure~\ref{ucbregret} shows the result of some simulations with $K = 10$ machines and fixed reward distributions. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[htb] \begin{center} \includegraphics[width=.75\textwidth]{regret_and_bound.pdf} \end{center} \caption{Regret and upper bound on regret given by Theorem~\ref{thmUCB}}\label{ucbregret} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{UCT} The UCT algorithm is the same as Monte Carlo planning except that UCB is used to select the action at each tree node. To do this we have to keep track of the average cost so far (i.e. up to time $t$) of taking action $a$ from state $s$, denoted $Q_t(s,a)$. We also track $N_{s}(t)$, the number of times that state $s$ has been visited up to time $t$. If we are at a tree node $s$ then the action is chosen using the UCB rule \begin{equation}\label{eqnuct} a^* = \underset{a \in A}{\operatorname{argmax}} \left\{ Q_t(s,a) + \alpha \sqrt{\frac{\ln N_{s}(t)}{N_{s,a}(t)}} \, \right\}. \end{equation} where $\alpha$ is a constant that has to be chosen empirically. A badly chosen $\alpha$ can have a huge effect on the convergence of the algorithm. See Section~\ref{pythonuct} for Python source code that uses UCT to solve the sailing problem.\footnote{I asked the authors of \cite{Kocsis06banditbased} for their full list of parameter values and source code but got no reply. If anyone improves my code I will be happy to update this document and give a reference or URL.} \subsection{UCT on the sailing problem} First we used value iteration to solve the sailing problem on lakes of size $5 \times 5$, $10 \times 10$, $20 \times 20$, and $30 \times 30$, with $\varepsilon = 0.01$ (solving just the $20 \times 20$ instance took about 5 hours on a 2.5Ghz Intel Xeon server). The optimal value and policy vectors $V^*$ and $\pi^*$ are available in the pickled Python format as {\texttt lake\_5.pkl}, {\texttt lake\_10.pkl}, {\texttt lake\_20.pkl}, and {\texttt lake\_30.pkl}. Next, we chose $20$ random points (excluding the target point). For each of these states $s$, we ran the Monte Carlo and UCT algorithm until the estimated cost for sailing from $s$ to the target state was within $0.1$ of the optimal value $V^*(s)$. A good measure of the efficiency of a Monte Carlo type of planning algorithm is the total number of samples taken from the underlying MDP's generative model. Figure~\ref{uctsailing} shows the number of samples required to get the estimated error to within the bounds specified. %%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{figure}[htb] \begin{center} \includegraphics[width=.75\textwidth]{nr_samples_uct_sailing.pdf} \end{center} \caption{Performance of plain Monte Carlo planning and UCT on the sailing problem.}\label{uctsailing} \end{figure} %%%%%%%%%%%%%%%%%%%%%%%%%%% \section{What else?} There are a few factors that have to be chosen by hand in {\texttt sailing\_uct.py}. If we are working on a huge problem, where value iteration isn't feasible, how do we know if we are converging to the optimal solution? The $0.1$ factor for the $Q$ value update step in Section~\ref{uctsec} is meant to decay to zero if we want convergence. But at what rate? Using something like $1/n$ works badly on the sailing domain with lake sizes up to $20 \times 20$. Also, there are many ways to improve the UCT algorithm, for example combining online and offline data to speed up convergence of the $Q$ values: \cite{citeulike:2976742}. \clearpage \section{iPod value iteration}\label{secipodpython} Here is \texttt{ipod\_mdp.py}, also available at\\ \url{http://carlo-hamalainen.net/stuff/mdpnotes} \input{ipod_mdp.tex} \clearpage \section{Value iteration for sailing}\label{pythonsailing} Here is \texttt{sailing.py}, also available at\\ \url{http://carlo-hamalainen.net/stuff/mdpnotes} \input{sailing.tex} \clearpage \section{UCB}\label{pythonucb} Here is \texttt{ucb.py}, also available at\\ \url{http://carlo-hamalainen.net/stuff/mdpnotes} \input{ucb.tex} \clearpage \section{Sailing MCT}\label{pythonmc} Here is \texttt{sailing\_mc.py}, also available at\\ \url{http://carlo-hamalainen.net/stuff/mdpnotes} \input{sailing_mc.tex} \clearpage \section{UCT}\label{pythonuct} Here is \texttt{sailing\_uct.py}, also available at\\ \url{http://carlo-hamalainen.net/stuff/mdpnotes} \input{sailing_uct.tex} \bibliographystyle{plain} \bibliography{mdp} \end{document}