First, as an initial value, consider the case where there are multiple options and one of them is randomly selected with the same probability.

def normalize_pi(theta):
    [m, n] = theta.shape
    pi = np.zeros((m, n))
    for i in range(0, m):
        pi[i, :] = theta[i, :] / np.nansum(theta[i, :])
        
    return np.nan_to_num(pi)

In this way, the initial value of pi`` pi_zero is obtained.

pi_zero = normalize_pi(theta_zero)

Let's visualize pi_zero.

visualize_matrix(pd.DataFrame(pi_zero, columns=location[:, 0]))

Note that pi_zero is not a symmetric matrix. For example, there is a side from Sapporo to Aomori and a side from Aomori to Sapporo, but the weights are not the same. This is because, in this graph, there is only one option from Sapporo to the next city, but there are multiple options from Aomori to the next city.

About "behavior" and "state"

Reinforcement learning treats "behavior" and "state" separately. As a result of selecting a certain action in a certain state, it transitions to a different state.

The examples covered in this article are special cases where it is not necessary to distinguish between "behavior" and "state". In other words, "state" refers to the "city at that time", "action" refers to the "next city", and that is the next "state". The code is simpler because you don't have to distinguish between actions and states.

Random walk

First, let's try a trip through Japan with this policy pi still pi_zero. Once you reach a city, it's a random walk that randomly decides which city to go to next.

Probabilistically choose the next action

Create a function that chooses the next city according to the probability indicated by policy pi.

def get_next(town, pi):
    return np.random.choice(range(len(pi)), p=pi[town, :])

For example, according to pi_zero, the next city in Tokyo (No. 12) is calculated as follows.

get_next(12, pi_zero)

13

Of course, the result will change with each run because it is decided randomly.

Random walk search

Now, let's search by random walk.

def explore(pi, start=0, goal=46): #Start is 0 (Sapporo), goal is 46 (Naha)
    route = [start] #List to put travel history
    town = start
    while True:
        next_town = get_next(town, pi) #Determine the next city
        if next_town in route: #If you visit the same city twice
            break #Game over

        route.append(next_town) #Add the next city to the history

        if next_town == goal: #Arrive at the goal
            break #Congratulations
        else: #Not a goal yet
            town = next_town #Set "next city" as "current position"

    return route

Start a random walk according to pi_zero

route = explore(pi_zero)
route

[0, 1, 2, 3, 4, 5]

It's random, so of course the result will change with each run. If it is difficult to imagine just the city number, it will be easier to imagine if you convert it to the city name as follows.

[location[i, 0] for i in route]

['Sapporo', 'Aomori', 'Morioka', 'Sendai', 'Akita', 'Yamagata']

In this example, the game is over because I walked from Sapporo to Yamagata and mistakenly selected the city I visited at the next destination.

By the way, can we really get from Sapporo to Naha with such a random walk?

Evaluation index of search results

Let's create a way to evaluate the quality of the search results. Calculate the straight-line distance from your current location to your destination (Naha) and the total distance traveled from your departure location (Sapporo) to your current location.

def evaluate(route, goal=46):
    dist_goal = dist_mat[route[-1], goal]
    len_route = sum([dist_mat[route[i], route[i + 1]] for i in range(len(route) - 1)])
    return [dist_goal, len_route]

dist_goal is the" straight line distance from your current location to your destination (Naha) ", and len_route is the "total distance traveled from your departure point (Sapporo) to your current location".

Then, it is an execution example.

route = explore(pi_zero)
[location[i, 0] for i in route], evaluate(route)

(['Sapporo', 'Aomori'], [19.597025248636594, 2.3205099949149073])

In this example, the game is over because I chose the route from Sapporo to Aomori and back from Aomori to Sapporo. The distance from the current location (Aomori) to the destination (Naha) is 19.597025248636594, and the distance from the departure point (Sapporo) to the current location (Aomori) is 2.3205099949149073.

The goal of a trip through Japan is to find a route where this dist_goal is zero, and then find the one with the smallest len_route. Shorter len_route is not a good thing.

From now on, I will try the game of moving from Sapporo to Naha many times, but when I get a certain result, let's create a function to judge whether it is the best one so far.

best_dist_goal = 1000000 #Initialize an impossibly large value
best_len_route = 1000000 #Initialize an impossibly large value

def is_best_ever():
    if best_dist_goal >= dist_goal: #Best dist_If equal to or less than the goal value
        if best_dist_goal > dist_goal: #Best dist of the past_If less than the goal value
            return True #Best ever
        elif best_len_route > len_route: #Best len of the past_If less than the route value
            return True #Best ever
        else:
            return False #Not the best
    else:
        return False #Not the best

Search execution

Now that we're ready, we'll repeat the random walk search 50,000 times.

%%time #Measure execution time
pi = pi_zero.copy() #Initialization of pi
theta = theta_zero.copy() #initialization of theta

best_dist_goal = 1000000 #Initialize an impossibly large value
best_len_route = 1000000 #Initialize an impossibly large value

dist_goal_history0 = [] # dist_history of goal
len_route_history0 = [] # len_route history
best_route0 = [] #Best route

for itera in range(50000): #Repeat 50,000 times
    route = explore(pi) #search
    dist_goal, len_route = evaluate(route) #Evaluation
    dist_goal_history0.append(dist_goal) #Record
    len_route_history0.append(len_route) #Record
    
    if is_best_ever(): #If the best of the past
        best_dist_goal = dist_goal #Best dist_goal
        best_len_route = len_route #Best len_route
        best_route0 = route #Best route

CPU times: user 10 s, sys: 118 ms, total: 10.2 s
Wall time: 11.8 s

Search results

Create a function to visualize the result. First, a function that visualizes the distribution of the obtained dist_goal`` len_route

def draw_histgrams(dist_goal_history, len_route_history):
    plt.figure(figsize=(max(len_route_history) / 4, max(dist_goal_history) / 4))
    x_max = max(dist_goal_history + len_route_history)
    ax = plt.subplot(2,1,1)
    ax.hist(dist_goal_history, label='dist_goal', bins=20)
    ax.set_xlim([0, x_max])
    ax.grid()
    ax.legend()
    ax = plt.subplot(2,1,2)
    ax.hist(len_route_history, label='len_route', bins=20)
    ax.set_xlim([0, x_max])
    ax.grid()
    ax.legend()

draw_histgrams(dist_goal_history0, len_route_history0)

Starting from Sapporo, it's hard to reach Naha. If you're lucky, you may get there. You can see that there are many cases where the game is over at the beginning of the search.

Next, let's look at the relationship of dist_goal`` len_route.

def draw_scatter(dist_goal_history, len_route_history):
    plt.figure(figsize=(max(len_route_history) / 4, max(dist_goal_history) / 4))
    plt.scatter(len_route_history, dist_goal_history, alpha=0.5)
    plt.ylabel('dist_goal')
    plt.xlabel('len_route')
    plt.ylim([0, max(dist_goal_history) + 1])
    plt.xlim([0, max(len_route_history) + 1])
    plt.grid()
    plt.show()

draw_scatter(dist_goal_history0, len_route_history0)

You can see that we have not reached Naha, and even if the game is over at the same place, the distance traveled to that point varies (various routes are selected).

Next, let's see how dist_goal`` len_route has changed.

def visualize_history(history, interval=50, window=50):
    plt.plot(range(0, len(history), interval),  #Overall maximum
             [np.array(history)[:i].max() for i in range(len(history)) if (i%interval) == 1], label='max')
    plt.plot(range(window, len(history)+window, interval), 
             [np.array(history)[i:i + window].mean() for i in range(len(history)) if (i%interval) == 1], 
             label='mean(recent)') #Average of recent interval times
    plt.plot(range(0, len(history), interval), 
             [np.array(history)[:i].mean() for i in range(len(history)) if (i%interval) == 1], label='mean') #Overall average
    plt.plot(range(0, len(history), interval), 
             [np.array(history)[:i].min() for i in range(len(history)) if (i%interval) == 1], label='min') #Overall minimum
    plt.legend()
    plt.grid()
    plt.show()

visualize_history(dist_goal_history0)

dist_goal sometimes updates the minimum value, but it does not reach Naha. And no matter how many times the search is repeated, the average value will not improve. It's natural because I haven't learned anything.

visualize_history(len_route_history0)

Similarly, len_route sometimes updates the maximum value, but it does not reach Naha, and the average value does not improve no matter how many times the search is repeated. It's natural because I haven't learned anything.

Finally, let's display the best route found in this random walk.

def draw_route(route):
    plt.figure(figsize=(11, 9))
    for i in route:
        plt.text(location[i, 2], location[i, 1], location[i, 0], alpha=0.4, size=12)
    plt.grid()
    plt.plot([location[i][2] for i in route], [location[i][1] for i in route])
    plt.ylabel('latitude')
    plt.xlabel('longitude')
    plt.show()

draw_route(best_route0)

I tried very hard, but I could only reach Saga. If you are lucky, you may reach Naha. But you don't learn it.

Policy gradient method

Now, it's time to start reinforcement learning </ b>. Reinforcement learning methods can be broadly divided into Policy gradient method </ b> and Value iteration method </ b>. First, let's do the policy gradient method.

def update_pi(pi, route, goal=46, eta=0.1): #Eta (η,eta) is the learning rate
    new_pi = pi.copy() #Copy numpy array
    for i in range(len(route) - 1): #For each side on the route
        town1 = route[i] #City of origin of side i
        town2 = route[i + 1] #City at the end of side i
        new_pi[town1, town2] +=  eta / (dist_mat[route[-1], goal] + 1)
        #The closer the end point of the route is to the goal, the higher the score
    
    return normalize_pi(new_pi) #Update pi

%%time
pi = pi_zero.copy()

best_dist_goal = 1000000
best_len_route = 1000000

dist_goal_history1 = []
len_route_history1 = []
best_route1 = []
for itera in range(50000):
    route = explore(pi)
    dist_goal, len_route = evaluate(route)
    dist_goal_history1.append(dist_goal)
    len_route_history1.append(len_route)
        
    pi = update_pi(pi, route)
    
    if is_best_ever():
        best_dist_goal = dist_goal
        best_len_route = len_route
        best_route1 = route

CPU times: user 59.9 s, sys: 340 ms, total: 1min
Wall time: 1min 1s

draw_histgrams(dist_goal_history1, len_route_history1)

draw_scatter(dist_goal_history1, len_route_history1)

visualize_history(dist_goal_history1)

visualize_history(dist_goal_history1[:5000])

visualize_history(len_route_history1)

visualize_history(len_route_history1[:5000])

draw_route(best_route1)

pi_pg = pi.copy()

visualize_matrix(pd.DataFrame(pi_pg, columns=location[:, 0]))

visualize_matrix(pd.DataFrame(pi_pg - pi_zero, columns=location[:, 0]))

def generate_Q(theta): #Generate unequal random probabilities from the choices you can take
    [m, n] = theta.shape
    Q = np.zeros((m, n))
    rand = np.random.rand(m, n)
    for i in range(m):
        for j in range(n):
            if theta[i, j] == 1:
                Q[i, j] = rand[i, j]
        
    return normalize_pi(Q)

Q_zero = generate_Q(theta_zero)

visualize_matrix(pd.DataFrame(Q_zero, columns=location[:, 0]))

visualize_matrix(pd.DataFrame(Q_zero - pi_zero, columns=location[:, 0]))

def get_next(town, Q, epsilon=0.1):
    if np.random.rand() < epsilon:
        return np.random.choice(range(len(Q)), p=Q[town, :])
    else:
        return np.nanargmax(Q[town, :])

def sarsa(town, Q, prev_t, next_t, reward, eta=0.1, gamma=0.9, goal=46):
    if reward == 1: #dist_mat[town, goal] == 0:
        Q[prev_t, town] = Q[prev_t, town] + eta * (reward - Q[prev_t, town])
    elif reward >= 0:
        Q[prev_t, town] = Q[prev_t, town] + eta * (reward - Q[prev_t, town] + gamma * Q[town, next_t])
    else:
        Q[prev_t, town] = Q[prev_t, town] - eta * Q[prev_t, town]
        
    return normalize_pi(Q)

def explore_epsilon_greedy(Q, epsilon=0.1, eta=0.1, gamma=0.9, start=0, goal=46, min_epsilon=0.01):
    epsilon = max(epsilon, min_epsilon)
    route = [start]
    town = get_next(start, Q, epsilon)
    prev_t = start
    while True:
        if town in route:
            reward = -1
            Q = sarsa(town, Q, prev_t, next_t, reward, eta, gamma)
            break
        elif town == goal:
            reward = 1
            route.append(town)
            Q = sarsa(town, Q, prev_t, next_t, reward, eta, gamma)
            break
        else:
            reward = 0
            route.append(town)
            next_t = get_next(town, Q, epsilon)
            Q = sarsa(town, Q, prev_t, next_t, reward, eta, gamma)
            prev_t = town
            town = next_t
    
    return [route, Q]

%%time
eta = 0.1 #Learning rate
gamma = 0.9 #Time discount rate
epsilon = 0.5
Q = Q_zero.copy()

best_dist_goal = 1000000
best_len_route = 1000000

best_route2 = []
dist_goal_history2 = []
len_route_history2 = []
for itera in range(50000):
    epsilon = epsilon * 0.99
    route, Q = explore_epsilon_greedy(Q, epsilon, eta, gamma)
    dist_goal, len_route = evaluate(route)
    dist_goal_history2.append(dist_goal)
    len_route_history2.append(len_route)

    if is_best_ever():
        best_dist_goal = dist_goal
        best_len_route = len_route
        best_route2 = route

CPU times: user 7min 50s, sys: 948 ms, total: 7min 50s
Wall time: 7min 53s

draw_histgrams(dist_goal_history2, len_route_history2)

draw_scatter(dist_goal_history2, len_route_history2) 

visualize_history(dist_goal_history2)

visualize_history(dist_goal_history2[:5000])

visualize_history(len_route_history2)

visualize_history(len_route_history2[:5000])

draw_route(best_route2)

Q_eg = Q.copy()

visualize_matrix(pd.DataFrame(Q_eg, columns=location[:, 0]))

def Q_learning(town, Q, prev_t, reward, eta=0.1, gamma=0.9, goal=46):
    if reward == 1: #dist_mat[town, goal] == 0:
        Q[prev_t, town] = Q[prev_t, town] + eta * (reward - Q[prev_t, town])
    elif reward >= 0:
        Q[prev_t, town] = Q[prev_t, town] + eta * (reward - Q[prev_t, town] + gamma * np.nanmax(Q[town, :]))
    else:
        Q[prev_t, town] = Q[prev_t, town] - eta * Q[prev_t, town]
        
    return normalize_pi(Q)

def explore_Q_learning(Q, epsilon=0.1, eta=0.1, gamma=0.9, start=0, goal=46):
    prev_t = start
    route = [start]
    town = get_next(start, Q, epsilon)
    while True:
        if town in route:
            reward = -1
            Q = Q_learning(town, Q, prev_t, reward, eta, gamma)
            break
        elif town == goal:
            reward = 1 
            route.append(town)
            Q = Q_learning(town, Q, prev_t, reward, eta, gamma)
            break
        else:
            reward = 0
            dist_goal, len_route = evaluate(route)
            if best_dist_goal > dist_goal:
                reward = 1
            route.append(town)
            next_t = get_next(town, Q, epsilon)
            Q = Q_learning(town, Q, prev_t, reward, eta, gamma)
            prev_t = town
            town = next_t
    
    return [route, Q]

%%time
eta = 0.1 #Learning rate
gamma = 0.9 #Time discount rate
epsilon = 0.5
Q = Q_zero.copy()

best_dist_goal = 1000000
best_len_route = 1000000

best_route3 = []
dist_goal_history3 = []
len_route_history3 = []
for itera in range(50000):
    epsilon = epsilon * 0.99
    route, Q = explore_Q_learning(Q, epsilon, eta, gamma)
    dist_goal, len_route = evaluate(route)
    dist_goal_history3.append(dist_goal)
    len_route_history3.append(len_route)

    if is_best_ever():
        best_dist_goal = dist_goal
        best_len_route = len_route
        best_route3 = route

CPU times: user 9min 50s, sys: 1.41 s, total: 9min 52s
Wall time: 9min 54s

draw_histgrams(dist_goal_history3, len_route_history3)

draw_scatter(dist_goal_history3, len_route_history3) 

visualize_history(dist_goal_history3)

visualize_history(dist_goal_history3[:5000])

visualize_history(len_route_history3)

visualize_history(len_route_history3[:5000])

draw_route(best_route3)

Q_qlearn = Q.copy()

visualize_matrix(pd.DataFrame(Q_qlearn, columns=location[:, 0]))