\documentclass[11pt]{article}
\usepackage{latexsym}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{epsfig}
\usepackage[tight]{subfigure}
\usepackage{comment}

\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{
      \hbox to 5.78in { {#1} \hfill #2 }
      \vspace{4mm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{2mm}
      \hbox to 5.78in { {\em #3 \hfill #4} }
    }
  }
  \end{center}
  \vspace*{4mm}
}

\newcommand{\lecture}[5]{\handout{#1}{#2}{#3}{#4}{#5}}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{fact}[theorem]{Fact}
\newtheorem{assumption}[theorem]{Assumption}

% 1-inch margins, from fullpage.sty by H.Partl, Version 2, Dec. 15, 1988.
\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 8.9in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in

\parindent 0in
\parskip 1.5ex

%\renewcommand{\baselinestretch}{1.25}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Document header: change the lecture number, date, scribe name, and
% lecture topic.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\lecture{Statistical Techniques in Robotics (16-831, F12)}
{Lecture\#20 (Monday November 12)}
{Lecturer: Drew Bagnell}
{Scribe: Venkatraman Narayanan\footnotemark}
{Gaussian Processes}

% Be sure to acknowledge all previous scribes.
\footnotetext{Some content adapted from previous scribes: Yamuna Krishnamurthy, Stephane Ross, and from \emph{Gaussian Processes for Machine Learning} by Carl Edward Rasmussen and Christopher K. I. Williams}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Document body: Scribe your notes here.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Applications of Gaussian Processes}

\begin{figure}[h]
\subfigure[Inverse Kinematics of a Robot Arm~\cite{JanPeters}]{
\includegraphics[width=0.5\textwidth]{Picture2.png}
\label{subfig:wam}
}
\subfigure[Predictive Soil Modeling~\cite{JuanPablo}]{
\includegraphics[width=0.5\textwidth]{Picture1.png}
\label{subfig:soil_model}
}\caption{Applications of Gaussian Processes \label{fig:examples}}
\end{figure}

Gaussian processes can be used as a supervised learning technique for classification as well as regression. An example of a classification task would be to recognize handwritten digits, whereas an example of a regression problem would be to learn the inverse 
dynamics of a robot arm (Fig.~\ref{subfig:wam}). For the latter, the task is to obtain a mapping from the state of the arm (given by the positions, velocities and accelerations of the joints) to the corresponding torques on the joints. Such a mapping can then be used to
compute the torques needed to move the arm along a given trajectory. 
Another example is predictive soil mapping (Fig.~\ref{subfig:soil_model}), where one is given a set of soil samples taken from some regions, and asked to predict the nature of soil in another region. A major benefit of using Gaussian processes to solve these problems is that they can provide confidence measures for the predictions. For instance, in the context of predictive soil mapping, one can use Gaussian processes to decide which regions should be given a higher priority for collecting soil samples, based on the uncertainty of the predictions.
The following sections provide a mathematical treatment of Gaussian processes, and their application to regression problems.

\begin{comment}
\section{How to go about solving them?}

In general the procedure is to denote the input as $x$, and the output (or target) as $y$. 
The input is usually represented as a vector $x$ as there are in general many input
variables-in the handwritten digit recognition example one may have a 256-
dimensional input obtained from a raster scan of a $16 \times 16$ image, and in the
robot arm example there are three input measurements for each joint in the
arm. The target $y$ may either be continuous (as in the regression case) or
discrete (as in the classification case). We have a dataset $D$ of $n$ observations,
$D = {(x_i, y_i)|i = 1, . . . , n}$.

Given this training data we wish to make predictions of the ouput $y$ for 
new inputs $x^*$ that we have not seen in the training set. Thus it is clear that the problem
at hand is inductive; we need to move from the finite training data $D$ to a
function $f$ that makes predictions for all possible input values. To do this we
must make assumptions about the characteristics of the underlying function,
as otherwise any function which is consistent with the training data would be
equally valid.

This is where the Gaussian process comes to our rescue. A Gaussian process is 
a generalization of the Gaussian probability distribution.
\end{comment}

%\newpage

\section{Gaussian Process (GP)} 
A Gaussian process can be thought of as a Gaussian distribution over functions (thinking of functions as infinitely long vectors containing the value of the function at every input). Formally let the input space be $\mathcal{X}$ and $f : \mathcal{X} \rightarrow \mathbb{R}$ be a function from the input space to the reals. We then say that $f$ is a Gaussian process if for any vector of inputs $\textbf{x} = [x_1, x_2, \dots, x_n]^T$ such that $x_i \in \mathcal{X}$ for all $i$, the vector of outputs $f(\textbf{x})$ = $[f(x_1), f(x_2), \dots, f(x_n)]^T$ is Gaussian distributed.

A Gaussian process is specified by a mean function $\mu : \mathcal{X} \rightarrow \mathbb{R}$, such that $\mu(x)$ is the mean of $f(x)$ and a covariance/kernel function $k: \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}$ such that $k(x_i,x_j)$ is the covariance between $f(x_i)$ and $f(x_j)$. We say $f \sim GP(\mu,k)$ if for any $x_1, x_2, \dots x_n \in \mathcal{X}$, $[f(x_1), f(x_2), \dots, f(x_n)]^T$ is Gaussian distributed with mean $[\mu(x_1), \mu(x_2), \dots, \mu(x_n)]^T$ and $n \times n$ covariance/kernel matrix $K_{\textbf{x}\textbf{x}}$:

\begin{displaymath}
K_{\textbf{x}\textbf{x}} = \left[ \begin{array}{cccc}
k(x_1,x_1) & k(x_1, x_2) & \dots & k(x_1,x_n) \\
k(x_2,x_1) & k(x_2, x_2) & \dots & k(x_2,x_n) \\
\vdots & \vdots & \vdots & \vdots \\
k(x_n,x_1) & k(x_n, x_2) & \dots & k(x_n,x_n) \\
\end{array}
\right]
\end{displaymath}

\begin{figure}
\center
\includegraphics[width=4 cm]{covariance.png}
\caption{Magnitude of values in covariance matrix, $K_{\textbf{x}\textbf{x}}$, decrease as the points are further apart.}
\label{fig:covariance}
\end{figure}

The kernel function must have the following properties:
\begin{itemize}
\item Be symmetric. That is, $k(x_i,x_j) = k(x_j,x_i)$\
\item Be positive definite. That is, the kernel matrix $K_{\textbf{x}\textbf{x}}$ induced by $k$ for any set of inputs is a positive definite matrix.
\end{itemize}

Examples of some kernel functions are given below:

\begin{itemize}
\item Squared Exponential Kernel (Gaussian/RBF): $k(x_i,x_j) = \exp(\frac{-(x_i-x_j)^2}{2\gamma^2})$ where $\gamma$ is the length scale of the kernel.
\item Laplace Kernel: $k(x_i,x_j) = \exp(\frac{-|x_i-x_j|}{\gamma})$.
\item Indicator Kernel: $k(x_i,x_j) = I(x_i=x_j)$, where $I$ is the indicator function.
\item Linear Kernel: $k(x_i,x_j) = x_i^Tx_j$.
\end{itemize}

More complicated kernels can be constructed by adding known kernel functions together, as the sum of 2 kernel functions is also a kernel function.

A Gaussian process is a random stochastic process where correlation is introduced between 
neighboring samples (think of a stochastic process as a sequence of random variables). The 
covariance matrix $K_{\textbf{x}\textbf{x}}$ has larger values, for points that are closer to each other, and smaller values 
for points further apart. This is illustrated in Fig.~\ref{fig:covariance}. The thicker 
the line, the larger the values. This is because the points are correlated by the difference in 
their means and their variances. If they are highly correlated, then their means are almost same 
and their covariance is high.

\subsection{Visualizing samples from a Gaussian process}

\begin{figure}
\centering
\subfigure[Gaussian kernel, $\gamma=0.5$]{
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0.5cm, clip=true]{gaussian_kernel.pdf}
\label{subfig:gaussian_kernel}
}
\subfigure[Laplace kernel, $\gamma=0.5$]{
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0.5cm, clip=true]{laplace_kernel.pdf}
\label{subfig:laplace_kernel}
}
\subfigure[Linear kernel]{
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0.5cm, clip=true]{linear_kernel.pdf}
\label{subfig:linear_kernel}
}
\subfigure[Indicator kernel]{
\includegraphics[width=0.4\textwidth, trim = 0cm 0cm 0cm 0.5cm, clip=true]{indicator_kernel.pdf}
\label{subfig:indicator_kernel}
}
\caption{Samples from a Gaussian process}
\label{GP_samples}
\end{figure}

To actually plot samples from a Gaussian process, one can adopt the following procedure:
\begin{enumerate}
\item Define the mean function and kernel for the GP. For instance, $\mu=0$, and $k(x_i,x_j) = \exp(\frac{-(x_i-x_j)^2}{2\gamma^2})$, with $\gamma=0.5$.
\item Sample inputs $x_i$; example: $x_i = \varepsilon\times i$, $i=0,1,\ldots,\frac{1}{\varepsilon}$
\item Compute the kernel matrix $\Sigma$. For the example with $\varepsilon=0.1$, we would have a $11\times 11$ matrix.
\item Sample from the multivariate Gaussian distribution $N({\bf 0},\Sigma)$
\item Plot the samples
\end{enumerate}
Fig.~\ref{GP_samples} shows examples of samples drawn from a Gaussian process, for different choice of kernels and $\varepsilon=0.01$.

\section{Inference}

Gaussian processes are useful as priors over functions for doing non-linear regression.
In Fig.~\ref{subfig:figInference_a}, we see a number of sample functions
drawn at random from a prior distribution over functions specified by a 
particular Gaussian process, which favours smooth functions. This prior is taken
to represent our prior belief over the kinds of functions we expect to observe,
before seeing any data. Note that $\frac{1}{N}\sum\limits_{i=1}^N f_i(x)\rightarrow \mu_f(x)=0$, as $N\rightarrow \infty$.

\begin{figure}
\subfigure[Samples from a zero-mean GP prior]{
\includegraphics[width=0.5\textwidth]{gp_prior.png}
\label{subfig:figInference_a}
}
\subfigure[Samples from the posterior after a few observations]{
\includegraphics[width=0.5\textwidth]{gp_posterior.png}
\label{subfig:figInference_b}
}
\caption{Gaussian process inference}
\label{figInference}
\end{figure}

Now, given a set of observed inputs and corresponding output values $(x_1, f(x_1))$, $(x_2, f(x_2))$, $\dots, (x_n, f(x_n))$, and a Gaussian 
process prior on $f$, $f \sim GP(\mu,k)$, we would like to compute 
the posterior over the value $f(x^*)$ at any query input $x^*$. 
Figure \ref{subfig:figInference_b} shows sample functions drawn from the posterior, given some observed $(x,y)$ pairs. We see that the sample functions from the posterior pass close to the observed values, but vary a lot in regions where there are no observations. 
This shows that uncertainty is reduced near the observed values.


\subsection{Computing the Posterior}

The posterior can be derived similarly to how the update equations for the Kalman filter were derived. First, we will find the joint distribution of $[f(x^*), f(x_1), f(x_2), \dots, f(x_n)]^T$, and then use the conditioning rules for a Gaussian to compute the conditional distribution of $f(x^*)|f(x_1), \dots, f(x_n)$.

Assume for now that the prior mean function $\mu = 0$. By definition of the Gaussian process, the joint distribution $[f(x^*), f(x_1), f(x_2), \dots, f(x_n)]^T$ is a Gaussian:

\begin{displaymath}
\begin{array}{rcl}
\left[ \begin{array}{c}
f(x^*)\\
f(x_1)\\
\vdots\\
f(x_n)\\
\end{array}
\right]
&
\sim
&
N \left(
\left[ \begin{array}{c}
0\\
0\\
\vdots\\
0\\
\end{array}
\right],
\left[ \begin{array}{cc}
k(x^*,x^*) & k(x^*,\textbf{x})^T\\
k(x^*,\textbf{x}) & K_{\textbf{x}\textbf{x}}\\
\end{array}
\right] \right)
\end{array}
\end{displaymath}

where $K_{\textbf{x}\textbf{x}}$ is the kernel matrix defined previously, and

\begin{align*}
k(x^*,\textbf{x}) = \left[ \begin{array}{c}
k(x^*,x_1) \\
k(x^*,x_2) \\
\vdots \\
k(x^*,x_n) \\
\end{array} \right]
\end{align*}

Using the conditioning rules we derived for a Gaussian,  the posterior for $f(x^*)$ is:

\begin{align*}
f(x^*)|f(\textbf{x}) & \sim N\left(k(x^*,\textbf{x})^T K_{\textbf{x}\textbf{x}}^{-1} f(\textbf{x}),  k(x^*, x^*) + k(x^*,\textbf{x})^T K_{\textbf{x}\textbf{x}}^{-1} k(x^*,\textbf{x})\right)\\
\end{align*}

The posterior mean $\mathbb{E}(f(x^*)|f(\textbf{x}))$ can be interpreted in two ways. We could group the last two terms $K_{\textbf{x}\textbf{x}}^{-1} f(\textbf{x})$ together and represent the posterior mean as linear combination of the kernel function values:

\begin{displaymath}
\mathbb{E}(f(x^*)|f(\textbf{x})) = \sum_{i=1}^n \alpha_i k(x^*,x_i)
\end{displaymath}

for $\alpha = K_{\textbf{x}\textbf{x}}^{-1} f(\textbf{x})$. This means we can compute the mean without explicitly inverting $K$, by solving $K \alpha = f(\textbf{x})$ instead. Similarly, by grouping the first two terms $k(x^*,\textbf{x})^T K_{\textbf{x}\textbf{x}}^{-1}$, the posterior mean can be represented as a linear combination of the observed function values:

\begin{displaymath}
\mathbb{E}(f(x^*)|f(\textbf{x})) = \sum_{i=1}^n \beta_i f(x_i)
\end{displaymath}

for $\beta = k(x^*,\textbf{x})^T K^{-1}$.

\subsection{Non-zero mean prior}

If the prior mean function is non-zero, we can still use the previous derivation by noting that if $f \sim GP(\mu,k)$, then the function $f' = f - \mu$ is a zero-mean Gaussian process $f' \sim GP(0,k)$. Hence, if we have observations from the values of $f$, we can subtract the prior mean function values to get observations of $f'$, do the inference on $f'$, and finally once we obtain the posterior on $f'(x^*)$, we can simply add back the prior mean $\mu(x^*)$ to the posterior mean, to obtain the posterior on $f$.

\subsection{Noise in observed values}

If instead of having noise-free observations of $f$, we observe $y(x) = f(x) + \varepsilon$, where $\varepsilon \sim N(0,\sigma^2)$ is some zero-mean Gaussian noise, then the joint distribution of $[f(x^*), y(x_1), \dots y(x_n)]^T$ is also Gaussian. Hence, we can apply a similar derivation to compute the posterior of $f(x^*)$. Specifically, if the prior mean function $\mu = 0$, we have that:

\begin{displaymath}
\begin{array}{rcl}
\left[ \begin{array}{c}
f(x^*)\\
y(x_1)\\
\dots\\
y(x_n)\\
\end{array}
\right]
&
\sim
&
N \left(
\left[ \begin{array}{c}
0\\
0\\
\dots\\
0\\
\end{array}
\right],
\left[ \begin{array}{cc}
k(x^*,x^*) + \sigma^2 & k(x^*,\textbf{x})^T\\
k(x^*,\textbf{x}) & K_{\textbf{x}\textbf{x}} + \sigma^2 I\\
\end{array}
\right] \right)
\end{array}
\end{displaymath}

The only difference with respect to the noise-free case is that the covariance matrix of the joint now has an extra $\sigma^2$ term on its diagonal. This is because the noise is independent for different observations, and also independent of $f$ (so no covariance between noise terms, and between $f$ and $\varepsilon$). So the posterior on $f(x^*)$ is:

\begin{align*}
f(x^*)|y(\textbf{x}) & \sim  N\left( k(x^*,\textbf{x})^T (K_{\textbf{x}\textbf{x}} + \sigma^2I)^{-1} y(\textbf{x}) , k(x^*, x^*) + \sigma^2 + k(x^*,\textbf{x})^T ( K_{\textbf{x}\textbf{x}} + \sigma^2I)^{-1} k(x^*,\textbf{x}) \right)\\
\end{align*}

\subsection{Choosing kernel length scale and noise variance parameters}

The kernel length scale ($\gamma$) and noise variance ($\sigma^2$) parameters are chosen such that they maximize the log likelihood of the observed data. Assuming a Gaussian kernel, we obtain the most likely parameters $\gamma$ and $\sigma$ by solving:

\begin{align*}
\max_{\gamma, \sigma} \log P(y(\textbf{x})|\gamma,\sigma) &= \max_{\gamma,\sigma} \left( - \frac{1}{2} y(\textbf{x})^T(K_{\textbf{x}\textbf{x}}+\sigma^2I)^{-1} y(\textbf{x}) - \frac{1}{2} \log( \textrm{det}(K_{\textbf{x}\textbf{x}} + \sigma^2 I))  - \frac{N}{2} \log(2\pi) \right)
\end{align*}

Here, the determinant term will be small when $K_{\textbf{x}\textbf{x}}$ is almost diagonal; thus this maximization favors smoother kernels (larger $\gamma$). Additionally $\sigma^2$ can be chosen to have a higher value to prevent overfitting, since larger values for $\sigma$ mean we trust observations lesser.

\subsection{Computational complexity}

One drawback of Gaussian processes is that it scales very badly with the number of observations $N$. Solving for the coefficients $\alpha$ that define the posterior mean function requires $O(N^3)$ computations. Note that Bayesian Linear Regression (BLR), which can be seen as a special case of GP with the linear kernel, has complexity of only $O(d^3)$ to find the mean weight vector, for a $d$ dimensional input space $\mathcal{X}$. Finally, to make a prediction at any point, Gaussian process requires $O(N\hat{d})$ (where $\hat{d}$ is the complexity of evaluating the kernel), while BLR only requires $O(d)$ computations.

\begin{thebibliography}{99}
\bibitem{JanPeters}
Botond Bocsi, Duy Nguyen-Tuong, Lehel Csató, Bernhard Schölkopf, Jan Peters, ``Learning inverse kinematics with structured prediction,'' in {\it IROS 2011: 698-703}
\bibitem{JuanPablo}
Juan Pablo Gonzalez, Simon Cook, Thomas Oberthur, Andrew Jarvis, J. Andrew (Drew) Bagnell, and M Bernardine Dias, ``Creating Low-Cost Soil Maps for Tropical Agriculture using Gaussian Processes,'' in {\it Workshop on AI in ICT for Development (ICTD) at the Twentieth International Joint Conference on Artificial Intelligence (IJCAI 2007), January, 2007.}
\end{thebibliography}

\end{document}
