czxttkl

In this post, I am going to give a practicum walk through on time series analysis. All related code is available in a python notebook.

The data we use is International airline passengers: monthly totals in thousands. which can be downloaded here as csv file (on the left panel, select Export->csv(,)). It is a univariate time series.

From intuition, we can see that there is pattern over here: (1) there exists an upward trend; (2) there exists seasonal fluctuation every once a while. However, this time series is not stationary.  

[3] gives an excellent explanation for the meaning of a stationary time series:

To put it mathematically, a stationary time series is:

1. $latex E(Y_1) = E(Y_2) = E(Y_3) = \cdots = E(Y_t) = \text{constant}$
2. $latex Var(Y_1) = Var(Y_2) = Var(Y_3) = \cdots = Var(Y_t)= \text{constant}$
3. $latex Cov(Y_1, Y_{1+k}) = Cov(Y_2, Y_{2+k}) = Cov(Y_3, Y_{3+k}) = \cdots = \gamma_k$ (depending only on lag $latex k$)

And also compare a stationary time series with white noise:

Let $latex \{\epsilon_t\}$ denote such a series then it has zero mean $latex E(\epsilon_t)=0$, $latex Var(\epsilon_t)=0$ for $latex \forall t$ and $latex E(\epsilon_t \epsilon_s)=0$ for $latex \forall t,s$. 

The typical steps to perform time series analysis are as follows:

1. check whether it is a white noise time series. If it is, just stop, because there is no way that you can predict noise. (In practice, however, it is hard to determine whether a time series is white noise just by observing with eyes. So, my suggestion is to go ahead assuming this is a valid time series, and use cross validation to check if there exists any model that can successfully predict the time series.)
2. if this is a valid but non-stationary time series, make it stationary
3. build predictive models on the stationary time series

Here is my understanding why we need to make a time series stationary before training a model. Suppose your time series varies a lot much in a specific time range. If we build a linear regression model based on the mean square error loss, the errors from that time range will dominate the mean square error and force the weight adjustment more biased to the time range. However, we want to build a predictive model that can predict the time series equally well at any time point. 

We can easily detect whether a time series is stationary using a helper function `test_stationarity`  from [3] and [4]. The usage of `test_stationarity` is that if the ‘Test Statistic’ is greater than the ‘Critical Value’, or ‘p-value’ is well below a convincing threshold (say 0.05),  then the time series is stationary. 

Base on the raw time series, we use test_stationarity and find it is not stationary:

The techniques that you can use to make a time series stationary: (1) logarithmic; (2) first difference; (3) seasonal difference. See [4] and [5] for more techniques to make time series stationary.

ARIMA

Autoregressive integrated moving average (ARIMA) model assumes each time step $latex Y_t$ is a linear function of previous time steps $latex Y_{t-1}, Y_{t-2},\cdots$ and error terms $latex \epsilon_{t}, \epsilon_{t-1}, \epsilon_{t-2}, \cdots$. An ARIMA model(p, q, d) first applies $latex d$-difference on itself to obtain $latex \{Y_{t}\}$ and then assume:

$latex Y_t =\beta_0 + \beta_1 Y_{t-1} + \beta_2 Y_{t-2} + \cdots + \beta_p Y_{t-p} + \epsilon_t+ \Phi_{1} \epsilon_{t-1} + \Phi_{2} \epsilon_{t-2} + \cdots + \Phi_{q} \epsilon_{t-q} &s=2$

Here, the error term $latex \{\epsilon_t\}$ is not the error term from the Autoregressive (AR) model. Instead $latex \{\epsilon_t\}$ are estimated purely from data [11].

To determine the lag parameter $latex p$ in AR model, we can use autocorrelation_plot from pandas.autocorrelation_plot, which is the same thing as Autocorrelation Function (ACF) from statsmodels:

Autocorrelation plots are often used for checking randomness in time series. This is done by computing autocorrelations for data values at varying time lags. If time series is random, such autocorrelations should be near zero for any and all time-lag separations. If time series is non-random then one or more of the autocorrelations will be significantly non-zero. Therefore, autocorrelation can be used as a tool to quickly check whether your time series is suitable to use machine learning models to forecast. After all, if your time series is just white noises, it is meaningless to do time series forecast. The horizontal lines displayed in the plot correspond to 95% and 99% confidence bands. The dashed line is 99% confidence band.

Note that the ACF from statsmodels provides a parameter unbiased. If True, the plot will automatically adjust for the length of the time series. [14]

To determine the lag parameter $latex q$ in MA model, we can use Partial Autocorrelation Function (PACF) in statsmodels. 

The rules to derive right $latex p$ and $latex q$ can be found from [12].

On the other hand, [4] points out that $latex p$ should be the value when the PACF chart crosses the upper confidence interval for the first time. $latex q$ should be the value when the ACF chart crosses the upper confidence interval for the first time. (Also see [15])

In my opinion, deriving $latex p$ and $latex q$ from observing ACF and PACF plots is a very ad-hoc, vague method. I’d prefer to do grid search on parameters as in [13] then filter out the best model according to AIC or similar metrics. 

From [13]: AIC (Akaike Information Criterion) value, which is conveniently returned with ARIMA models fitted using statsmodels. The AIC measures how well a model fits the data while taking into account the overall complexity of the model. A model that fits the data very well while using lots of features will be assigned a larger AIC score than a model that uses fewer features to achieve the same goodness-of-fit. 

Seasonality plot allows you to separate trend, seasonality and residual from a time series. According to [10], a time series consists of four components:

1. level, the average value in the series
2. trend, the increasing or decreasing value in the series
3. seasonality, the repeating short-term cycle in the series
4. noise, the random variation in the series

Below is the seasonality plot based on our data:

The residual is thus the combination of level and noise. One thing that I am not sure is how to use seasonality plot. One possible way is to fit ARIMA on residual because that’s the thing de-trended and de-seasonalized (see one person asked in a comment [8]). However, since the residual is already the level (seen as a constant?) + noise, is there any point to fit ARIMA? Also, even if we can fit an ARIMA on the residual, when we use it to predict the future, how can we transform the predicted value back to the raw scale? [15] actually supports this direction and also mentions the same concern: “Now that you have stationarized your time series, you could go on and model residuals (fit lines between values in the plot). However, as the patterns for the trend and seasonality information extracted from the series that are plotted after decomposition are still not consistent and cannot be scaled back to the original values, you cannot use this approach to create reliable forecasts“. Another direction is to use exogenous variables to fit the residuals, treating them as something not predictable by the time series itself [3]. However, we also face the same issue about how to transfer the predicted residuals back to the raw scale.

There are also other resources introducing seasonality decomposition, such as [9].

Based on what I’ve read so far, I’d say using seasonality plot as an exploration tool. But do not rely on it to train forecast model. After all, modern python ARIMA APIs integrates fitting trend and seasonality at the same time so that you don’t need to separate them out manually [7].

Reference:

[1] http://machinelearningmastery.com/time-series-prediction-with-deep-learning-in-python-with-keras/

[2] http://pandasplotting.blogspot.com/2012/06/autocorrelation-plot.html

[3] http://www.seanabu.com/2016/03/22/time-series-seasonal-ARIMA-model-in-python/

[4] https://www.analyticsvidhya.com/blog/2016/02/time-series-forecasting-codes-python/

[5] http://people.duke.edu/~rnau/whatuse.htm

[6] http://www.statsmodels.org/dev/generated/statsmodels.tsa.arima_model.ARIMA.html

[7] http://www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html

[8] https://stats.stackexchange.com/questions/223117/forecasting-residuals-from-seasonal-decomposition-appears-to-be-highly-auto-cor

[9] http://www.cbcity.de/timeseries-decomposition-in-python-with-statsmodels-and-pandas

[10] http://machinelearningmastery.com/decompose-time-series-data-trend-seasonality/

[11] https://stats.stackexchange.com/questions/26024/moving-average-model-error-terms/74826#74826

[12] https://www.youtube.com/watch?v=Aw77aMLj9uM

[13] https://www.digitalocean.com/community/tutorials/a-guide-to-time-series-forecasting-with-arima-in-python-3

[14] https://stats.stackexchange.com/questions/164526/is-an-auto-correlation-plot-suitable-for-determining-at-what-point-time-series-d

[15] https://datascience.ibm.com/exchange/public/entry/view/815137c868b916821dec777bdc23013c

https://www.google.com/search?q=episode+experience+replay&oq=episode+experience+replay&aqs=chrome..69i57.26352j0j7&sourceid=chrome&ie=UTF-8

https://www.google.com/search?q=experience+replay&oq=experience+rep&aqs=chrome.0.0j69i57j0l4.2129j0j7&sourceid=chrome&ie=UTF-8

In this post, I am going to talk my understanding in tricks in training deep neural network.

ResNet [1]

Why does ResNet network work? https://www.quora.com/How-does-deep-residual-learning-work

Here is my answer:

It is hard to know the desired depth of a deep network. If layers are too deep, errors are hard to propagate back correctly. if layers are too narrow, we may not learn enough representation power.

However, in deep residual network, it is safe to train very deep layers in order to get enough learning power without worrying about the degradation problem too much, because in the worst case, blocks in those “unnecessary layers” can learn to be an identity mapping and do no harm to performance. This is achieved by the solver driving weights of ReLu layers close to zeros thus only the shortcut connection is active and acts as an identity mapping. Though not proved theoretically, adjusting weights close to zeros might be an easier task to the solver than adjusting weights to effective representation all at once.

The authors observe empirically (Fig. 7) that ResNet has smaller magnitude of layer responses on average than plain networks, suggesting that many blocks are just learning little, incremental information.

To conclude, the core idea of ResNet is providing shortcut connection between layers, which make it safe to train very deep network to gain maximal representation power without worrying about the degradation problem, i.e., learning difficulties introduced by deep layers.

All my answer is based on empirical observation and intuition. I’d like to know more theories behind ResNet.

Now according to [2], all the following techniques are regularization methods.

Data Augmentation

augment the training set via domain-specific transformations. For image data, commonly used transformations include random cropping, random perturbation of brightness, saturation, hue and contrast.

Early Stopping

Early stopping was shown to implicitly regularize on some convex learning problems (Yao et al., 2007; Lin et al., 2016)

Dropout 

mask out each element of a layer output randomly with a given dropout probability. [1] says that We adopt batch normalization (BN) [16] right after each convolution and before activation, following [16].  but We
do not use dropout [14], following the practice in [16].

Weight Decay 

Why is it equivalent to $latex \ell_2$ norm regularizer on the weights? Also equivalent to a hard constraint of the weights to an Euclidean ball, with the radius decided by the amount of weight decay?

Batch Norm

https://theneuralperspective.com/2016/10/27/gradient-topics/

http://leix.me/2017/03/02/normalization-in-deep-learning/

ReLU, Leaky ReLU & MaxOut

Comparison between different activation functions

https://en.wikipedia.org/wiki/Activation_function#Comparison_of_activation_functions

Reference

[1] Deep Residual Learning for Image Recognition

[2] Understand Deep Learning Requires Rethinking Generalization

[3] Note for [2]: https://danieltakeshi.github.io/2017/05/19/understanding-deep-learning-requires-rethinking-generalization-my-thoughts-and-notes

[4] https://datascience.stackexchange.com/questions/14349/difference-of-activation-functions-in-neural-networks-in-general

In this post, I am doing a brief code walkthrough for the code written in https://medium.com/emergent-future/simple-reinforcement-learning-with-tensorflow-part-8-asynchronous-actor-critic-agents-a3c-c88f72a5e9f2

The code implements A3C algorithm (Asynchronous Methods for Deep Reinforcement Learning). It follows the pseudocode given in supplemental part in the paper:

The structure of this model is:

For LSTM structure detail, refer to http://colah.github.io/posts/2015-08-Understanding-LSTMs/. I am using the same notation in the model structure diagram as in the link.

Below is the code: 

"""
Asynchronous Advantage Actor-Critic algorithm

See:
https://medium.com/emergent-future/simple-reinforcement-learning-with-tensorflow-part-8-asynchronous-actor-critic-agents-a3c-c88f72a5e9f2
https://arxiv.org/pdf/1602.01783.pdf
"""

import copy
import threading
import multiprocessing
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import tensorflow.contrib.slim as slim
import scipy.signal
from helper import *
from vizdoom import *

from random import choice
from time import sleep
from time import time


# Copies one set of variables to another.
# Used to set worker network parameters to those of global network.
def update_target_graph(from_scope,to_scope):
    from_vars = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES, from_scope)
    to_vars = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES, to_scope)

    op_holder = []
    for from_var,to_var in zip(from_vars,to_vars):
        op_holder.append(to_var.assign(from_var))
    return op_holder

# Processes Doom screen image to produce cropped and resized image.
def process_frame(frame):
    s = frame[10:-10,30:-30]
    s = scipy.misc.imresize(s,[84,84])
    s = np.reshape(s,[np.prod(s.shape)]) / 255.0
    return s

# Discounting function used to calculate discounted returns.
def discount(x, gamma):
    return scipy.signal.lfilter([1], [1, -gamma], x[::-1], axis=0)[::-1]

#Used to initialize weights for policy and value output layers
def normalized_columns_initializer(std=1.0):
    def _initializer(shape, dtype=None, partition_info=None):
        out = np.random.randn(*shape).astype(np.float32)
        out *= std / np.sqrt(np.square(out).sum(axis=0, keepdims=True))
        return tf.constant(out)
    return _initializer


class AC_Network():
    def __init__(self, s_size, a_size, scope, trainer):
        with tf.variable_scope(scope):
            # Input and visual encoding layers
            self.inputs = tf.placeholder(shape=[None, s_size], dtype=tf.float32)
            self.imageIn = tf.reshape(self.inputs, shape=[-1, 84, 84, 1])
            self.conv1 = slim.conv2d(activation_fn=tf.nn.elu,
                                     inputs=self.imageIn, num_outputs=16,
                                     kernel_size=[8, 8], stride=[4, 4], padding='VALID')
            self.conv2 = slim.conv2d(activation_fn=tf.nn.elu,
                                     inputs=self.conv1, num_outputs=32,
                                     kernel_size=[4, 4], stride=[2, 2], padding='VALID')
            hidden = slim.fully_connected(slim.flatten(self.conv2), 256, activation_fn=tf.nn.elu)

            # Recurrent network for temporal dependencies
            lstm_cell = tf.contrib.rnn.BasicLSTMCell(256, state_is_tuple=True)
            self.lstm_cell_c_size = lstm_cell.state_size.c
            self.lstm_cell_h_size = lstm_cell.state_size.h
            c_init = np.zeros((1, lstm_cell.state_size.c), np.float32)
            h_init = np.zeros((1, lstm_cell.state_size.h), np.float32)
            self.state_init = [c_init, h_init]

            c_in = tf.placeholder(tf.float32, [1, lstm_cell.state_size.c])
            h_in = tf.placeholder(tf.float32, [1, lstm_cell.state_size.h])
            self.state_in = (c_in, h_in)
            rnn_in = tf.expand_dims(hidden, [0])
            step_size = tf.shape(self.imageIn)[:1]
            state_in = tf.contrib.rnn.LSTMStateTuple(c_in, h_in)
            lstm_outputs, lstm_state = tf.nn.dynamic_rnn(
                lstm_cell, rnn_in, initial_state=state_in, sequence_length=step_size,
                time_major=False)
            lstm_c, lstm_h = lstm_state
            self.state_out = (lstm_c[:1, :], lstm_h[:1, :])
            rnn_out = tf.reshape(lstm_outputs, [-1, 256])

            # Output layers for policy and value estimations
            self.policy = slim.fully_connected(rnn_out, a_size,
                                               activation_fn=tf.nn.softmax,
                                               weights_initializer=normalized_columns_initializer(0.01),
                                               biases_initializer=None)
            self.value = slim.fully_connected(rnn_out, 1,
                                              activation_fn=None,
                                              weights_initializer=normalized_columns_initializer(1.0),
                                              biases_initializer=None)

            # Only the worker network need ops for loss functions and gradient updating.
            if scope != 'global':
                self.actions = tf.placeholder(shape=[None], dtype=tf.int32)
                self.actions_onehot = tf.one_hot(self.actions, a_size, dtype=tf.float32)
                self.target_v = tf.placeholder(shape=[None], dtype=tf.float32)
                self.advantages = tf.placeholder(shape=[None], dtype=tf.float32)

                self.responsible_outputs = tf.reduce_sum(self.policy * self.actions_onehot, [1])

                # Loss functions
                self.value_loss = 0.5 * tf.reduce_sum(tf.square(self.target_v - tf.reshape(self.value, [-1])))
                self.entropy = - tf.reduce_sum(self.policy * tf.log(self.policy))
                self.policy_loss = -tf.reduce_sum(tf.log(self.responsible_outputs) * self.advantages)
                self.loss = 0.5 * self.value_loss + self.policy_loss - self.entropy * 0.01

                # Get gradients from local network using local losses
                local_vars = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES, scope)
                self.gradients = tf.gradients(self.loss, local_vars)
                self.var_norms = tf.global_norm(local_vars)
                grads, self.grad_norms = tf.clip_by_global_norm(self.gradients, 40.0)

                # Apply local gradients to global network
                global_vars = tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES, 'global')
                self.apply_grads = trainer.apply_gradients(zip(grads, global_vars))

    def state_in_init(self):
        c_init = np.zeros((1, self.lstm_cell_c_size), np.float32)
        h_init = np.zeros((1, self.lstm_cell_h_size), np.float32)
        self.state_init = [c_init, h_init]


class Worker():
    def __init__(self, game, name, s_size, a_size, trainer, model_path, global_episodes):
        self.name = "worker_" + str(name)
        self.number = name
        self.model_path = model_path
        self.trainer = trainer
        self.global_episodes = global_episodes
        self.increment = self.global_episodes.assign_add(1)
        self.episode_rewards = []
        self.episode_lengths = []
        self.episode_mean_values = []
        self.summary_writer = tf.summary.FileWriter("train_" + str(self.number))

        # Create the local copy of the network and the tensorflow op to copy global paramters to local network
        self.local_AC = AC_Network(s_size, a_size, self.name, trainer)
        self.update_local_ops = update_target_graph('global', self.name)

        # The Below code is related to setting up the Doom environment
        game.set_doom_scenario_path("basic.wad")  # This corresponds to the simple task we will pose our agent
        game.set_doom_map("map01")
        game.set_screen_resolution(ScreenResolution.RES_160X120)
        game.set_screen_format(ScreenFormat.GRAY8)
        game.set_render_hud(False)
        game.set_render_crosshair(False)
        game.set_render_weapon(True)
        game.set_render_decals(False)
        game.set_render_particles(False)
        game.add_available_button(Button.MOVE_LEFT)
        game.add_available_button(Button.MOVE_RIGHT)
        game.add_available_button(Button.ATTACK)
        game.add_available_game_variable(GameVariable.AMMO2)
        game.add_available_game_variable(GameVariable.POSITION_X)
        game.add_available_game_variable(GameVariable.POSITION_Y)
        game.set_episode_timeout(300)
        game.set_episode_start_time(10)
        game.set_window_visible(False)
        game.set_sound_enabled(False)
        game.set_living_reward(-1)
        game.set_mode(Mode.PLAYER)
        game.init()
        self.actions = np.identity(a_size, dtype=bool).tolist()
        # End Doom set-up
        self.env = game

    def train(self, rollout, sess, gamma, bootstrap_value):
        rollout = np.array(rollout)
        observations = rollout[:, 0]
        actions = rollout[:, 1]
        rewards = rollout[:, 2]
        next_observations = rollout[:, 3]
        values = rollout[:, 5]

        # Here we take the rewards and values from the rollout, and use them to
        # generate the advantage and discounted returns.
        # The advantage function uses "Generalized Advantage Estimation"
        self.rewards_plus = np.asarray(rewards.tolist() + [bootstrap_value])
        discounted_rewards = discount(self.rewards_plus, gamma)[:-1]
        self.value_plus = np.asarray(values.tolist() + [bootstrap_value])
        advantages = rewards + gamma * self.value_plus[1:] - self.value_plus[:-1]
        advantages = discount(advantages, gamma)

        # Update the global network using gradients from loss
        # Generate network statistics to periodically save
        rnn_state = self.local_AC.state_init
        feed_dict = {self.local_AC.target_v: discounted_rewards,
                     self.local_AC.inputs: np.vstack(observations),
                     self.local_AC.actions: actions,
                     self.local_AC.advantages: advantages,
                     self.local_AC.state_in[0]: rnn_state[0],
                     self.local_AC.state_in[1]: rnn_state[1]}
        v_l, p_l, e_l, g_n, v_n, _ = sess.run([self.local_AC.value_loss,
                                               self.local_AC.policy_loss,
                                               self.local_AC.entropy,
                                               self.local_AC.grad_norms,
                                               self.local_AC.var_norms,
                                               self.local_AC.apply_grads],
                                              feed_dict=feed_dict)
        return v_l / len(rollout), p_l / len(rollout), e_l / len(rollout), g_n, v_n

    def work(self, max_episode_length, gamma, sess, coord, saver):
        episode_count = sess.run(self.global_episodes)
        total_steps = 0
        print("Starting worker " + str(self.number))
        with sess.as_default(), sess.graph.as_default():
            while not coord.should_stop():
                # copy global network parameter to self parameter
                sess.run(self.update_local_ops)
                # when a new episode starts, C_0 & h_0 of LSTM are reset to zero
                self.local_AC.state_in_init()
                episode_buffer = []
                episode_values = []
                episode_frames = []
                episode_reward = 0
                episode_step_count = 0
                d = False

                self.env.new_episode()
                s = self.env.get_state().screen_buffer
                episode_frames.append(s)
                s = process_frame(s)
                rnn_state = self.local_AC.state_init

                while self.env.is_episode_finished() == False:
                    # Take an action using probabilities from policy network output.
                    # after this step, a_dist shape (1,3), v shape (1,1),
                    # rnn_state: first array (C): (1, 256), second array (h): (1, 256)
                    a_dist, v, rnn_state = sess.run(
                        [self.local_AC.policy, self.local_AC.value, self.local_AC.state_out],
                        feed_dict={self.local_AC.inputs: [s],
                                   self.local_AC.state_in[0]: rnn_state[0],
                                   self.local_AC.state_in[1]: rnn_state[1]})
                    a = np.random.choice(a_dist[0], p=a_dist[0])
                    a = np.argmax(a_dist == a)     # return the best action index

                    # see http://vizdoom.cs.put.edu.pl/tutorial and
                    # https://github.com/mwydmuch/ViZDoom/blob/3bdb16935668aa42bb14cc38ac397b8954999cca/doc/DoomGame.md
                    # for reward description
                    # the agent gets rewards for each action which is -1 for each time tick,
                    # -6 if he shots but misses, and 100 if he kills the monster
                    r = self.env.make_action(self.actions[a]) / 100.0     # make_action returns reward

                    # In this example, the episode finishes after 300 tics or when the monster gets killed
                    # ref: http://www.cs.put.poznan.pl/visualdoomai/tutorial.html#basic - Game Runtime
                    d = self.env.is_episode_finished()
                    if d == False:
                        s1 = self.env.get_state().screen_buffer
                        episode_frames.append(s1)
                        s1 = process_frame(s1)
                    else:
                        s1 = s

                    episode_buffer.append([s, a, r, s1, d, v[0, 0]])
                    episode_values.append(v[0, 0])

                    episode_reward += r
                    s = s1
                    total_steps += 1
                    episode_step_count += 1

                    # If the episode hasn't ended, but the experience buffer is full, then we
                    # make an update step using that experience rollout.
                    if len(episode_buffer) == 30 and d != True and episode_step_count != max_episode_length - 1:
                        # Since we don't know what the true final return is, we "bootstrap" from our current
                        # value estimation.
                        v1 = sess.run(self.local_AC.value,
                                      feed_dict={self.local_AC.inputs: [s],
                                                 self.local_AC.state_in[0]: rnn_state[0],
                                                 self.local_AC.state_in[1]: rnn_state[1]})
                        v1 = v1[0, 0]
                        v_l, p_l, e_l, g_n, v_n = self.train(episode_buffer, sess, gamma, v1)
                        episode_buffer = []
                        sess.run(self.update_local_ops)
                        # original code does not update state_init:
                        # in train function, rnn_state is always set to self.state_init which are two zero numpy arrays
                        self.local_AC.state_init = copy.deepcopy(rnn_state)
                    if d == True:
                        break

                self.episode_rewards.append(episode_reward)
                self.episode_lengths.append(episode_step_count)
                self.episode_mean_values.append(np.mean(episode_values))

                # Update the network using the experience buffer at the end of the episode.
                if len(episode_buffer) != 0:
                    v_l, p_l, e_l, g_n, v_n = self.train(episode_buffer, sess, gamma, 0.0)

                # Periodically save gifs of episodes, model parameters, and summary statistics.
                if episode_count % 5 == 0 and episode_count != 0:
                    if self.name == 'worker_0' and episode_count % 25 == 0:
                        time_per_step = 0.05
                        images = np.array(episode_frames)
                        make_gif(images, './frames/image' + str(episode_count) + '.gif',
                                 duration=len(images) * time_per_step, true_image=True, salience=False)
                    if episode_count % 250 == 0 and self.name == 'worker_0':
                        saver.save(sess, self.model_path + '/model-' + str(episode_count) + '.cptk')
                        print("Saved Model")

                    mean_reward = np.mean(self.episode_rewards[-5:])
                    mean_length = np.mean(self.episode_lengths[-5:])
                    mean_value = np.mean(self.episode_mean_values[-5:])
                    summary = tf.Summary()
                    summary.value.add(tag='Perf/Reward', simple_value=float(mean_reward))
                    summary.value.add(tag='Perf/Length', simple_value=float(mean_length))
                    summary.value.add(tag='Perf/Value', simple_value=float(mean_value))
                    summary.value.add(tag='Losses/Value Loss', simple_value=float(v_l))
                    summary.value.add(tag='Losses/Policy Loss', simple_value=float(p_l))
                    summary.value.add(tag='Losses/Entropy', simple_value=float(e_l))
                    summary.value.add(tag='Losses/Grad Norm', simple_value=float(g_n))
                    summary.value.add(tag='Losses/Var Norm', simple_value=float(v_n))
                    self.summary_writer.add_summary(summary, episode_count)

                    self.summary_writer.flush()
                if self.name == 'worker_0':
                    sess.run(self.increment)
                episode_count += 1




max_episode_length = 300
gamma = .99 # discount rate for advantage estimation and reward discounting
s_size = 7056 # Observations are greyscale frames of 84 * 84 * 1
a_size = 3 # Agent can move Left, Right, or Fire
load_model = False
model_path = './model'

tf.reset_default_graph()

if not os.path.exists(model_path):
    os.makedirs(model_path)

# Create a directory to save episode playback gifs to
if not os.path.exists('./frames'):
    os.makedirs('./frames')

with tf.device("/cpu:0"):
    global_episodes = tf.Variable(0, dtype=tf.int32, name='global_episodes', trainable=False)
    trainer = tf.train.AdamOptimizer(learning_rate=1e-4)
    master_network = AC_Network(s_size, a_size, 'global', None)  # Generate global network
    num_workers = multiprocessing.cpu_count()  # Set workers ot number of available CPU threads
    workers = []
    # Create worker classes
    for i in range(num_workers):
        workers.append(Worker(DoomGame(), i, s_size, a_size, trainer, model_path, global_episodes))
    saver = tf.train.Saver(max_to_keep=5)

with tf.Session() as sess:
    coord = tf.train.Coordinator()
    if load_model == True:
        print('Loading Model...')
        ckpt = tf.train.get_checkpoint_state(model_path)
        saver.restore(sess, ckpt.model_checkpoint_path)
    else:
        sess.run(tf.global_variables_initializer())

    # This is where the asynchronous magic happens.
    # Start the "work" process for each worker in a separate threat.
    worker_threads = []
    for worker in workers:
        worker_work = lambda: worker.work(max_episode_length, gamma, sess, coord, saver)
        t = threading.Thread(target=(worker_work))
        t.start()
        sleep(0.5)
        worker_threads.append(t)
    coord.join(worker_threads)

Note that Worker.train implements the loop $latex for \; i \in \{t-1, \cdots,t_{start}\}$ in the pseudocode.

Reinforcement learning algorithms can be divided into many families. In model-free temporal difference methods like Q-learning/SARSA, we try to learn action value for any state-action pair, either by recording (“memorizing”) exact values in a tabular or learning a function to approximate it. Under -greedy, the action to be selected at a state will therefore be  but there is also a small constant chance to be any non-optimal action.

Another family is called policy gradient methods which directly map states to actions. To select actions, they do not need to consult a value table or a value function. Instead, each action can be selected with a probability determined by a parameterized policy function , where is the policy function’s parameters.

The advantages of policy gradient methods over Q-learning/SARSA using greedy are mainly two:

1. in some situations the optimal approximate policy must be stochastic.  An example from [1]: in card games with imperfect information the optimal play is often to do two different things with specific probabilities, such as when bluffing in Poker. Action-value methods have no natural way of finding stochastic optimal policies.
2. problems vary in the complexity of their policies and action-value functions. In some problems, the policy is a much simpler function to approximate than action-values so it will be faster to learn

The general update form of policy gradient methods is , where is performance measure with respect to the policy weights.

Now, the policy gradient theorem can be briefly described as [5]:

In episodic case, the policy gradient theorem derives as follows (according to [1]):

The last expression is an exact formula of the gradient of . It can also be seen as an expectation of over the probability distribution of landing on in steps. Please note for any fixed , . Therefore, we can rewrite the above last expression as an expectation form:

where . is an unbiased estimator of in the last two steps since we do not have estimation for . We use to replace , meaning that the sequence are generated following the policy and the transition probability . Sometimes we can also write because the policy is parameterized by . We can also write as . In other words,

What we do in reality is to use these collected state, actions and rewards as samples to approximate the exact expectation of :

This kind of policy gradient method is called REINFORCE, by Ronald J. Williams from Northeastern University. The original paper [2] is hard to read in my opinion. It directly tells you what is the update rule of by construction, and then tells you the expected update aligns in the same direction as the performance gradient. What I wish to be told is how he derived the update rule of in the first place.

(Updated 09/18/2017: The same derivative procedure of REINFORCE is illustrated more clearly in [10])

(Updated 01/31/2020: The derivative procedure of REINFORCE in continuous state/action space is illustrated in [15])

One important extension of REINFORCE is to offsetting by a baseline , a function of . Intuitively, we need not how good an action is, but how much better this action is compared to average actions. For example, if is uniformly larger than zero for either good or bad actions, will always get updated to encourage either kind of actions. We need to calibrate such that it can differentiate good or bad actions better. Mathematically, adding a baseline can reduce the variance of while still keeping it as unbiased estimator.

First, let’s look at why offsetting by still makes an unbiased estimator (reference from [8]):

To ensure is an unbiased estimate of , should only be a function of , but not . Otherwise is not zero.

It is less obvious that adding can reduce the variance of .

From here, we can see if has large enough covariance with to outweigh its own variance, then the variance is reduced. Unrealistically, if , then variance will be zero, although this is impossible because is only a function of without magic forecast ability to 100% approach .

(sidenote: I cannot follow [8]’s deduction on variance reduction part.)

One way is to train a predictor on with parameter as a baseline:

REINFORCE is an on-policy algorithm because in the gradient update depends on , the returns generated by following the current policy . The specific one-step actor-critic algorithm we just described is also an on-policy algorithm because depends on the next state which is the result of applying at the current state .  There also exists off-policy actor-critics, see an overview of on-policy and off-policy policy gradient methods at [14].

A more recent advance in baseline is Generalized Advantage Estimation (GAE) [6]. They introduce two parameters, and , in an un-discounted reward problem to help estimate with little introduced bias and reduced variance. (Note that how discounted reward problems can be transformed into an un-discounted problem: “But the discounted problem (maximizing can be handled as an instance of the undiscounted problem in which we absorb the discount factor into the reward function, making it time-dependent.”)

They introduce their notations:

Note that is a biased estimator of but as they claim previous works have studied to “reduce variance by downweighting rewards corresponding to delayed effects, at the cost of introducing bias”.

The paper’s goal is to find a good estimator of which is called .

In other word, if is -just, then it helps to construct an unbiased estimator of . Equation (8) just uses the property that the expectation of a sum equals to the sum of expectations.

Now, what other property does  have? If we know any property of , it will help us find more easily. The paper proposes one property:

Sketch proof of proposition 1:

First of all, to understand the notations clearly, think that and are functions with input as the whole trajectory . Similarly, think as a function with the input as the former part of the trajectory .

Now, suppose a satisfied such that where for all , then:

Next, they find a good candidate function for , which is . Only when we know is -just. Otherwise is a biased estimator of . However, if we can take average of different variations of (equation 11, 12, 13, 14, 15, and 16), then we might get a low bias, low variance estimator, which is called .

is -just regardless of the accuracy of (again, this is because ). However is believed (I don’t know how to prove that) to have high variance due to the long sum of rewards. On the other extreme end, has low variance but since we are estimating the value function ,  must be a biased estimator of . with would make a trade-off between variance and bias.

Update 2018-11-08

Policy gradient is better illustrated in several recent posts: [11] and [12]

Reference

[1] Reinforcement learning: An introduction

[2] Simple Statistical Gradient-Following Algorithms for Connectionist Reinforcement Learning

[3] https://medium.com/emergent-future/simple-reinforcement-learning-with-tensorflow-part-8-asynchronous-actor-critic-agents-a3c-c88f72a5e9f2

[4] Asynchronous Methods for Deep Reinforcement Learning

[5] http://www.breloff.com/DeepRL-OnlineGAE/

[6] HIGH-DIMENSIONAL CONTINUOUS CONTROL USING GENERALIZED ADVANTAGE ESTIMATION

[7] Variance Reduction Techniques for Gradient Estimates in Reinforcement Learning

[8] https://danieltakeshi.github.io/2017/03/28/going-deeper-into-reinforcement-learning-fundamentals-of-policy-gradients/

[9] https://danieltakeshi.github.io/2017/04/02/notes-on-the-generalized-advantage-estimation-paper/

[10] https://www.youtube.com/watch?v=kUiR0RLmGCo

[11] https://lilianweng.github.io/lil-log/2018/04/08/policy-gradient-algorithms.html#soft-actor-critic

[12] https://spinningup.openai.com/en/latest/spinningup/rl_intro3.html#expected-grad-log-prob-lemma

[13] Supplement material of DeepMimic: Example-Guided Deep Reinforcement Learning of Physics-Based Character Skills

[14] Unifying On-Policy and Off-Policy Learning

[15] http://web.stanford.edu/class/cme241/lecture_slides/PolicyGradient.pdf

latex for policy gradient theorem:

latex for expectation form rewritten:

latex for adding baseline is still unbiased estimator:

latex for sketch proof of proposition 1:

Why TD($latex lambda$)?

Why actor-critic?

Why eligibility trace?

Why contextual regret minimization?

Importance sampling is a way to reduce variance of your estimation on integration over a region for an integrand.

Let’s first see how traditional Monte Carlo method is used to estimate integration [2]. To estimate $latex \int_a^b f(x) dx$, one can think of reshaping the area to be integrated as a rectangle, whose width is $latex b-a$ and height is $latex E_{\mathcal{U}[a,b]}(f(x))$, a virtual, expected height.

Therefore, to use Monte Carlo estimation, we can randomly draw $latex N$ numbers from the uniform distribution $latex \mathcal{U}[a,b]$, and use the following result as the approximation to $latex \int_a^b f(x) dx$:

$latex \frac{(b-a)}{N}\sum\limits_{i=1}^N f(x_i) &s=2$

Of course, this estimation will vary from run to run (every run is an estimation based on $latex N$ samples). Therefore, Monte Carlo estimation for integration has variance. To reduce this variance, we use importance sampling [1]. The idea is to convert the integration formula into an expectation formula with an auxiliary, known distribution $latex h(x)$:

$latex \int_a^b f(x) d(x) = \int_a^b f(x) \frac{h(x)}{h(x)} dx = \int_a^b \frac{f(x)}{h(x)} h(x) dx = \mathbb{E}_h \big(\frac{f(X)}{h(X)} \big) &s=2$

Using Monte Carlo method to estimate this expectation, we have:

$latex \mathbb{E}_h \big(\frac{f(X)}{h(X)} \big) \approx \frac{1}{n} \sum\limits_{i=1}^n \frac{f(X_i)}{h(X_i)} ,\quad X_i \sim h(x) &s=2$ 

When $latex \alpha h(x) = f(x)$  (assuming $latex f(x)>0$), the importance sampling variance will become zero because $latex \mathbb{E}_h \big(\frac{f(X)}{h(X)}\big) = \alpha$ which is constant no matter of $latex x_i$. Though a nice property, it is impossible to achieve zero variance in reality. If you know $latex h(x)$ is proportional $latex f(x)$, then you know the density of $latex h(x)$, therefore you induce the density of $latex f(x)$ just timing $latex h(x)$ by $latex \alpha$, therefore you do not need Monte Carlo method in the first place to calculate the integration. The common solution is to pick a known distribution $latex h(x)$ that appears close, in terms of the shape, to $latex f(x)$. The closer $latex h(x)$ is to $latex f(x)$, the smaller the variance is. If $latex h(x)$ is a uniform distribution in the range of $latex (a,b)$, it is equivalent to the traditional method of Monte Carlo method we introduced at the beginning of this article.

Here is the code to verify that importance sampling really reduces the variance for estimating integration. 

# estimate integration over a normal distribution (mean=0, std=1) in the range of (a, b)
from scipy.stats import norm
import numpy
from scipy import stats

T = 200  # num of rounds
N = 10   # num of samples per round
a, b = -2., 2.   # lower bound, upper bound
fx_loc, fx_scale = 0, 1   # f(x) is a normal distribution (mean=0, std=1)

class hx1_dist_gen(stats.rv_continuous):
    def _pdf(self, x):
        return 1. / b - 1. / (b**2) * numpy.abs(x)

class hx2_dist_gen(stats.rv_continuous):
    def _pdf(self, x):
        return 1. / (b**2) * numpy.abs(x)

hx1_dist = hx1_dist_gen(a=a, b=b)
hx2_dist = hx2_dist_gen(a=a, b=b)

# monte carlo estimation
mc_res = []
for t in range(T):
    print("MC round %d" % t)
    x = numpy.random.uniform(low=a, high=b, size=N)
    fxs = norm.pdf(x, loc=fx_loc, scale=fx_scale)
    mc_res.append((b - a) * numpy.mean(fxs))

# importance sampling estimation using h1(x)
is1_res = []
for t in range(T):
    print("Importance Sampling h1(x) round %d" % t)
    fx_div_hxs = [] 
    # approximate to let x be drawn in the range (-1, 1)
    x = hx1_dist.rvs(size=N)
    fx = norm.pdf(x, loc=fx_loc, scale=fx_scale)
    hx = hx1_dist.pdf(x)
    fx_div_hxs = fx / hx
    is1_res.append(numpy.mean(fx_div_hxs))

# importance sampling estimation using h2(x)
is2_res = []
for t in range(T):
    print("Importance Sampling h2(x) round %d" % t)
    fx_div_hxs = [] 
    # approximate to let x be drawn in the range (-1, 1)
    x = hx2_dist.rvs(size=N)
    fx = norm.pdf(x, loc=fx_loc, scale=fx_scale)
    hx = hx2_dist.pdf(x)
    fx_div_hxs = fx / hx
    is2_res.append(numpy.mean(fx_div_hxs))

print("exact: {0}".format(norm.cdf(b) - norm.cdf(a)))
print("MC estimation (mean/std): {0}/{1}".format(numpy.mean(mc_res), numpy.std(mc_res)))
print("Importance Sampling h1 estimation (mean/std): {0}/{1}".format(numpy.mean(is1_res), numpy.std(is1_res)))
print("Importance Sampling h2 estimation (mean/std): {0}/{1}".format(numpy.mean(is2_res), numpy.std(is2_res)))

Output on my machine:

exact: 0.9544997361036416
MC estimation (mean/std): 0.958336124778972/0.15458260669929705
Importance Sampling h1 estimation (mean/std): 0.9515730240850715/0.07416369320580925
Importance Sampling h2 estimation (mean/std): 1.0496523604436687/0.7489580081170375

I set $latex f(x)$ as a normal distribution with mean=0 and std=1. I use the traditional Monte Carlo method as well as two kinds of distributions, $latex h_1(x)$ and $latex h_2(x)$, to do importance sampling. All distributions can be represented in the following chart:

From the output we can see that, estimation using importance sampling with $latex h_1(x)$ has the smallest std. This is because $latex h_1(x)$ (the blue triangle) looks the closest to $latex f(x)$. $latex h_2(x)$ (the red shape formed by two triangles) has the largest std because it looks the least close to $latex f(x)$. The traditional MC method （equivalently, importance sampling with a uniform distribution) has a medium std.

Reference

[1] https://www.dropbox.com/s/75vb40b2vfnnlcm/Monte%20Carlo%20Methods%20and%20Importance%20Sampling%5Bai%20classicML%5D.pdf?dl=0

[2] http://ta.twi.tudelft.nl/mf/users/oosterle/oosterlee/lec8-hit-2009.pdf

In my rough understanding, inverse reinforcement learning is a branch of RL research in which people try to perform state-action sequences resembling given tutor sequences.

There are two famous works on inverse reinforcement learning. One is Apprenticeship Learning via Inverse Reinforcement Learning [1], and the other is Maximum Margin Planning [2].

Maximum Margin Planning

In MMP learn a cost function for estimating the cost at each state-action pair. The desired cost function should emit low costs to those state-action pairs appearing in “tutor” sequences (examples) and return high costs to those deviating tutor sequences. If we manage to recover a desired cost function, we can assign the cost for any state-action pair and use path-finding algorithms (e.g., A*, Dijkstra, etc.) to find a sequence on a new environment which has the minimum accumulative sum of costs, that is, a sequence deemed resembling given tutor sequences.

The baisc idea is given in [4]:

As you see, in line 6, you need to adjust cost function to increase on the current “not-perfect” minimum cost path and decrease the example path. However, you can see that the algorithm doesn’t know how much it should increase the cost for minimum cost path. If one path just deviates from the example path a little bit and another path deviates far more abnormally, how should the adjustment of the cost function reflect that?

So the author proposes an improvement called “loss augmentation”. Suppose the original cost at a state-action is , then we say its augmented cost , where is a loss function with large values for those deviated state-action pairs and with close-to-zero values for those resembling tutor sequences. Why do we want the augmented cost to be even smaller for those defiant state-action pairs? Because in this way, we make those deviated pairs even harder to be distinguished with the desired pairs. Therefore, we force the algorithm to try hard to update the cost function until it can differentiate augmented costs from different state-action pairs. In other words, the algorithm will learn the cost function such that it differentiates any a desired state-action pair with a less desired state-action pair  by at least a margin .

Here is the visualization showing how deviated paths get high rewards (+1) and desired paths get low rewards (-1).

And here is the revised algorithm based on Algorithm 1:

Now, let’s look at how to learn the cost function using a framework called maximum margin structured classification.

Before that, we introduce preliminaries. In this paper, we are given examples, each example is a MDP and a tutor sequence. Different examples may come from different MDPs so this model can learn the cost function generalized to new MDPs which are never seen in the training examples. However, for simplicity, we just assume all training examples are different tutor sequences from the same MDP with state and . For simplicity purpose also, an important assumption made throughout [4] is that a tutor sequence is acyclic. This means, any state-action pair is visited at most once in a tutor sequence. Then, it is easy to represent a tutor sequence by , a state-action frequency count vector. If the assumption about tutor sequence being acyclic were not made, which satisfies certain flow constraints.

Now remember that the idea is to approximate everything linearly, such as the loss function or the cost function. Each training example has a loss function, which quantifies the difference between an arbitrary sequence and -th example tutor sequence . Since we want loss function to be a linear function of , we let . The more deviates from , the higher should be. A straightforward way to define  is to make , with 0 denoting a particular state-action is visited by both or neither of and and 1 otherwise. Of course more advanced loss function can be used but the bottom line is to assume . Now as for cost function, cost function should also factor over states and actions. This means, the cost of a sequence should be the sum of the costs at each state-action pairs visited in the sequence. Suppose we have a feature matrix , whose columns are feature vectors at every state-action pair. We also suppose the cost function has a linear weight . Then, the cost at one state-action pair is just where is the column corresponding to the state-action pair.  Finally, the cost of the whole sequence is . You can think as the accumulated sum of feature vectors on visited state-action pairs.

Since we assume all training examples come from the same MDP, we can write the data set as .

Now, the problem is how to obtain . Recall that we want to be able to differentiate desired paths ( from ) from tentative paths (arbitrary , where is all possible sequence in the MDP) by a margin . Therefore, the constraints we would like to impose are: for each , for each , . This actually introduces a large number of constraints when space is large. Therefore, we reduce the number of constraints by rewriting as: for each , . Intuitively, this means if can successfully distinguish from the loss-augmented cost of  the hardest tentative path , then it can distinguish from any other tentative paths. Also, we would like to regularize on the norm of because if were without norm regularization, it can always scale arbitrarily large to  overcome margin difference constraints, voiding the original purpose of adding margins. Moreover, we want to introduce a slack term for each training example. This allows some margin difference constraints to not be satisfied because in real data you may not find a which perfectly satisfies all margin difference constraints.

Putting all we describe above, we come to the optimization objective:

We note that resides in both the objective function and the constraints. We want as small as possible. Therefore, at the minimizer the slack variables will always exactly equal the constraint violation: . Therefore, we can rewrite our constrained objective function into another non-constrained objective function:

Now, is a non-differentiable (because of ), convex objective. The author in [4] proposes to optimize using subgradient method.

In my mind, subgradient method is really similar to normal gradient descend method. The gradient of a differentiable convex function at point satisfies: . This means the growth of line from point is always equal or slower than [7]:

Subgradient equal to gradient at the part being differentiable. At the non-differentiable part, subgradient can be any values such that :

Now let’s look at the part  in . This part can be rewritten as . This can be seen as a max over affine functions , each shown as dashed line the figure below:

The bold line is . The red line is the subgradient at specific value of marked by the vertical dashed line. Therefore, the subgradient of is:

where .

Now, can be updated online, i.e., being updated after seeing example one by one:

where is projection operator that maps to feasible region. For example, sometimes we even require all costs should be positive. How to do projection is something I don’t fully understand in [4]. I hope I will illustrate once I am crystal clear about that.

Reference

[1] Apprenticeship Learning via Inverse Reinforcement Learning

[2] Maximum Margin Planning

[3] Boosting Structured Prediction for Imitation Learning

[4] Learning to Search: Structured Prediction Techniques for Imitation Learning (Ratliff’s thesis)

[5] Max Margin Learning (Youtube video)

[6] Max Margin Planning Slides

[7] Subgradient slides

2018.4.8

Seems like this is a good tutorial (https://thinkingwires.com/posts/2018-02-13-irl-tutorial-1.html) for explaining [8]

[8] Algorithms for Inverse Reinforcement Learning