Run on Google Colab to test out the code!

The purpose of this post is to introduce some basic concepts from the reinforcement learning (RL) literature using the setting and notation that is familiar to economists studying infinite horizon dynamic decision problems. In particular, I make use of Rust (1987)'s bus engine replacement problem as a demonstrative example.

Note: the methods discussed here are mainly concerned with solving the dynamic problem, and not with estimation.

For those unfamiliar with RL, here are some recent breakthrough examples of these methods

• AlphaGo beating human Go champion
• Playing Atari games by only observing raw image input
• Robotic manipulation of arbitrary objects

What is notable about these applications?

They all involve notable complexity:

1. A complex/large state space, as in the case of Go
2. Complex transition dynamics, which are analytically intractible, as in the case of Atari games
3. Complex policy functions, as is the case with the large number of motors involved in robotic manipulation

The reader should keep these more complex applications in mind. This notebook focuses on a simple case for expositional purposes, with some particular features that RL struggles with. The real power of the latter methods is only apparent in much more complex situations, which I hope to revisit in future posts.

---

Setup and Notation¶

The methods presented here are appropriate for finite Markov Decision Processes, characterized by the following:

States: $x \in \mathcal{X}$, finite

Actions: $a \in \mathcal{A}(x)$, finite

Transition probabilities: $p(x'|x,a)$

Flow utility: $u(x,a)$

Discount rate: $\beta$

Policy function: $a^*(x) \in \mathcal{A}(x)$

The agent's problem: $$ \begin{align} \max_{a^*} \mathbb{E}_{x_1,x_2,\dots} &\Big[ \sum_{t=0}^\infty \beta^t u(x,a) | x_0 \Big] \\ &s.t. \\ a_t &= a^*(x_t) \\ Pr(x_{t+1}|x_t) &= p(x_{t+1}|x_t,a^*(x_t)) \\ \end{align} $$

We can write the problem in recursive form using the Bellman equation

$$ V(x) = \max_{a \in \mathcal{A}(x)} u(x,a) + \beta \sum_{x' \in \mathcal{X}} V(x') p(x'|x,a) $$

Bus Engine replacement¶

Recall the problem studied by Rust (1987):

1. Bus engines get worn over time and require increasing amounts of maintenance each month.
2. Alternatively, the entire engine can be replaced at a high fixed cost

Here, the state variable $x$ represents the mileage of the engine. The actions space at every period is $\mathcal{A} = \{0,1\}$, where $1$ implies replacing the engine.

Utilities are:

• $u(x,0) = -c(x)$
• $u(x,1) = -R - c(0)$

with $c(x) = \theta x$. For the remainder, we fix $\theta = 8.6*10^{-4}$ and $R = 8000$.

The transition process for $x$ is s.t.:

• Conditional on $a_t = 0$, $x_{t+1} - x_t \sim$ Exponential($\lambda$)
• Conditional on $a_t = 1$, $x_{t+1} \sim$ Exponential($\lambda$)

This process will be discretized into a 201-point grid for $x \in [0, 300000]$ with a spacing of $1500$ miles in between points.

Rust (1987) solution¶

Rust's approach uses the contraction mapping property of the Bellman equation to iterate over choice-specific value functions (which we will call Q-functions).

Let: $$ Q(x,a) \equiv u(x,a) + \beta \sum_{x' \in \mathcal{X}} V(x') p(x'|x,a) $$

Then we have the following contraction mapping: $$ Q(x,a) = u(x,a) + \beta \sum_{x \in \mathcal{X}} p(x'|x,a) \arg \max_{y \in \mathcal{Y}(x)} Q_{x'y} $$

This is almost the contraction mapping from Rust (1987), simply without the random utility component.

This can be iterated on $Q(x,a)$ to solve the problem.

Algorithm¶

Denote the current iteration by $k$. Initialize $k=0$, $Q^0$

For every $k=0,1,\dots$:

• $Q^{k+1}(x,a) = u(x,a) + \beta \sum_{x \in \mathcal{X}} p(x'|x,a) \arg \max_{a \in \mathcal{A}(x)} Q^k(x',a), \quad \forall (x,a)$
• Stop when $|Q^k - Q^{k-1}| < \delta$

Key point: notice that at every iteration $k$, the update is performed over $|\mathcal{X}|\cdot|\mathcal{A}|$ points.

Computational performance¶

In terms of compute time, the classic method is extremely fast, despite having to perform what seems like a large (>200K) number of updates. Partly this can be attributed to the ease of executing the algorithm in parallel with very little overhead.

Benchmark solution¶

By leveraging the contraction mapping property of the Bellman equation and by iterating over all $(x,a)$ pairs equally, we can be sure that the solution we have obtained is correct up to our specified error tolerance.

To benchmark the other solutions, we can compare:

1. The difference in payoffs between replacing and not replacing as implied by the computed Q-functions (upper plot).
2. The resulting policy function or argmax (lower plot).

---

Q-learning¶

Q-learning, introduced by Watkins (1989) is another way of solving such an MDP.

An intuitive description of Q-learning is "learning by doing". It makes use of simulation to iterate on the equation. This affects both which $(x,a)$ pairs are updated at each iteration and how they are updated.

We can break down this approach into two elements: asynchronous updating and stochastic iteration

---

Asynchronous (and stochastic) update path¶

Whereas Rust (1987) updates every $(x,a)$ at every iteration, the basic Q-learning iteration only updates a single $(x,a)$ per iteration. The choice of $(x,a)$ is determined by simulation, so $y$ is chosen according to the current policy, while the new $x$ is then chosen according to the transition probabilities.

Algorithm¶

Denote the current iteration by $k$. Initialize $k=0$, $x^0$, $Q^0$

For every $k=0,1,\dots, K$:

• $a^k = \arg \max Q^k_(x^k,a)$
• Draw $x^{k+1}$ according to $p(x^{k+1}|x^k,a)$
• Update:

  • If $(x,a) = (x^k,a^k)$: $Q^{k+1}(x,a) = u(x,a) + \beta \sum_{x \in \mathcal{X}} p(x'|x,a) \arg \max_{y \in \mathcal{Y}(x)} Q^k(x',a)$
  • Otherwise: $Q^{k+1}_{xy} = Q^k(x,a)$

• Non-standard: check for convergence by keeping a running history of the size of function updates

Key point: note that, in contrast to the previous solution, the update is only performed at one point in $(x,a)$ space at every iteration. In principle, this leads to a more efficient use of information, since the update of one $(x,a)$ pair can immediately be used in the update of another $(x',a')$ pair instead of using 'outdated' information. This is a general feature of asynchronous iterative methods and usually comes at the cost of reduced parallelizability.

1st issue: Exploration¶

Picking $a^k = \arg \max Q^k(x^k,a)$ has the advantage of naturally improving our estimate of $Q$ at the points that occurr most frequently in reality.

However, in the beginning, our guess for $Q$ may be off substantially, and we may fail to pick a $\tilde{a}$ because the current guess for $Q(x,\tilde{a})$ is too low. Since we never pick it, the value can never be updated, and the iteration gets stuck in a suboptimal path.

To fix this, we need to incoporate some form of exploration.

The simplest form of this is to pick $a^k$ according to the $\epsilon$-greedy policy:

• With probability $1-\epsilon$: $a^k = \arg \max Q^k(x^k,a)$
• With probability $\epsilon$: pick $a^k$ randomly

Depending on the transition matrix, some exploration in terms of the state space $\mathcal{X}$ may also be needed. This is apparent in the bus engine replacement problem, since higher-mileages can only be observed after passing through the lower mileages, so estimates at these states can be very noisy. A simple way to accomplish this is by periodically resetting the state variable randomly.

Solution comparison¶

Comparing the solution with the benchmark, we see that the solution generally matches very well, except for a few points. This demonstrates the stochastic nature of the update process and the consequent difficulty of verifying convergence.

Furthermore, notice that the total number of function updates is comparable to the synchronous and parallel update method, while being much slower due to lack of parallelization and general overhead from executing a loop.

The lackluster improvement in terms of the number of function updates has partly to do with the nature of the problem, which results in heavily skewed sampling of the $(x,a)$ space.

Here, we re-run the algorithm but set the reset frequency to 1, such that the state $x$ is randomized at every iteration, yields a much stronger improvement in the number of function updates required.

---

Stochastic iteration¶

Now reconsider the term: $p(x'|x,a) \arg \max_{y \in \mathcal{Y}(x)} Q^k_{x'y}$. In some cases, this can be very expensive to compute, in particular when the transition probabilities $p(x'|x,a)$ are intractable or take on a complex form.

Recall that we are already drawing the next state $x^{k+1}$ from $p(x'|x,a)$ anyways. We can use this to get a very noisy estimate of the above expectation, and use averaging over many periods to smooth out the noise, sing the same principle as stochastic optimization.

Essentially, when we visit $(x^k,a^k)$ in iteration $k$ and get a draw of the following state $x^{k+1}$, this gives a noisy estimate of $Q(x^k,a^k)$:

$$\tilde{Q}(x^k,a^k) = u(x^k,a^k) + \beta \arg \max_{a'} Q^k(x^{k+1},a')$$

Because this estimate is noisy, we want to average it over multiple visits of $(x,a)$, so we do not fully update $Q$ to the new value.

Algorithm¶

Initialize $k=0$, $x^0$, $Q^0$.

For every $k=0,1,\dots$:

• $a^k = \arg \max Q^k(x^k,a)$
• Draw $x^{k+1}$ according to $p(x'|x^k,a^k)$
• Get the target value: $\tilde{Q}(x^k,a^k) = u(x^k,a^k) + \beta \arg \max_{a'} Q^k(x^{k+1},a')$
• Update:

  • $Q^{k+1}(x,a) = Q^{k}(x,a) + \alpha^k (\tilde{Q}(x^k,a^k) - Q^k(x,a))$ if $(x,a) = (x^k,a^k)$
  • $Q^{k+1}(x,a) = Q^k(x,a)$ otherwise

Key point: note that the updates do not require integrating over the future possible states $x'$. In this example with a pre-computed Markov transition matrix and a small state space, this provides little to no advantage. However, in situations with a very large state space and expensive to compute transition probabilities, this can be a significant advantage.

2nd issue: Convergence¶

With the new form of updating, we no longer have the convergence guarantees of the contraction mapping.

Instead, $Q^k$ will behave more like a random variable. If $\alpha^k$ is held fixed, it will never completely converge, only bounce around in the neighborhood of the true $Q$.

However, there are conditions on the sequence of $\alpha^k$ which should guarantee convergence:

1. $\sum_{k=0}^\infty \alpha^k = \infty$
2. $\sum_{k=0}^\infty \alpha^k < \infty$

Sutton and Barto (2020), p.33

In practice, these conditions are not strictly observed, particularly because they can lead to very slow convergence rates. However, for econometric applications, guaranteed convergence may be more important.

Solution comparison¶

We can see here that the computed Q-functions approximate the benchmark fairly closely, though the noise due to stochastic iteration is apparent, even leading to an incorrect policy function in some states.

However, one should also note that this was accomplished with roughly 5 times as many updates, even though we were previously integrating over a state space of $|\mathcal{X}| = 201$. In cases where the computation of transition probabilities is the main computational burden, this represents a speedup factor of roughly 40!

---

Deep Q-learning¶

So far, the algorithms presented use a large table to represent the Q-values (they are all tabular algorithms). In cases with large state and actions spaces, this implies that (a) a large table is needed to store the values (poor scaling in space complexity) and (b) a large number of updates is needed to reach convergence everywhere.

The recent breakthroughs linked at the beginning of this notebook were only made possible by combining Q-learning with state of the art deep neural networks, which can represent the Q-values much more sparsely and are generally flexible enough to approximate whatever underlying structure the Q-values have.

In general, any universal function approximator will work, as we typically expect the Q-values to be smooth enough in $(x,a)$ such that neighboring pairs will have similar Q-values. Using a function approximator will then also speed up convergence, since a single update can affect a larger region of the $(x,a)$ space.

Finally, the use of such an approximation also allows one to dispense with the discretization of a continuous state space.

In the following, let $Q^\theta$ denote an approximation to $Q$ parametrized by $\theta$.

Algorithm¶

Initialize $k=0$, $x^0$, $\theta^0$.

For every $k=0,1,\dots$:

• $a^k = \arg \max Q^{\theta^k}(x^k,a)$
• Draw $x^{k+1}$ according to $p(x'|x^k,a^k)$
• Get the target value: $\tilde{Q}^k = u(x^k,a^k) + \beta \arg \max_{a'} Q^{\theta^k}(x^{k+1},a') $
• Update:

  • $\theta^{k+1} = \theta^k + \alpha^k \nabla_\theta \big( \tilde{Q}^k - Q^{\theta^k}(x^k,a^k) \big)$

Notice: the update step corresponds to gradient descent of the problem: $\arg \min_\theta \big(\tilde{Q}^k - Q^{\theta^k}_{x^ka^k}\big)^2 $

Initial parameters¶

The parameters are initialized by the pytorch library. For context, here are the Q-functions given by the initial neural network.

Solution comparison¶

We can immediately see that the approximation excessively smooths out the value functions. In principle this can be remedied with a larger number of iterations and smaller step-sizes, at the cost of increased computing time. This essentially echoes the bias-variance trade-off that is found throughout machine learning and non-parametric estimation.

However, note that even with a clear bias in the value functions, the resulting policy function is not far off, with replacement occurring a bit earlier. This demonstrates that large biases in the value functions need not translate into the same level of biases in the policy functions.

Finally, note that this solution was obtained after only 30,000 function updates. When compared to the tabular Q-learning solution obtained after 1 million iterations, this represents a speedup factor of 30!

---

Conclusion¶

The objective of this notebook was to give a brief but concrete introduction to some of the methods used in the modern reinforcement learning literature.

The key advantages are:

• Less iterations required to reach an adequate solution
• Far less computations of transition probabilities
• Can be used with a simulation model without analytical transition probabilities
• Naturally combines well with function approximation to deal with large state-action spaces

The key disadvantages are:

• Convergence is no longer guaranteed and can be difficult to verify
• Parallelization is possible but less trivial
• Function approximation introduces bias

As a reminder, the example problem here is used for its familiarity to economists, but certainly does a poor job of showing the advantage of reinforcement learning techniques. In future notebooks, I hope to provide other examples that better showcase the power of these tools.

---

References¶

Setup and notation:

• Rust, J. (1987). Optimal Replacement of GMC Bus Engines: An Empirical Model of Harold Zurcher. Econometrica, 55(5), 999-1033.

Papers in economics:

• Iskhakov, F., Rust, J., & Schjerning, B. (2020). Machine learning and structural econometrics: contrasts and synergies. The Econometrics Journal, 23(3), S81-S124.
• Igami, M. (2020). Artificial intelligence as structural estimation: Deep Blue, Bonanza, and AlphaGo. The Econometrics Journal, 23(3), S1-S24.

Reference textbook:

• Sutton, R. S., & Barto, A. G. (2018). Reinforcement learning: An introduction. MIT press.

Web resources:

• CS 285 at Berkeley, taught by Sergey Levine
• Spinning up by OpenAI