This code is provided under Creative Commons Attribution-ShareAlike 4.0 International (CC BY-SA 4.0) License.
AIMA chapter 17 is about sequential decision problems where the agent’s utility depends on a sequence of decisions. We will implement the key concepts using R for the AIMA 3x4 grid world example.
More details can be found in Reinforcement Learning: An Introduction (RL) by Sutton and Barto (2020) in Chapter 4: Dynamic Programming.
Note that the code in this notebook defines explicit functions matching the textbook definitions and is for demonstration purposes only. Efficient implementations for larger problems use fast vector multiplications instead.
MDPs are sequential decision problems with
MDPs are defined by:
A policy \(\pi = \{\pi(s_0), \pi(s_1), \dots\}\) defines for each state which action to take. If we assume that under policy \(\pi\), the agent will go through the state sequence \(s_0, s_1, ..., s_\infty\), then the expected utility of being in state \(s_0\) can be calculated as a sum. To incorporate that earlier rewards are more important, a discount factor \(\gamma\) is used.
\(U^\pi = \mathbb{E}\left[\sum_{t=0}^\infty \gamma^t R(s_t, \pi(s_t), s_{t+1})\right]\)
The goal of solving a MDP is to find an optimal policy that maximizes the expected future utility.
\(\pi^*(s) = \mathrm{argmax}_\pi U^\pi(s)\) for all \(s \in S\)
The used example is the simple 4x3 grid world described in AIMA Figure 12.1 and used again in 22.1 as:
In the following, we implement states, actions, the reward function,and the transition model. We also show how to represent a policy and how to estimate the expected utility using simulation.
The MDP will be defined as the following global variables/functions:
S: set of states.A: set of actions.actions(s): returns the available actions for state \(s\).P(sp, s, a): a function returning the transition probability \(P(s' | s, a)\).R(s, a, s_prime): reward function for the transition from \(s\) to \(s`\) with action a.Policies are represented as:
Other useful functions:
sample_transition(s, a): returns a \(s'\) sampled using the transition model.simulate_utilities(pi, s0 = 1, N = 1000, max_t = 100): Estimates the utility of following policy \(\pi\), starting \(N\) episodes in state \(s_0\). The maximal episode length is max_t to ensure that the function finishes also for policies that do not lead to a terminal state.We define the atomic state space \(S\) by labeling the states \(1, 2, ...\). We convert coordinates (rows, columns) to the state label.
# I use capitalized variables as global constants
COLS <- 4
ROWS <- 3
S = seq_len(ROWS * COLS)
LAYOUT <- matrix(S, nrow = ROWS , ncol = COLS)
LAYOUT [,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12Note that the rows are displayed upside-down compared to the text book, so we use a function to display them in reverse order.
show_layout <- function(x) {
x <- matrix(x, ncol = COLS, nrow = ROWS,
dimnames = list(row = seq_len(ROWS), col = seq_len(COLS)))
x[rev(seq_len(ROWS)), ]
}
show_layout(LAYOUT) col
row 1 2 3 4
3 3 6 9 12
2 2 5 8 11
1 1 4 7 10Convert between coordinates and state labels.
rc_to_s <- function(rc)
LAYOUT[rbind(rc)]
s_to_rc <-
function(s)
drop(which(LAYOUT == s, arr.ind = TRUE, useNames = FALSE))
rc_to_s(c(3, 4))[1] 12[1] 3 4Start state
Define terminal states.
The complete set of actions is \(A = \{\mathrm{'Up', 'Right', 'Down', 'Left', 'None'}\}\). Not all actions are available in every state. Also, action None is added as the only possible action in an absorbing state.
A = c('Up', 'Right', 'Down', 'Left', 'None')
actions <- function(s) {
# absorbing states
if(s == 11 || s == 12) return('None')
# illegal state
if(s == 5) return('None')
c('Up', 'Right', 'Down', 'Left')
}
lapply(S, actions)[[1]]
[1] "Up" "Right" "Down" "Left"
[[2]]
[1] "Up" "Right" "Down" "Left"
[[3]]
[1] "Up" "Right" "Down" "Left"
[[4]]
[1] "Up" "Right" "Down" "Left"
[[5]]
[1] "None"
[[6]]
[1] "Up" "Right" "Down" "Left"
[[7]]
[1] "Up" "Right" "Down" "Left"
[[8]]
[1] "Up" "Right" "Down" "Left"
[[9]]
[1] "Up" "Right" "Down" "Left"
[[10]]
[1] "Up" "Right" "Down" "Left"
[[11]]
[1] "None"
[[12]]
[1] "None"\(P(s' | s, a)\) is the probability of going from state \(s\) to \(s'\) by when taking action \(a\). We will create a matrix \(P_a(s' | s)\) for each action.
calc_transition <- function(s, action) {
action <- match.arg(action, choices = A)
if(length(s) > 1) return(t(sapply(s, calc_transition, action = action)))
# deal with absorbing and illegal state
if(s == 11 || s == 12 || s == 5 || action == 'None') {
P <- rep(0, length(S))
P[s] <- 1
return(P)
}
action_to_delta <- list(
'Up' = c(+1, 0),
'Down' = c(-1, 0),
'Right' = c(0, +1),
'Left' = c(0, -1)
)
delta <- action_to_delta[[action]]
dr <- delta[1]
dc <- delta[2]
rc <- s_to_rc(s)
r <- rc[1]
c <- rc[2]
if(dr != 0 && dc != 0)
stop("You can only go up/down or right/left!")
P <- matrix(0, nrow = ROWS, ncol = COLS)
# UP/DOWN
if(dr != 0) {
new_r <- r + dr
if(new_r > ROWS || new_r < 1) new_r <- r
## can't got to (2, 2)
if(new_r == 2 && c == 2) new_r <- r
P[new_r, c] <- .8
if(c < COLS & !(r == 2 & (c + 1) == 2))
P[r, c + 1] <- .1 else P[r, c] <- P[r, c] + .1
if(c > 1 & !(r == 2 & (c - 1) == 2))
P[r, c - 1] <- .1 else P[r, c] <- P[r, c] + .1
}
# RIGHT/LEFT
if(dc != 0) {
new_c <- c + dc
if(new_c > COLS || new_c < 1) new_c <- c
## can't got to (2, 2)
if(r == 2 && new_c == 2) new_c <- c
P[r, new_c] <- .8
if(r < ROWS & !((r + 1) == 2 & c == 2))
P[r + 1, c] <- .1 else P[r, c] <- P[r, c] + .1
if(r > 1 & !((r - 1) == 2 & c == 2))
P[r - 1, c] <- .1 else P[r, c] <- P[r, c] + .1
}
as.vector(P)
}Try to go up from state 1 (this is (1,1), the bottom left corner). Note: we cannot go left so there is a .1 chance to stay in place.
[1] 0.1 0.8 0.0 0.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 col
row 1 2 3 4
3 0.0 0.0 0 0
2 0.8 0.0 0 0
1 0.1 0.1 0 0Try to go right from (2,1). Since right is blocked, there is a .8 probability of staying in place.
Calculate transitions for each state to each other state. Each row represents a state \(s\) and each column a state \(s'\) so we get a complete definition for \(P_a(s' | s)\). Note that the matrix is stochastic (all rows add up to 1).
Create a matrix for each action.
P_matrices <- lapply(A, FUN = function(a) calc_transition(S, a))
names(P_matrices) <- A
str(P_matrices)List of 5
$ Up : num [1:12, 1:12] 0.1 0 0 0.1 0 0 0 0 0 0 ...
$ Right: num [1:12, 1:12] 0.1 0.1 0 0 0 0 0 0 0 0 ...
$ Down : num [1:12, 1:12] 0.9 0.8 0 0.1 0 0 0 0 0 0 ...
$ Left : num [1:12, 1:12] 0.9 0.1 0 0.8 0 0 0 0 0 0 ...
$ None : num [1:12, 1:12] 1 0 0 0 0 0 0 0 0 0 ...Create a function interface for \(P(s' | s, a)\).
\(R(s, a, s')\) define the reward for the transition from \(s\) to \(s'\) with action \(a\).
For the textbook example we have:
Note that once you are in an absorbing state (11 or 12), then the problem is over and there is no more reward!
The solution to an MDP is a policy \(\pi\) which defines which action to take in each state.
It can be shown that all MPDs have an optimal deterministic policy with one action per state. We represent deterministic policies as a vector of actions. I make up a policy that always goes up and then to the right once the agent hits the top.
[1] "Up" "Up" "Right" "Up" "Up" "Right" "Up" "Up" "Right"
[10] "Up" "Up" "Up" col
row 1 2 3 4
3 "Right" "Right" "Right" "Up"
2 "Up" "Up" "Up" "Up"
1 "Up" "Up" "Up" "Up"We can also create a random policy by randomly choosing from the available actions for each state.
create_random_deterministic_policy <-
function()
structure(sapply(
S,
FUN = function(s)
sample(actions(s), 1L)
), names = S)
set.seed(1234)
pi_random <- create_random_deterministic_policy()
pi_random 1 2 3 4 5 6 7 8 9 10
"Left" "Left" "Right" "Right" "None" "Left" "Down" "Up" "Up" "Right"
11 12
"None" "None" col
row 1 2 3 4
3 "Right" "Left" "Up" "None"
2 "Left" "None" "Up" "None"
1 "Left" "Right" "Down" "Right"Stochastic policies use probabilities of actions in each state. This is useful to create policies that explore by trying different actions in the same state.
We use as simple table with probabilities where each row is a state and the columns are the actions. Here we create a random \(\epsilon\)-soft policy. Each available has at least a probability of \(\epsilon\).
We can make a deterministic policy soft.
make_policy_soft <-
function(pi, epsilon = 0.1) {
if(!is.vector(pi))
stop("pi is not a deterministic policy!")
p <-
matrix(0,
nrow = length(S),
ncol = length(A),
dimnames = list(S, A))
for (s in S) {
p[s, actions(s)] <- epsilon / length(actions(s))
p[s, pi[s]] <- p[s, pi[s]] + (1 - epsilon)
}
p
}
make_policy_soft(pi_random) Up Right Down Left None
1 0.025 0.025 0.025 0.925 0
2 0.025 0.025 0.025 0.925 0
3 0.025 0.925 0.025 0.025 0
4 0.025 0.925 0.025 0.025 0
5 0.000 0.000 0.000 0.000 1
6 0.025 0.025 0.025 0.925 0
7 0.025 0.025 0.925 0.025 0
8 0.925 0.025 0.025 0.025 0
9 0.925 0.025 0.025 0.025 0
10 0.025 0.925 0.025 0.025 0
11 0.000 0.000 0.000 0.000 1
12 0.000 0.000 0.000 0.000 1Or we can create a completely random soft policy
create_random_epsilon_soft_policy <-
function(epsilon = 0.1) {
# split total randomly into n numbers that add up to total
random_split <- function(n, total) {
if (n == 1)
return(total)
bordersR <- c(sort(runif(n - 1)), 1)
bordersL <- c(0, bordersR[1:(n - 1)])
(bordersR - bordersL) * total
}
p <-
matrix(0,
nrow = length(S),
ncol = length(A),
dimnames = list(S, A))
for (s in S)
p[s, actions(s)] <- epsilon / length(actions(s)) +
random_split(n = length(actions(s)), 1 - epsilon)
p
}
set.seed(1234)
pi_random_epsilon_soft <- create_random_epsilon_soft_policy()
pi_random_epsilon_soft Up Right Down Left None
1 0.12733307 0.47101419 0.0367222 0.36493054 0
2 0.58604150 0.04023805 0.2235443 0.15017615 0
3 0.03354618 0.22574927 0.4151799 0.32552462 0
4 0.48782603 0.05265132 0.1587548 0.30076784 0
5 0.00000000 0.00000000 0.0000000 0.00000000 1
6 0.27946023 0.03362403 0.5930059 0.09390986 0
7 0.26513870 0.04246225 0.5209651 0.17143393 0
8 0.19305051 0.06595281 0.1009479 0.64004879 0
9 0.06099633 0.13214508 0.1542826 0.65257597 0
10 0.22191959 0.30120821 0.2814109 0.19546130 0
11 0.00000000 0.00000000 0.0000000 0.00000000 1
12 0.00000000 0.00000000 0.0000000 0.00000000 1The value function given a policy \(\pi\) is defined as the expected utility:
\(U^\pi = E\left[\sum_{t=0}^\infty \gamma^t R(s_t, \pi(s_t), s_{t+1})\right]\)
We need to define the discount factor.
We can directly estimate the expected utility of a state using a Monte Carlo simulation that follows the policy. For the stochastic transition model, we need to be able to sample the state \(s'\) the system transitions to when using action \(a\) in state \(s\).
sample_transition <- function(s, a)
sample(S, size = 1, prob = P_matrices[[a]][s,])
sample_transition(1, 'Up')[1] 1
1 2 4
9 84 7 We can now simulate the utility for one episode. Note that we use the cutoff max_t in case a policy does not end up in a terminal state before that.
simulate_utility <- function(pi, s0 = 1, max_t = 100) {
s <- s0
U <- 0
t <- 0
while (TRUE) {
## get action from policy (matrix means it is a stochastic policy)
if (!is.matrix(pi))
a <- pi[s]
else
a <- sample(A, size = 1, prob = pi[s, ])
## sample a transition given the action from the policy
s_prime <- sample_transition(s, a)
##
U <- U + GAMMA ^ t * R(s, a, s_prime)
s <- s_prime
## reached an absorbing state?
if (s == 11 || s == 12 || s == 5)
break
t <- t + 1
if (t >= max_t)
break
}
U
}Simulate the N episodes.
The expected utility for starting from state \(s_0 = 1\) is.
Compare with the utility of the random policy.
utility_random
-4
1000 The random policy performs really poorly. It most likely always stumbles around for max_t moves at a cost of .04 each. The manually created policy should obviously do much better.
We can use simulation to estimate the expected utility for starting from each state following the policy.
It is convenient for solving MDPs using state-action values so we first define the Q-function.
Let’s define the expected utility of a state given that the agent always chooses the optimal action (which it hopefully will if it is rational).
\(U(s) = \max_{a \in A(s)}\sum_{s'} P(s'|s,a) [R(s,a,s') + \gamma U(s')]\)
This equation is called the Bellman equation resulting in an equation system with one equation per state \(s\). This system of equations is hard to solve for all \(U(s)\) values because of the nonlinear \(\max()\) operator.
Lets define a function for the expected utility of any possible action \(a\) (not just the optimal one) in a given state \(s\). This is called the (Q-function or state-action value function):
\(Q(s,a) = \sum_{s'} P(s'|s,a) [R(s,a,s') + \gamma U(s')]\)
This function is convenient for solving MDPs and can easily be implemented.
The issue is that we need to know \(U\) representing the expected utility of a state given optimal decisions is needed. Value iteration uses a simple iterative algorithm to solve this problem by successively updating Q and U.
Note that \(U(s) = \max_a Q(s, a)\) holds and we get:
\(Q(s,a) = \sum_{s'} P(s'|s,a) [R(s,a,s') + \gamma \max_{a'} Q(s', a')]\)
The goal is to find the unique utility function \(U\) (a vector of utilities, one for each state) for the MDP and then derive the implied optimal policy \(\pi^*\).
Algorithm: Start with a \(U(s)\) vector of 0 for all states and then update (Bellman update) the vector iteratively until it converges. This procedure is guaranteed to converge to the unique optimal solution. The pseudocode from AIMA Figure 17.6.:
Stopping criterion: \(||U^\pi - U||_\infty\) is called the policy loss (i.e., the most the agent can loose by using policy \(\pi\) instead of the optimal policy \(\pi^*\) implied in \(U\)). The max-norm \(||x||_\infty\) is defined as the largest component of a vector \(x\).
It can be shown that if \(||U_{i+1} - U_i||_\infty < \epsilon(1-\gamma)/\gamma\) then \(||U_{i+1} - U||_\infty < \epsilon\). This can be used as a stopping criterion with guarantee of a policy loss of less than \(\epsilon\).
value_iteration <- function(eps, verbose = FALSE) {
U_prime <- rep(0, times = length(S))
i <- 1L
while (TRUE) {
if(verbose) cat("Iteration:", i)
#cat("U:", U_prime, "\n")
U <- U_prime
delta <- 0
for (s in S) {
U_prime[s] <- max(sapply(
actions(s),
FUN = function(a)
Q_value(s, a, U)
))
delta <- max(delta, abs(U_prime[s] - U[s]))
}
if(verbose) cat(" -> delta:", delta, "\n")
if (delta <= eps * (1 - GAMMA) / GAMMA)
break
i <- i + 1L
}
cat("Iterations needed:", i, "\n")
U
}Iterations needed: 58 col
row 1 2 3 4
3 0.8515582 0.9078082 0.9578082 0.0000000
2 0.8015582 0.0000000 0.7002740 0.0000000
1 0.7453082 0.6953082 0.6514155 0.4279249For the optimal policy, we choose in each state the action that maximizes the expected utility. This is called the maximum expected utility (MEU) policy. The action that maximizes the utility can be found using the Q-function.
\(\pi^*(s) = \mathrm{argmax}_a Q(s, a)\)
For state 1, 'Up' is the best move
Up Right Down Left None
0.7453082 0.6709332 0.7003082 0.7109332 NA Calculate the Q-function for all \(S \times A\).
Q_value_vec <- Vectorize(Q_value, vectorize.args = c("s", "a"))
QVs <- outer(S, A, FUN = function(s, a) Q_value_vec(s, a, U = U))
colnames(QVs) <- A
QVs Up Right Down Left None
[1,] 0.7453082 0.6709332 0.7003082 0.7109332 NA
[2,] 0.8015582 0.7609332 0.7165582 0.7609332 NA
[3,] 0.8171832 0.8515582 0.7771832 0.8065582 NA
[4,] 0.6559189 0.6201941 0.6559189 0.6953082 NA
[5,] NA NA NA NA 0
[6,] 0.8671832 0.9078082 0.8671832 0.8228082 NA
[7,] 0.6325425 0.4375089 0.5934557 0.6514155 NA
[8,] 0.7002740 -0.6470776 0.4551598 0.6811416 NA
[9,] 0.9210274 0.9578082 0.7150000 0.8520548 NA
[10,] -0.7000660 0.2491324 0.4102740 0.4279249 NA
[11,] NA NA NA NA 0
[12,] NA NA NA NA 0The optimal policy is the greedy policy that always picks the action with the largest Q-value.
col
row 1 2 3 4
3 "Right" "Right" "Right" "None"
2 "Up" "None" "Up" "None"
1 "Up" "Left" "Left" "Left"Estimate the expected utility using simulation.
[1] 0.73692
Compare three policies.
random manual opt
-4.00000 0.69896 0.73692 Since we know that utility_opt is very close to \(U\), we can estimate the policy loss (i.e., the most the agent can loose by using \(\pi\) instead of \(\pi*\)) of the other policies given by:
\(||U^\pi - U||_\infty\)
Here is the policy loss for the manual policy. The maximum norm is the component with the largest difference. First, we calculate the absolute difference for each state.
col
row 1 2 3 4
3 0.008681781 0.01136822 0.005608219 0.000000
2 0.007558219 0.04000000 0.016086027 0.000000
1 0.027828219 0.33630822 0.135135525 1.294925The maximum is:
The policy loss is driven by the bad action taken in state 10 which is at coordinate (1, 4).
Policy iteration tries to directly find the optimal policy. It alternates between two steps:
The pseudocode from AIMA Figure 17.9.:
For policy evaluation, we need to solve:
\(U_i(s) = \sum_{s'} P(s'|s, \pi_i(s))[R(s, \pi_i(s), s') + \gamma U_i(s')]\)
This is slightly simpler than the general Bellman equation, since the action in each state is fixed by the policy and there is no non-linear \(\max()\) operator. For small state spaces this can be solved fast using a LP in \(O(n^3)\).
For large state spaces, we can do approximate policy evaluation by performing only a few iterations of a simplified Bellman update:
\(U_{i+1}(s) \leftarrow \sum_{s'} P(s'|s, \pi_i(s))[R(s, \pi_i(s), s') + \gamma U_i(s')]\)
We implement here approximate policy evaluation with N iterations.
approx_policy_evaluation <- function(pi, U = NULL, N = 10) {
# start with all 0s if no previous U is given
if (is.null(U))
U <- rep(0, times = length(S))
for (i in seq_len(N)) {
for (s in S) {
U[s] = sum(sapply(
S,
FUN = function(s_prime) {
P(s_prime, s, pi[s]) * (R(s, pi[s], s_prime) + GAMMA * U[s_prime])
}
))
}
}
U
} [1] -0.42447359 -0.48544307 -0.62180170 -0.60470718 0.00000000 -0.65671466
[7] -0.64166748 -0.06952390 0.07027771 -0.88579732 0.00000000 0.00000000 [1] 0.6772416 0.7999334 0.8512886 0.1822976 -0.4000000 0.9078008
[7] 0.4536453 0.7002716 0.9578078 -0.8474310 0.0000000 0.0000000We will implement modified policy iteration. Modified means that we use the approximate policy evaluation.
policy_iteration <- function(N = 10) {
U <- rep(0, times = length(S))
pi <- create_random_deterministic_policy()
while (TRUE) {
U <- approx_policy_evaluation(pi, U, N)
unchanged <- TRUE
for (s in S) {
actns <- actions(s)
a <- actns[which.max(sapply(actns, FUN = function(a) Q_value(s, a, U)))]
if (Q_value(s, a, U) > Q_value(s, pi[s], U)) {
pi[s] <- a
unchanged <- FALSE
}
}
if(unchanged) break
}
pi
}Rewriting the Bellman equations as an LP formulation requires replacing the non-linear \(\max()\) operation using additional constraints. The LP can be solved in polynomial time. In practice this is too slow for larger problems. The dynamic programming solution above is typically more efficient, but it is also restricted to small problems.
Reinforcement learning is discussed in AIMA Chapter 22.