<commit_before>\documentclass[twoside]{article}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{amssymb,amsmath,amsthm}
\usepackage{graphicx}
\usepackage{preamble}
\usepackage{natbib}
\usepackage{hyperref}
\usepackage{color}
\usepackage{wasysym}
\usepackage{subfigure}
\usepackage{tabularx}
\usepackage{booktabs}
\usepackage{bm}
\newcommand{\theHalgorithm}{\arabic{algorithm}}
\definecolor{mydarkblue}{rgb}{0,0.08,0.45}
\hypersetup{ %
    pdftitle={},
    pdfauthor={},
    pdfsubject={},
    pdfkeywords={},
    pdfborder=0 0 0,
    pdfpagemode=UseNone,
    colorlinks=true,
    linkcolor=mydarkblue,
    citecolor=mydarkblue,
    filecolor=mydarkblue,
    urlcolor=mydarkblue,
    pdfview=FitH}

\newcolumntype{x}[1]{>{\centering\arraybackslash\hspace{0pt}}m{#1}}
\newcommand{\tabbox}[1]{#1}


\setlength{\marginparwidth}{0.6in}
\input{include/commenting.tex}

\usepackage{format/icml2013}
%\usepackage[left=1.00in,right=1.00in,bottom=0.25in,top=0.25in]{geometry} %In case we want larger margins for commenting purposes

%% For submission, make all render blank.
%\renewcommand{\LATER}[1]{}
%\renewcommand{\fLATER}[1]{}
%\renewcommand{\TBD}[1]{}
%\renewcommand{\fTBD}[1]{}
%\renewcommand{\PROBLEM}[1]{}
%\renewcommand{\fPROBLEM}[1]{}
%\renewcommand{\NA}[1]{#1}  %% Note, NA's pass through!
    
\begin{document}



\twocolumn[
\icmltitle{Nonparametric Structure Discovery in Supervised Regression through Compositional Kernel Search}
%\icmlauthor{Anonymous Authors}
%\icmladdress{ Unknown Institution}
%%%%FIXME - Make sure to update this!
\icmlkeywords{nonparametrics, gaussian process, machine learning, ICML, structure learning, extrapolation, regression, kernel learning, equation learning, supervised learning}
\vskip 0.3in
]

\begin{abstract}
%Gaussian process (GP) models are used widely and successfully.
%However, their effictiveness depends critically on choosing an appropriate family of kernels.
%This aspect of GP modeling has been sorely underdeveloped.
%In this paper, we introduce a procedure for automatically and efficiently searching through a large space of GP models.
The effectiveness of nonparametric regression models depends heavily on the choice of kernel.
We introduce a marginal-likelihood-based search over composite kernel structures which automatically constructs a structured Gaussian process model appropriate for the dataset.
We further demonstrate that such kernels often allow the posterior to be automatically decomposed into a sum of interpretable components, and in some cases allows long-range extrapolation.
We demonstrate this technique on several real datasets, and achieve state-of-the-art predictive performance.
\end{abstract}

\section{Introduction}

Supervised learning problems, such as classification and regression, learn a function ${\function : \InputSpace \to \OutputSpace}$ from some input (predictor) variables, $\InputVar$, to some output (response) variables, $\outputVar$.
Kernel-based nonparametric models, such as support vector machines and Gaussian processes, have been one of the dominant paradigms for supervised machine learning over the last 20 years.
These methods depend on defining a kernel or covariance function, $\kernel : \InputSpace \times \InputSpace \to \Reals$ which specifies how similar or correlated outputs $\outputVar$ and $\outputVar'$ are at two inputs $\inputVar$ and $\inputVar'$, respectively.

A key challenge for kernel based methods is learning an appropriate kernel from data, and a great many papers have been written on this important topic\TBD{ cite kernel learning for SVMs and GP literature}.
In this paper, we pose the problem of kernel learning as a problem of structure discovery from data.
Specifically, we focus on a Bayesian setting where the kernel specifies a covariance function for Gaussian process (GP) regression.

Traditionally, Gaussian process kernels have simple forms, such as squared-exponential or Matern, which depend on few free hyperparameters which can be optimised (or inferred) from data.
However, one can also straightforwardly form composite kernels by combining simple base kernels through operations that such as summation, multiplication, and change of coordinates, which all preserve positive definiteness.
We use this insight to define a simple grammar over composite kernels, and we develop an automated search algorithm over the (exponentially large) space of kernels that can be derived from this grammar.
The search criterion is the marginal likelihood of the data for a given kernel, which is the driving term for Bayesian model comparison and selection in analogous structure discovery tasks such as graphical model learning \cite{heckerman1995learning}.

%\paragraph{Benefits of kernel structure}  In many regression problems, or time series analysis, or when modeling dynamical systems, we often think we are observing the superposition of multiple distinct causal processes.
%For example, an underlying linear function plus Gaussian noise in classical linear regression, or an underlying linear dynamics plus Gaussian observation noise in a Kalman filter.
%Our approach can be seen as a general way to learn functions like this, where the product of kernels can be used to capture richer causal proceses (e.g., nonlinear rather simply linear, amplifying perioidic (as in the second layer of the airline example) rather than simply periodic) and the sum of kernels corresponds to the superposition of these causal processes.

Our experimental findings demonstrate that this structure discovery algorithm can automatically recover known kernel structure from data, and that on real data the learned structure of the kernel often corresponds to natural, interesting and interpretable decompositions of the unknown function.
Our hope is that the algorithm developed in this paper will help replace the current and often opaque art of kernel engineering with a more transparent science of automated kernel discovery\fTBD{maybe $\to$ conclusion}.

%Kernel-based nonparametric regression models have been widely and succesfully used over the last 20 years. [citations needed]
%As part of the development of this model class, a rich set of kernels have been developed for both continuous and structured data.
%Because different kernels reflect different properties of the function being modeled, the choice of kernel can have a large impact on the sorts of structure that can be expressed or discovered by the model.

%When approaching a new dataset, it is often not clear which kernel family is most appropriate, especially for non-experts.
%Often, a standard kernel family (typically squared-exponential) is chosen out of convenience, and kernel parameters are set by cross-validation or maximum marginal likelihood. [citations?]

%Sometimes the argument is made that since most kernels are universal (i.e. consistent in the limit), that the choice of kernel is not important[citation?].  However, our experiments demonstrate that for datasets of moderate size, predictive performance, extrapolation and interpretability all depend heavily on the choice of kernel. \TBD{RBG: Do people actually say the choice of kernel isn't important?  This sounds like a straw man...}

%\paragraph{Main contribution}
%In this paper, we introduce an automated search over composite kernels.
%We define the search space through a simple grammar consisting of a set of one-dimensional base kernels (squared-exp, periodic, linear, etc.) and operations which combine those kernels (add and multiply).
%Our objective function is the marginal likelihood of a Gaussian process model with the given kernel, conditioned on the data.

%We examine the composite kernels discovered by this search on several datasets, and demonstrate that they can recover known structure, and result in state-of-the-art predictive performance.  We further demonstrate that in some cases, the structured posterior can be automatically decomposed into a sum of interpretable components, e.g. in Figure \ref{fig:airline_decomp}. \TBD{RBG: Need to clarify that the contribution isn't that we can decompose it (that was in GPML), but that the decomposition is often interesting.}

%\section{Motivation}
%\TBD{RBG: Do we need a separate motivation section?  The introduction should motivate the work.}

%There is a large and mature literature on automatic structure discovery in unsupervised settings [citations including Josh].
%However, relatively little attention has been paid to automatic structure discovery in supervised learning.

%Similar searches over large model classes have been succesfully used in machine vision [cite Cox + Pinto].
%  This work respresents a step towards more data-based modeling, as opposed to proposing the model beforehand.
%In general, learning the model class from data seems superior proposing the model beforehand.
%In high dimensional problems, it is also hard for a practitioner to propose an appropriate model even after examining a dataset closely.
% Choosing a kernel family is also a stumbling block for non-experts who wish to use Gaussian Process models.

%Machine learning can be more data-driven, analogous to the high-thoughput approaches being used in biology. 

%\subsection{Convenience Priors}
%One of the main advantage of kernel methods, and Gaussian process regression in particular, is that the kernel can be used to specify interesting prior structure to the model.
%However, in practice, GPs are typically used with only the squared-exp kernel.
%[Apparently, this was one of Sholkopf's big dissapointment with GPs... is there a quote from him on this? -David]

%A frequent criticism of the Bayesian framework is that priors are often chosen purely for convenience. 
%Gaussian process priors may be a good examples of a model chosen for its tractability.
%However, specifying a compound kernel typically does not change the order of computational complexity of inference in Gaussian process models.
%The broad use of simple kernels may be simply due to the difficulty of determining which kinds of structure exist in the data.
%As well, even if one understands the structure present in a dataset, it is often non-obvious how to express that through a covariance function. 
%\TBD{RBG: Not sure we need this in the motivation. Is this basically just to answer the objection of why we don't consider the tradeoff between expressiveness and tractability?  Also, the prior can affect computational tractability if we're using approximate algorithms, so this makes it sound like we're only thinking about small datasets.}

%Automating the search over kernel structures should be especially useful in allowing pracitioners from a wide variety of backgrounds to use appropriate models for their data.
%Indeed, the authors were sometimes surprised by the kernels chosen by the automated search, and the elegant ways in which complex structure was expressed through simple combination of kernels.

\section{Expressing structure through kernels}
\label{sec:Structure}

Kernel functions $\kernel : \InputSpace \times \InputSpace \to \Reals$ can be used to define a measure of similarity between two points $\inputVar, \inputVar'$ in some space $\InputSpace$.
In the case of Gaussian process regression, the kernel defines a covariance function ${\KernelMatrix_{\inputVar \inputVar'} := \textrm{Cov}(\outputVar, \outputVar') = \kernel(\inputVar,\inputVar')}$.
%The kernel, $\kernel$, must define a valid covariance function\fTBD{expand me}; when this is the case $\kernel$ is said to be positive semi-definite (PSD).

Examples of commonly used kernels include squared exponential (SE), periodic (Per) and linear kernels (Lin) defined below\fTBD{If we add RQ then we have everything in one place}
\begin{eqnarray}
\kernel_\textrm{SE}(\inputVar, \inputVar') = & \sigma^2\exp\left(-\frac{(\inputVar - \inputVar')^2}{2\ell^2}\right) \\
\kernel_\textrm{PE}(\inputVar, \inputVar') = & \sigma^2\exp\left(-\frac{2\sin^2(\pi|\inputVar - \inputVar'|/p)}{\ell^2}\right) \\
\kernel_\textrm{LN}(\inputVar, \inputVar') = & \sigma_1^2 + \sigma_2^2(\inputVar - \ell)(\inputVar' - \ell).
\end{eqnarray}

When used within the context of GP regression, these kernels define priors on functions which place mass only on smooth, periodic or linear functions respectively.
\TBD{Reference a diagram - even if a composite kernel we can point out what the base kernel is doing.}

Valid kernels are closed under addition and multiplication with other valid kernels;
% \ie any algebraic composition of PSD kernels will define a PSD kernel.\fTBD{Cite theorem}
this allows one to create richly structured kernels from well understood base components.

\fTBD{Can probably kill simple kernel table. Focus on one table.}
%\input{tables/simple_kernels_table.tex}
\input{tables/example_structures_table.tex}

\paragraph{Summation}

Summation of kernels corresponds to the summation of independent functions in the following sense.
Suppose ${\function_1 \dist \GP(0, \kernel_1)}$, and ${\function_2 \dist \GP(0, \kernel_2)}$ and are independent.
Then ${\function := \function_1 + \function_2 \dist \GP(0, \kernel_1 + \kernel_2)}$\footnotemark.
\footnotetext{Additionally, the posterior of the component functions conditioned on observations of the sum are analytically tractable. We make use of this in section\TBD{ add reference} and give derivations in the supplementary material.}

In one dimension, sums of kernels can express structure such as variational at multiple scales, or of summation of different function types.  Some examples can be found in Figure \ref{fig:kernels}.

In higher dimensions, summing kernels can correspond to decompositions of the form
\begin{equation}
\function((\inputVar_1, \inputVar_2)) = \function_1(\inputVar_1) + \function_2(\inputVar_2)
\end{equation}
\ie `independent' variation in different dimensions.
%
\input{tables/multi_d_additive_fig}
%
Figure \ref{fig:multi_d_additivity} demonstrates the decomposition across dimensions.

\paragraph{Interpretation}
Loosely speaking, a sum of kernels can be understood as an 'OR' operation: two points are considered similar if either kernel has a high value.
More precisely, we can interpret a sum of kernels as a concatenation of their features in the corresponding RKHS.

\paragraph{Multiplication}

Multiplication of kernels does not correspond to a simple composition of functions; here, we give some examples of the properties that arise from multiplication of kernels.

%but allows for the construction of different structures and is best understood on a case by case basis\fTBD{Better language?}.
%
%For example, multiplying base kernels of the same form that act upon different dimensions of $\InputSpace$ can create the multidimensional generalisations of these kernels; this is true of \eg SE and LN\TBD{ reference a picture}.
%
%Multiplying kernels of the same form in one dimension is often redundant \eg for SE kernels this just produces SE kernels with different parameters.
%\TBD{What is PE $\times$ PE?}
%However, multiplying kernels of different forms can give rise to interesting structures.

Multiplying a kernel $\kernel$ by the SE kernel allows the structures implied by $\kernel$ to vary smoothly according the lengthscale of the SE kernel.  This is shown in row 2 of Figure \ref{kernels}.
Multiplying a kernel by LN results in the amplitude of structures implied by $\kernel$ growing or shrinking linearly. This is shown in rows 4 and 5 of Figure \ref{kernels}.
%\fTBD{Could perhaps perversely claim that multiplying by periodic results in wiggles being added to previous structure?}
%
%\fTBD{Can we do more than just give examples}

%In general, multiplying kernels can be thought of as an `AND' operation when viewing kernels as a measure of similarity.
This arises since two points, $\inputVar,\inputVar'$, can only be considered similar by a product kernel function (\ie $\kernel_\textrm{prod}(\inputVar,\inputVar')$ is relatively large) if the points are considered similar by each component kernel function.

\paragraph{Interpretation}
Loosely speaking, multiplying kernels can be understood as an 'AND' operation: two points are considered similar only if both kernels have a high value.
More precisely, we can interpret a product of kernels as a concatenation of the cartesian product of both kernels' features in the corresponding RKHS.

To develop an intuition for the effect of multiplying two stationary kernels, we can note that the posterior mean function must be a weight sum of kernel functions.  Thus the third column of rows 1-2 of Figure \ref{kernels} show the functions that can be combined linearly to produce the posterior means in the 5th column.

%\section{Gaussian Processes Priors}

%Gaussian processes are a flexible and tractable prior over functions, useful for solving regression and classification tasks\cite{rasmussen38gaussian}.
%The kind of structure which can be captured by a GP model is mainly determined by its \emph{kernel}: the covariance function.
%One of the main difficulties in specifying a Gaussian process model is in choosing a kernel which can represent the structure present in the data.
%For small to medium-sized datasets, the kernel has a large impact on modeling efficacy.
%\TBD{Note: The above paragraph is plagarized from my additive GP paper. -David} 

%
%The technique of constructing composite kernels using sums and products of existing kernels is not new \cite{rasmussen38gaussian} [more cites, <NAME>'s astronomy work?].  
%However, the main contribution of this paper is to automate the search over kernel structures.

\section{Searching over structures}

\paragraph{A grammar for kernels}

The discussion above demonstrates that rich structure can be captured with summation and multiplication of simple base kernels.
We therefore consider searching over all kernel structures that can be expressed as sums and products of base kernels.

\begin{center}
\begin{tabular}{rccc}
\textrm{Replacement} & $\kernel_i$ & $\to$ & $\kernel'_i$\\% & $\forall\, \kernel' $\\
\textrm{Addition} & $\kernel_i$ & $\to$ & $\kernel_i + \kernel'_j$\\% & $\forall\, j,\kernel' $\\
\textrm{Multiplication} & $\kernel_i$ &  $\to$ & $\kernel_i \times \kernel'_j$\\% & $\forall\, j,\kernel'$\\
\end{tabular}
\end{center}

Details of the search algorithm are given in the supplementary material \TBD{For now}.

\paragraph{A greedy search algorithm}
In order to find the model with highest marginal likelihood, the base kernels and production rules generate a graph structure over kernel expressions and we follow a greedy search to traverse this graph.
This requires a method of scoring a kernel expression, which we achieve by using the Bayesian Information Criterion (BIC)\TBD{ cite me} as an estimate of the marginal likelihood of some data given the kernel expression\footnotemark.
\footnotetext{We also experimented with using a Laplace approximation to the marginal likelihood\TBD{ cite} but this was found to be less numerically stable and could not be computed if an optimiser failed to find a true optimum.}
Specifically, we optimise the marginal likelihood of the data given specific kernel parameters\TBD{ using LBFGS - cite?}, randomly restarting previously unoptimised kernel parameters.
This is then adjusted by the number of parameters optimised\footnotemark to compute the BIC\fTBD{formula?}.
\footnotetext{Adjustments were made for obvious parameter redundancy \eg multiplicative terms in products of kernels.}

The search then proceeds by scoring all base kernels applied to all input dimensions.
The highest scoring kernel is then expanded by applying all production rules, and these kernel expressions are scored.
The production rules are then applied to the highest scoring kernel expression and the search continues iteratively until some fixed depth.

\paragraph{Example expressions}

\TBD{Draw a picture to show what our model includes \ie polynomial regression, nonparametric `Fourier' decomposition, GAM, additive kernels.}

\section{Related Work}

\paragraph{Composite kernels in GP models} The technique of constructing composite kernels using sums and products of existing kernels was demonstrated in detail in Chapter 5 of \cite{rasmussen38gaussian}, where the resulting posterior mean was also decomposed into a sum of component-wise means, although the posterior variance was not.  Our work can be viewed as an extension and automation of this approach.

\paragraph{Generalized Additive Models} \cite{hastie1990generalized} are models for which the posterior mean can be expressed $\expect[f(\vx)] = \sum_{d=1}^D f_d(x_d)$.  These models have a limited compositional form, but one which is very interpretable.  We compare against this model in our multi-dimensional experiments.

\cite{plate1999accuracy} constructs a \gp{} with a composite kernel, summing an SE kernel along each dimension, with an SE-ARD kernel on all dimensions.  This model is motivated by the desire to specify the trade off between the interpretability of additive models and the flexibility of multiplicative models.

\paragraph{Compositional model search for unsupervised learning} \cite{grosse2012exploiting} performed a greedy search over a compositional model class for unsupervised learning.  Their models were composed of sums and element-wise products of random variables defined on matrices.  This model class contained a large number of existing unsupervised models as special cases.  Earlier work in unsupervised Bayseian structure discovery was done in \cite{kemp2008discovery}.

\paragraph{ANOVA Kernels}

Support vector regression with ANOVA decomposition kernels \cite{stitson1999support}

%\subsubsection{Smoothing spline ANOVA models}
A closely related procedure from the statistics literature is smoothing-splines ANOVA (SS-ANOVA)\cite{wahba1990spline, gu2002smoothing}.
This model is a weighted sum of splines along each dimension, plus a sum of splines over all pairs of dimensions, all triplets, etc, with each individual interaction term having a separate weighting parameter.
Learning in SS-ANOVA is usually done via penalized-maximum likelihood with a fixed sparsity hyperparameter.
Because the number of terms to consider grows exponentially in the order, in practice, only terms of first and second order are usually considered.  This model is similar to ours, if we were to restict our attention to spline kernels

\paragraph{Additive Gaussian Processes} \cite{duvenaud2011additive11} is a Gaussian process model whose kernel implicitly sums over all $2^D$ possible products of one-dimensional base kernels.  
However, this model's parameterization is not as flexible as the kernels expressible in our grammar.

\paragraph{Hierarchical Kernel Learning}

%In "High-Dimensional Non-Linear Variable Selection through Hierarchical Kernel Learning", Bach
\cite{DBLP:journals/corr/abs-0909-0844} uses a regularized optimization framework to learn a weighted sum over an exponential number of kernels which can be computed in polynomial time.
The subsets of kernels considered by this method are restricted to be a \textit{hull} of kernels. In the setting we are considering in this paper, a hull can be defined as a subset of all terms such that if term $\prod_{j \in J} k_j(\bf x, x')$ is included in the subset, then so are all terms $\prod_{j \in J / i} k_j(\bf x, x')$, for all $i \in J$.
%For details, see \cite{DBLP:journals/corr/abs-0909-0844}.
%Given each dimension's kernel, and a pre-defined weighting over all terms, HKL performs model selection by searching over hulls of interaction terms.
% \subsubsection{All-subsets kernel with uniform weightings}
%In \cite{DBLP:journals/corr/abs-0909-0844}, Bach also fixes the relative weighting between orders of interaction with a single term $\alpha$, computing the sum over all orders by:
%\begin{equation}
%\label{eqn:uniform}
%k_{a}({\bf x, x'}) = v_D^2 \prod_{d=1}^D \left(1 + \alpha k_{d}(x_{d}, x_{d}') \right)
%\end{equation}
%which has computational complexity $O(D)$.
%However, this formulation forces the weight of all $n$th order terms to be weighted by $\alpha^n$.
%
The main difficulty with this approach is that hyperparameters are hard to set other than by cross-validation. 
%In contrast, our method optimizes the hyperparameters of each dimension's base kernel, as well as the relative weighting of each order of interaction. 



%A related functional ANOVA GP model\cite{kaufman2010bayesian} decomposes the \emph{mean} function into a weighted sum of GPs.
%However, the effect of a particular degree of interaction cannot be quantified by that approach.
%Also, computationally, the Gibbs sampling approach used in \cite{kaufman2010bayesian} is disadvantageous.

\paragraph{Genetic Searches over Kernel Structures}
A similar method was developed by \cite{diosan2007evolving} using SVMs in place of Gaussian processes, cross-validation error instead of marginal likelihood, and a genetic algorithm instead of a greedy search.  
In \cite{bing2010gp}, a genetic programming construction and optimization method was used to select a kernel for relevance vector machines.

\paragraph{Equation Learning}
A related body of work with similar goals are the various approaches to equation learning.  
For example, \cite{dzeroski1999equation} and \cite{washio1999discovering} attempt to learn parametric forms of equations to describe time series, or relations between quantities.  The advantage of a \gp{} approach is that 

\paragraph{Semiparametric Regression}
Semiparametric regression \TBD{add cites} attempts to combine interpretability with flexibility by building  a composite model out of an interpretable, parametric part (such as linear regression) and a ``catch-all'' nonparametric part (such as a \gp{} with an SE kernel).
Structure search has the same goal, but improves upon semiparametric regression in two ways.
First, construcing a compound model through structure in the kernel allows us to integrate over the parameters of the implicit components of the model.
Second, by automatically searching over components, we remove the limitation that the practictioner must be able to specify in advance the parametric form of any part of the model.


\paragraph{Multiple Kernel Learning}
There is a large body of work attempting to construct a rich kernel through a weighted sum of base kernels, including in \gp{} models.  
\cite{christoudias2009bayesian}
\TBD{cite more} 
However, these methods typically construct all kernels in advance, and simply reweight them.  This sort of approach would be inapplicable in our case, since our search space is exponential in our search depth, and we also must estimate the hyperparameters of each kernel.%One could imagine expanding our search for some fixed depth, instantiating a sum of all
%In the \gp{} framework,  showed how mixtures of kernels can be learnt by gradient descent in the Gaussian process framework.
%They call this \emph{Bayesian localized multiple kernel learning}.

\paragraph{Nonparametric kernel functions}

Since the complexity of our kernel is bounded only by the search depth, we could consider the structure search to be a non-parametric kernel model - in a sense, a doubly nonparametric model.  However, there has been some similar work along these lines in \cite{salakhutdinov2008using}, where a deep neural network was use to learn the kernel of a \gp{}.

%Hyperkernels \cite{ong2002hyperkernels}

\section{Decomposing a function}

Sums of kernels were shown in section~\ref{sec:Structure} to correspond to sums of independent functions.
By splitting kernel expressions found by our search algorithm into additive components, we can also decompose the \gp{} posterior over functions into a joint distribution over functions that are summed togther.
This potentially allows us to probabilistically separate causally the independent processes that give rise to the data.

Moreover, the posterior distribution over additive components is analytically tractable, allowing us to plot the individual components.
%\fTBD{If kept here, make notation consistent through paper}
The conditional distribution of a Gaussian $\function_1$ conditioned on its sum with another Gaussian\fTBD{Make sure notation strictly correct and defined} $\vf = \vf_1 + \vf_2$ where $\vf_1 \sim \GP( \vmu_1, k_1)$ and $\vf_2 \sim \GP( \vmu_2, k_2)$ is given by
\begin{align}
\vf_1 | \vf \sim \mathcal{N} \big( & \vmu_1 + \vk_1\tra (\vK_1 + \vK_2)\inv \left( \vf - \vmu_1 - \vmu_2 \right), \nonumber \\
& \vk_1 - \vk_1\tra (\vK_1 + \vK_2)\inv \vk_1 \big).
\end{align}
%and the covariance between the two components, conditioned on their sum is given by:
%\begin{align}
%\cov(\vf_1^\star, \vf_2^\star) | \vf = \vk_1^{\star\tra} (\vK_1 + \vK_2)\inv \vk_2^\star
%\end{align}
%Derivations can be found in the supplementary material.


We demonstrate the decompositions produced by our method on two time series.
We ran the search to a depth of 10 and used sq-exp, periodic, linear and rational quadratic kernels.

% \input{"tables/kernel_descriptions.tex"}

%The sorts of structure implied by complex kernel might not be immediately interpretable.
%However, there is one property of GPs that allows us to break down the signal into interpretable parts.
%Specifically, a function $f \sim \GP(0, k_1(x, x') + k_2(x,x')$, whose kernel is a sum of kernels, has the same distribution as a sum of two functions $\vf = \vf_1 + \vf_2$ where $\vf_1 \sim \GP( \vmu_1, k_1)$ and $\vf_2 \sim \GP( \vmu_2, k_2)$.
%In addition, the posterior distributions over individual components after conditioning on data are also analytically tractable.

\paragraph{Mauna Loa Atmospheric Carbon Dioxide}

We revisited a dataset explored in \cite{rasmussen38gaussian}, pages 120-126, where a kernel was hand-tailored to fit a GP model to the dataset.

Figures \ref{fig:mauna_grow} show the posterior of the increasingly appropriate model as the search depth increases on this dataset.
Note in particular that, while the data can be smoothly interpolated by a single base kernel model, the extrapolations improve dramatically as the search depth increases.

Figure \ref{fig:mauna_decomp} shows the complete posterior of the maximum depth model together with the posteriors of the additive components.
The plot shows a very plausible extrapolation and highly interpretable components i.e.~a long-term trend, annual periodicity, medium-term deviations and short-term serially correlated `noise'.\fTBD{Remember to update}

\begin{figure*}[h!]
\centering
\newcommand{\wmg}{5.5cm}  % width maunu growth
\newcommand{\hmg}{4cm}  % height maunu growth
\begin{tabular}{c}
 \includegraphics[width=\wmg,height=\hmg]{../figures/decomposition/03-mauna2003-s_max_level_0/03-mauna2003-s_all_small} 
 \includegraphics[width=\wmg,height=\hmg]{../figures/decomposition/03-mauna2003-s_max_level_1/03-mauna2003-s_all_small}
 \includegraphics[width=\wmg,height=\hmg]{../figures/decomposition/03-mauna2003-s_max_level_2/03-mauna2003-s_all_small}
\end{tabular}
\caption{Posterior mean and variance for different depths of kernel search.  In the first row, the function is only modeled as a locally smooth function, and the extrapolation is extemely poor.  In the second, a periodic component is added, and the extrapolation is better.  At depth 3, the kernel can capture most of the relevant structure, and is able to extrapolate reasonably. \TBD{RBG: (1) I think we somehow need to visualize the lengthscales, to make it obvious that the SE kernels really mean different things. (2) Why isn't SE + PE the correct answer?}}
\label{fig:mauna_grow}
\end{figure*}

\begin{figure}[h!]
\newcommand{\wmgd}{9.5cm}  % width mauna growth decomp
\newcommand{\hmgd}{3.1cm}  % height mauna growth decomp
\newcommand{\mdrd}{../figures/decomposition/11-Feb-03-mauna2003-s}  % mauna decomp results dir
\begin{tabular}{c}
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_all} \\ = \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_1} \\ + \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_2_zoom} \\ + \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_3} \\ + \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_resid}
\end{tabular}
\caption{First row: The posterior on the Mauna Loa dataset after a search of depth 8.  Subsequent rows: automatic decomposition of Mauna Loa data.  The signal has been automatically decomposed into long-term, yearly periodic, medium-term anomaly and short-term noise components, respectively.}
\end{figure}
\label{fig:mauna_decomp}

\paragraph{Airline passenger data}

\TBD{Give an intro to this dataset.}
Figures \ref{fig:airline_grow} and \ref{fig:airline_decomp} show the corresponding plots for the airline data set (\TBD{describe me!}).
\TBD{Make similar observations of good performance and perhaps draw attention to the selection of a linear kernel}.

\begin{figure*}[h!]
\centering
\newcommand{\wag}{4.8cm}  % width airline growth
\newcommand{\hag}{4cm}  % height airline growth
\begin{tabular}{c}
\includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_0/01-airline-s_all_small}
\hspace{-1cm} \includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_1/01-airline-s_all_small}
\hspace{-1cm} \includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_2/01-airline-s_all_small} 
\hspace{-1cm} \includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_3/01-airline-s_all_small}
\end{tabular}
\caption{Posterior mean and variance for different levels of kernel search on the airline dataset.}
\label{fig:airline_grow}
\end{figure*}


\begin{figure}[h!]
\centering
\newcommand{\wagd}{9.5cm}  % width airline decomp
\newcommand{\hagd}{3.5cm}  % height airline decomp
%\newcommand{\ard}{../figures/decomposition/11-Feb-01-airline-s}  % airline results dir
\newcommand{\ard}{../figures/decomposition/31-Jan-v201-airline-months}  % airline results dir
\begin{tabular}{c}
\hspace{-1cm}\includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_all} \\
 = \\ 
%\hspace{-1cm} \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-s_4} \\
% + \\
\hspace{-1cm} \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_1} \\
 + \\
\hspace{-1cm}  \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_2} \\
 + \\
\hspace{-1cm}   \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_3} \\
 + \\
\hspace{-1cm} \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_resid}
\end{tabular}
\caption{First row:  The airline dataset and posterior after a search of depth 8.  Rows 2-4: Additive decomposition of posterior into long-term smooth trend, yearly variation, and short-term noise.  Note that the variance of the noise grows over time, making this a heteroskedastic model.}
\label{fig:airline_decomp}
\end{figure}

%To evaluate the effectiveness of our method we performed two types of experiments.
%First, we examined the ability of our algorithm to discover useful structure in one-dimensional datasets.  Second, we examined the predictive performance of our model on multi-dimensional datasets.
% testing the ability of the algorithm to infer both an accurate regression function estimate and to produce an interpretable decomposition of a regression function.
%The accuracy experiments show that our method consistently matches or beats previous state of the art regression methods (\TBD{not yet statistically significantly within each experiment - more experiments currently running\ldots}).
%The decomposition experiments produce highly interpretable decompositions of time series data and highly plausible extrapolations; a particularly rare property of naively applied (linear) smoothers (\TBD{justify comment - e.g.~local linear is just linear, GP Sq-Exp nose dives towards the mean}).

\paragraph{Solar Irradiance Data}
\cite{lean1995reconstruction}
\TBD{Discuss failure mode}

\subsection{Multidimensional decomposition}

\begin{figure*}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.3\textwidth]{../figures/additive_multi_d/interpretable_1st_order1.pdf} &
\includegraphics[width=0.3\textwidth]{../figures/additive_multi_d/interpretable_1st_order2.pdf}& 
\includegraphics[width=0.3\textwidth]{../figures/additive_multi_d/interpretable_2nd_order1.pdf}\\
\end{tabular}
\caption{One-dimensional decompositions of the concrete dataset posterior.  Left, Centre:  Green points indicate the original data, blue points are data after the mean contribution from the other dimensions' first-order terms has been subtracted.  The black line is the posterior mean of a GP with only one term in its kernel.  Right:  The posterior mean of a GP with only one second-order term in its kernel.}
\label{fig:interpretable functions}
\end{figure*}
%
Figure \ref{fig:interpretable functions} demonstrates that an additive posterior can be decomposed into a sum of functions across dimensions.


\section{Empirical evaluation}

In addition to producing highly interpretable decompositions of regression functions, our method also produces state of the art predictions in both extrapolation and interpolation tasks.

\subsection{Extrapolation}

\begin{figure}[h!]
\centering
\begin{tabular}{c}
\hspace{-0.5cm} \includegraphics[width=\columnwidth]{../figures/extrapolation_curves/01-airline-s-ex-curve.pdf}
\end{tabular}
\caption{Extrapolation performance.  Here, we show the average test set MSE as a function of the fraction of the dataset used for training.}
\end{figure}

Figure \ref{extrapolation} shows the learning curves of our method compared against simple linear regression and \gp{} regression with a SE kernel.  The amount of data in the 1D airline timeseries was increased by increments of 10\%. Both linear and SE regression plateau, while unbounded complexity of our method allows the model to continue to capture structure. 

\subsection{Interpolation}

To demonstrate the interpolation accuracy of our method we extend the comparison of \cite{duvenaud2011additive11} to include our method.
5 comparsion algorithms were evaluated on 5 datasets in terms of mean-squared (MSE) error and predictive likelihood using 10 fold cross validation.
Results are presented in tables \ref{tbl:Regression Mean Squared Error} and \ref{tbl:Regression Negative Log Likelihood}.
%
\input{tables/regression_results_mse.tex}
\input{tables/regression_results_nll.tex}
%
Our method outperforms the other methods in all tests.

\TBD{Other versions of tests / learning curves}

\subsubsection{Data sets}


\paragraph{Bach Synthetic Dataset}
\fTBD{Too much text? Move to supplementary}
In addition to standard UCI repository datasets, we generated a synthetic dataset following the same recipe as \cite{DBLP:journals/corr/abs-0909-0844}.
% From a covariance matrix drawn from a Wishart distribution with 1024 degrees of freedom, we select 8 variables.
%We then construct the non-linear function $f(X) = \sum_{i=1}^4 \sum_{j=1+1}^4 X_i X_j + \epsilon$, which sums all 2-way products of the first 4 variables, and adds Gaussian noise $\epsilon$.
This dataset is one which can be predicted well by a kernel which is a sum of two-way interactions over the first 4 variables, ignoring the extra 4 noisy copies.
%
This dataset was designed by \cite{DBLP:journals/corr/abs-0909-0844} to demonstrate the advantages of HKL over a \gp{} with a SE-ARD kernel. 

%If the dataset is large enough, HKL can construct a hull around only those subsets of cross terms that are optimal for predicting the output.
%GP-ARD, in contrast, can only learn to ignore the noisy copy variables, but cannot learn to ignore the higher-term interactions between the predictive variables.
%However, a GP with an additive kernel can learn both to ignore irrelevant variables, and to ignore certain orders of interaction.
%In this example, the additive GP is able to recover the relevant structure.

[TODO: Describe concrete, pumadyn, servo and housing]

\subsection{Details of methods}

Linear regression was included as a baseline, with noise variance set to the empirical variance on the training residuals.
%
Generalized Additive Models (GP GAM) were implemented through a \gp{} whose kernel is a sum of sq-exp kernels across dimensions.
This model was included to demonstrate that the gain in predictive performance was not simply due to the inclusion of additive structure.
%
Hierarchical Kernel Learning (HKL)
%\footnote{Code for HKL available at \texttt{http://www.di.ens.fr/\textasciitilde fbach/hkl/}} 
was run using the all-subsets kernel, which corresponds to the same set of kernels as considered by the additive GP with a squared-exp base kernel.

Additive Gaussian Processes was run using SE as a base kernel, and at most ten degrees of interaction.

\paragraph{source code}All of the experiments in this paper were performed using the standard GPML toolbox\footnote{Available at \texttt{http://www.gaussianprocess.org/gpml/code/}}; code to perform all experiments is available on github.

%$k = 1$, $D = 8$, kernels are SE and RQ (\TBD{currently running other experiments that may be more canonical}).
%We have extended the comparison of \cite{duvenaud2011additive11} to include our method.

%Some points for discussion
%\begin{itemize}
%\item Experiments just using SE kernel can outperform additive kernel surprisingly. This is presumably a regularisation effect of using a finite depth search and/or BIC. We could make this a more Bayesian result (i.e.~more a property of the model) by placing a prior on kernels that depends on the number of components.
%\item Need to discuss design choices e.g.~$k$, depth of search, base kernels.
%\end{itemize}



\section{Discussion}

There is a large and mature literature on automatic structure discovery in unsupervised settings [citations including Josh].  However, relatively little attention has been paid to automatic structure discovery in supervised learning.

Similar searches over large model classes have been succesfully used in machine vision \cite{pinto2009high}
This work respresents a step towards more data-based modeling, as opposed to proposing the model beforehand.


Machine learning can be more data-driven, analogous to the high-thoughput approaches being used in biology. 

[Blessing of abstraction, doing lots with a small amount of data.]

[Extrapolation]

Choosing a kernel family is also a stumbling block for non-experts who wish to use Gaussian Process models.  Automating the search over kernel structures should be especially useful in allowing pracitioners from a wide variety of backgrounds to use appropriate models for their data.

%Indeed, the authors were sometimes surprised by the kernels chosen by the automated search, and the elegant ways in which complex structure was expressed through simple combination of kernels.



\section{Conclusion}

In this paper, we introduced an automated search over composite kernels, defining the search space through a simple grammar consisting of a set of one-dimensional base kernels and operations which combine those kernels.

We examined the composite kernels discovered by this search on several datasets, and demonstrate that they can recover known structure, and result in state-of-the-art predictive performance.  We further demonstrated an automated decomposition of the posterior into a sum of interpretable components.

%\section{Future Work}

%\TBD{RBG: Do we need this section?  Bill hates future work sections since they basically just give the reviewers a list of things to complain about.}

%\paragraph{Support Vector Machines} Structure search could be used to select kernels for support vector machines, using cross-validated prediction error as an objective function.

%However, it would be presumably quite difficult to learn the moderate number of kernel hyperparameters present in composite kernels using only cross-validation.

%\paragraph{Noise Models and Mean Functions}
%The models used in this paper all assumed Gaussian noise, as well a constant mean function for the Gaussian process.
%A more sophisticated approach might also define grammars over composite mean functions and noise models.
%However, in most cases, the covariance function can express structure that might be expressed in the mean function, and often in a more flexible way.
%The same can be said about noise models - for instance, a t-noise can be approximated by a sum of sq-exp kernels with very small lengthscales, but output variances which vary widely.


%\paragraph{Bayes model averaging}
%A natural extension of this work is to average over models in the grammar with a Bayes model averaging approach.
%This would mean augmenting our grammar with a prior over kernel structures, and approximately integrating over both structures and kernels.  

%In the first few levels of the search, where the difference between likelihoods of different kernel structures is usually high, model averaging might not make a substantial difference in terms of prediction.
%However, at deeper levels of the kernel search, integrating over structures and parameters would address overfitting concerns.

%\paragraph{Parsimony}
%Integrating over kernel structures raises the problem of reporting interpretable structures to the user.
%Currently, this is achieved through the use of the BIC complexity penalty, and reporting only the highest-scoring structure found.
%In experiments which used the more sophisticated Laplace approximation as a complexity penalty, the kernels found were more complicated.

%\paragraph{Bayesian Optimization} 
%The greedy search strategy used in this paper is rudimentary.  A natural extension of this work would be a model-based search over composite kernels, where the results of previous searches are used to propose informative directions.  This would be an example of ``meta-learning''.
%\cite{snoek2012practical} [also cite Lizotte, de Frietas]
%
%\paragraph{A Kernel between kernels}
%Placing a \gp{} model on the log-likelihood of different kernels would require constructing a kernel on kernels, as in \cite{ong2002hyperkernels}.

\subsubsection*{Acknowledgements}

\bibliographystyle{format/icml2013}
\bibliography{gpss}
\end{document}<commit_msg>Minor updates to text, and added solar decomposition fig<commit_after>\documentclass[twoside]{article}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{amssymb,amsmath,amsthm}
\usepackage{graphicx}
\usepackage{preamble}
\usepackage{natbib}
\usepackage{hyperref}
\usepackage{color}
\usepackage{wasysym}
\usepackage{subfigure}
\usepackage{tabularx}
\usepackage{booktabs}
\usepackage{bm}
\newcommand{\theHalgorithm}{\arabic{algorithm}}
\definecolor{mydarkblue}{rgb}{0,0.08,0.45}
\hypersetup{ %
    pdftitle={},
    pdfauthor={},
    pdfsubject={},
    pdfkeywords={},
    pdfborder=0 0 0,
    pdfpagemode=UseNone,
    colorlinks=true,
    linkcolor=mydarkblue,
    citecolor=mydarkblue,
    filecolor=mydarkblue,
    urlcolor=mydarkblue,
    pdfview=FitH}

\newcolumntype{x}[1]{>{\centering\arraybackslash\hspace{0pt}}m{#1}}
\newcommand{\tabbox}[1]{#1}


\setlength{\marginparwidth}{0.6in}
\input{include/commenting.tex}

\usepackage{format/icml2013}
%\usepackage[left=1.00in,right=1.00in,bottom=0.25in,top=0.25in]{geometry} %In case we want larger margins for commenting purposes

%% For submission, make all render blank.
%\renewcommand{\LATER}[1]{}
%\renewcommand{\fLATER}[1]{}
%\renewcommand{\TBD}[1]{}
%\renewcommand{\fTBD}[1]{}
%\renewcommand{\PROBLEM}[1]{}
%\renewcommand{\fPROBLEM}[1]{}
%\renewcommand{\NA}[1]{#1}  %% Note, NA's pass through!
    
\begin{document}



\twocolumn[
\icmltitle{Nonparametric Structure Discovery in Regression through Compositional Kernel Search}
%\icmlauthor{Anonymous Authors}
%\icmladdress{ Unknown Institution}
%%%%FIXME - Make sure to update this!
\icmlkeywords{nonparametrics, gaussian process, machine learning, ICML, structure learning, extrapolation, regression, kernel learning, equation learning, supervised learning}
\vskip 0.3in
]

\begin{abstract}
%Gaussian process (GP) models are used widely and successfully.
%However, their effictiveness depends critically on choosing an appropriate family of kernels.
%This aspect of GP modeling has been sorely underdeveloped.
%In this paper, we introduce a procedure for automatically and efficiently searching through a large space of GP models.
The effectiveness of nonparametric regression models depends heavily on the choice of kernel.
We introduce a marginal-likelihood-based search over composite kernel structures which automatically constructs a structured Gaussian process model appropriate for the dataset.
We further demonstrate that such kernels often allow the posterior to be automatically decomposed into a sum of interpretable components, and in some cases allows long-range extrapolation.
We demonstrate this technique on several real datasets, and achieve state-of-the-art predictive performance.
\end{abstract}

\section{Introduction}

Supervised learning problems, such as classification and regression, learn a function ${\function : \InputSpace \to \OutputSpace}$ from some input (predictor) variables, $\InputVar$, to some output (response) variables, $\outputVar$.
Kernel-based nonparametric models, such as support vector machines and Gaussian processes (\gp{}), have been one of the dominant paradigms for supervised machine learning over the last 20 years.
These methods depend on defining a kernel or covariance function, $\kernel : \InputSpace \times \InputSpace \to \Reals$ which specifies how similar or correlated outputs $\outputVar$ and $\outputVar'$ are at two inputs $\inputVar$ and $\inputVar'$, respectively.

A key challenge for kernel based methods is learning an appropriate kernel from data, and a great many papers have been written on this important topic\TBD{ cite kernel learning for SVMs and GP literature}.
In this paper, we pose the problem of kernel learning as a problem of structure discovery from data.
Specifically, we focus on a Bayesian setting where the kernel specifies a covariance function for Gaussian process (GP) regression.

Traditionally, Gaussian process kernels have simple forms, such as squared-exponential or Matern, which depend on few free hyperparameters which can be optimised (or inferred) from data.
However, one can also straightforwardly form composite kernels by combining simple base kernels through operations that such as summation, multiplication, and change of coordinates, which all preserve positive definiteness.
We use this insight to define a simple grammar over composite kernels, and we develop an automated search algorithm over the (exponentially large) space of kernels that can be derived from this grammar.
The search criterion is the marginal likelihood of the data for a given kernel, which is the driving term for Bayesian model comparison and selection in analogous structure discovery tasks such as graphical model learning \cite{heckerman1995learning}.

%\paragraph{Benefits of kernel structure}  In many regression problems, or time series analysis, or when modeling dynamical systems, we often think we are observing the superposition of multiple distinct causal processes.
%For example, an underlying linear function plus Gaussian noise in classical linear regression, or an underlying linear dynamics plus Gaussian observation noise in a Kalman filter.
%Our approach can be seen as a general way to learn functions like this, where the product of kernels can be used to capture richer causal proceses (e.g., nonlinear rather simply linear, amplifying perioidic (as in the second layer of the airline example) rather than simply periodic) and the sum of kernels corresponds to the superposition of these causal processes.

Our experimental findings demonstrate that this structure discovery algorithm can automatically recover known kernel structure from data, and that on real data the learned structure of the kernel often corresponds to natural, interesting and interpretable decompositions of the unknown function.
Our hope is that the algorithm developed in this paper will help replace the current and often opaque art of kernel engineering with a more transparent science of automated kernel discovery\fTBD{maybe $\to$ conclusion}.


\section{Expressing structure through kernels}
\label{sec:Structure}

Kernel functions $\kernel : \InputSpace \times \InputSpace \to \Reals$ can be used to define a measure of similarity between two points $\inputVar, \inputVar'$ in some space $\InputSpace$.
In the case of Gaussian process regression, the kernel defines a covariance function ${\KernelMatrix_{\inputVar \inputVar'} := \textrm{Cov}(\outputVar, \outputVar') = \kernel(\inputVar,\inputVar')}$.
%The kernel, $\kernel$, must define a valid covariance function\fTBD{expand me}; when this is the case $\kernel$ is said to be positive semi-definite (PSD).

Examples of commonly used kernels include squared exponential (SE), periodic (Per) and linear kernels (Lin) defined below\fTBD{If we add RQ then we have everything in one place}
\begin{eqnarray}
\kernel_\textrm{SE}(\inputVar, \inputVar') = & \sigma^2\exp\left(-\frac{(\inputVar - \inputVar')^2}{2\ell^2}\right) \\
\kernel_\textrm{Per}(\inputVar, \inputVar') = & \sigma^2\exp\left(-\frac{2\sin^2(\pi|\inputVar - \inputVar'|/p)}{\ell^2}\right) \\
\kernel_\textrm{Lin}(\inputVar, \inputVar') = & \sigma_1^2 + \sigma_2^2(\inputVar - \ell)(\inputVar' - \ell).
\end{eqnarray}

When used within the context of GP regression, these kernels define priors on functions which place mass only on smooth, periodic or linear functions respectively.
\TBD{Reference a diagram - even if a composite kernel we can point out what the base kernel is doing.}

Positive semidefinite kernels (\ie those which define valid covariance functions) are closed under addition and multiplication.
% \ie any algebraic composition of PSD kernels will define a PSD kernel.\fTBD{Cite theorem}
This allows one to create richly structured kernels from well understood base components.

\fTBD{Can probably kill simple kernel table. Focus on one table.}
%\input{tables/simple_kernels_table.tex}
\input{tables/example_structures_table.tex}

\paragraph{Summation}

Summation of kernels corresponds to the summation of independent functions in the following sense.
Suppose ${\function_1 \dist \GP(0, \kernel_1)}$, and ${\function_2 \dist \GP(0, \kernel_2)}$ and are independent.
Then ${\function := \function_1 + \function_2 \dist \GP(0, \kernel_1 + \kernel_2)}$\footnotemark.
\footnotetext{Additionally, the posterior of the component functions conditioned on observations of the sum are analytically tractable. We make use of this in section\TBD{ add reference} and give derivations in the supplementary material.}

In one dimension, sums of kernels can express structure such as variation at multiple scales, or of superposition of different function types.  Some examples can be found in Figure \ref{fig:kernels}.

In higher dimensions, summing kernels can correspond to decompositions of the form
\begin{equation}
\function((\inputVar_1, \inputVar_2)) = \function_1(\inputVar_1) + \function_2(\inputVar_2)
\end{equation}
\ie `independent' variation in different dimensions.
%
\input{tables/multi_d_additive_fig}
%
Figure \ref{fig:multi_d_additivity} demonstrates a decomposition across dimensions.

\paragraph{Interpretation}
\fTBD{Not sure about this section}
Loosely speaking, a sum of kernels can be understood as an 'OR' operation: two points are considered similar if either kernel has a high value.
More precisely, we can express a sum of kernels as a concatenation of their features in the corresponding RKHS.

\paragraph{Multiplication}

Multiplication of kernels does not correspond to a simple composition of functions; here, we give some examples of the properties that arise from multiplication of kernels.

%but allows for the construction of different structures and is best understood on a case by case basis\fTBD{Better language?}.
%
%For example, multiplying base kernels of the same form that act upon different dimensions of $\InputSpace$ can create the multidimensional generalisations of these kernels; this is true of \eg SE and LN\TBD{ reference a picture}.
%
%Multiplying kernels of the same form in one dimension is often redundant \eg for SE kernels this just produces SE kernels with different parameters.
%\TBD{What is PE $\times$ PE?}
%However, multiplying kernels of different forms can give rise to interesting structures.

Multiplying a kernel $\kernel$ by the SE kernel allows the structures implied by $\kernel$ to vary smoothly according the lengthscale of the SE kernel.  This is shown in row 2 of Figure \ref{kernels}.
Multiplying a kernel by LN results in the amplitude of structures implied by $\kernel$ growing or shrinking linearly. This is shown in rows 4 and 5 of Figure \ref{kernels}.
%\fTBD{Could perhaps perversely claim that multiplying by periodic results in wiggles being added to previous structure?}
%
%\fTBD{Can we do more than just give examples}

%In general, multiplying kernels can be thought of as an `AND' operation when viewing kernels as a measure of similarity.
%This arises since two points, $\inputVar,\inputVar'$, can only be considered similar by a product kernel function (\ie $\kernel_\textrm{prod}(\inputVar,\inputVar')$ is relatively large) if the points are considered similar by each component kernel function.

\paragraph{Interpretation}
Loosely speaking, multiplying kernels can be understood as an 'AND' operation: two points are considered similar only if both kernels have a high value.
More precisely, we can express a product of kernels as a concatenation of the cartesian product of both kernels' features in the corresponding RKHS.

To develop an intuition for the effect of multiplying two stationary kernels, we can note that the posterior mean function must be a weight sum of kernel functions.  Thus the third column of rows 1-2 of Figure \ref{kernels} show the functions that can be combined linearly to produce the posterior means in the 5th column.

%\section{Gaussian Processes Priors}

%Gaussian processes are a flexible and tractable prior over functions, useful for solving regression and classification tasks\cite{rasmussen38gaussian}.
%The kind of structure which can be captured by a GP model is mainly determined by its \emph{kernel}: the covariance function.
%One of the main difficulties in specifying a Gaussian process model is in choosing a kernel which can represent the structure present in the data.
%For small to medium-sized datasets, the kernel has a large impact on modeling efficacy.
%\TBD{Note: The above paragraph is plagarized from my additive GP paper. -David} 

%
%The technique of constructing composite kernels using sums and products of existing kernels is not new \cite{rasmussen38gaussian} [more cites, Phil Hennig's astronomy work?].  
%However, the main contribution of this paper is to automate the search over kernel structures.

\section{Searching over structures}

\paragraph{A grammar for kernels}

The discussion above demonstrates that rich structure can be captured with summation and multiplication of simple base kernels.
We therefore consider searching over all kernel structures that can be expressed as sums and products of base kernels.

%\begin{center}
%\begin{tabular}{rccc}
%\textrm{Replacement} & $\kernel_i$ & $\to$ & $\kernel'_i$\\% & $\forall\, \kernel' $\\
%\textrm{Addition} & $\kernel_i$ & $\to$ & $\kernel_i + \kernel'_j$\\% & $\forall\, j,\kernel' $\\
%\textrm{Multiplication} & $\kernel_i$ &  $\to$ & $\kernel_i \times \kernel'_j$\\% & $\forall\, j,\kernel'$\\
%\end{tabular}
%\end{center}

Details of the search algorithm are given in the supplementary material \TBD{For now}.

\paragraph{A greedy search algorithm}
In order to find the model with highest marginal likelihood, the base kernels and production rules generate a graph structure over kernel expressions and we follow a greedy search to traverse this graph.
This requires a method of scoring a kernel expression, which we achieve by using the Bayesian Information Criterion (BIC)\TBD{ cite me} as an estimate of the marginal likelihood of some data given the kernel expression\footnotemark.
\footnotetext{We also experimented with using a Laplace approximation to the marginal likelihood\TBD{ cite} but this was found to be less numerically stable and could not be computed if an optimiser failed to find a true optimum.}
Specifically, we optimise the marginal likelihood of the data given specific kernel parameters\TBD{ using LBFGS - cite?}, randomly restarting previously unoptimised kernel parameters.
This is then adjusted by the number of parameters optimised\footnotemark to compute the BIC\fTBD{formula?}.
\footnotetext{Adjustments were made for obvious parameter redundancy \eg multiplicative terms in products of kernels.}

The search then proceeds by scoring all base kernels applied to all input dimensions.
The highest scoring kernel is then expanded by applying all production rules, and these kernel expressions are scored.
The production rules are then applied to the highest scoring kernel expression and the search continues iteratively until some fixed depth.

\paragraph{Example expressions}

\TBD{Draw a picture to show what our model includes \ie polynomial regression, nonparametric `Fourier' decomposition, GAM, additive kernels.}

\section{Related Work}

\paragraph{Composite kernels in GP models} The technique of constructing composite kernels using sums and products of existing kernels was demonstrated in detail in Chapter 5 of \cite{rasmussen38gaussian}, where the resulting posterior mean was also decomposed into a sum of component-wise means, although the posterior variance was not.  Our work can be viewed as an extension and automation of this approach.

\NA(\paragraph{Generalized Additive Models} \cite{hastie1990generalized} are models for which the posterior mean can be expressed $\expect[f(\vx)] = \sum_{d=1}^D f_d(x_d)$.  ) These models have a limited compositional form, but one which is very interpretable.  We compare against this model in our multi-dimensional experiments.

\cite{plate1999accuracy} constructs a \gp{} with a composite kernel, summing an SE kernel along each dimension, with an SE-ARD kernel on all dimensions.  This model is motivated by the desire to specify the trade off between the interpretability of additive models and the flexibility of multiplicative models.

\paragraph{Compositional model search for unsupervised learning} \cite{grosse2012exploiting} performed a greedy search over a compositional model class for unsupervised learning.  Their models were composed of sums and element-wise products of random variables defined on matrices.  This model class contained a large number of existing unsupervised models as special cases.  Earlier work in unsupervised Bayseian structure discovery was done in \cite{kemp2008discovery}.

\paragraph{ANOVA Kernels}

Support vector regression with ANOVA decomposition kernels \cite{stitson1999support}

%\subsubsection{Smoothing spline ANOVA models}
A closely related procedure from the statistics literature is smoothing-splines ANOVA (SS-ANOVA)\cite{wahba1990spline, gu2002smoothing}.
This model is a weighted sum of splines along each dimension, plus a sum of splines over all pairs of dimensions, all triplets, etc, with each individual interaction term having a separate weighting parameter.
Learning in SS-ANOVA is usually done via penalized-maximum likelihood with a fixed sparsity hyperparameter.
Because the number of terms to consider grows exponentially in the order, in practice, only terms of first and second order are usually considered.  This model is similar to ours, if we were to restict our attention to spline kernels

\paragraph{Additive Gaussian Processes} \cite{duvenaud2011additive11} is a Gaussian process model whose kernel implicitly sums over all $2^D$ possible products of one-dimensional base kernels.  
However, this model's parameterization is not as flexible as the kernels expressible in our grammar.

\paragraph{Hierarchical Kernel Learning}

%In "High-Dimensional Non-Linear Variable Selection through Hierarchical Kernel Learning", Bach
\cite{DBLP:journals/corr/abs-0909-0844} uses a regularized optimization framework to learn a weighted sum over an exponential number of kernels which can be computed in polynomial time.
The subsets of kernels considered by this method are restricted to be a \textit{hull} of kernels. In the setting we are considering in this paper, a hull can be defined as a subset of all terms such that if term $\prod_{j \in J} k_j(\bf x, x')$ is included in the subset, then so are all terms $\prod_{j \in J / i} k_j(\bf x, x')$, for all $i \in J$.
%For details, see \cite{DBLP:journals/corr/abs-0909-0844}.
%Given each dimension's kernel, and a pre-defined weighting over all terms, HKL performs model selection by searching over hulls of interaction terms.
% \subsubsection{All-subsets kernel with uniform weightings}
%In \cite{DBLP:journals/corr/abs-0909-0844}, Bach also fixes the relative weighting between orders of interaction with a single term $\alpha$, computing the sum over all orders by:
%\begin{equation}
%\label{eqn:uniform}
%k_{a}({\bf x, x'}) = v_D^2 \prod_{d=1}^D \left(1 + \alpha k_{d}(x_{d}, x_{d}') \right)
%\end{equation}
%which has computational complexity $O(D)$.
%However, this formulation forces the weight of all $n$th order terms to be weighted by $\alpha^n$.
%
The main difficulty with this approach is that hyperparameters are hard to set other than by cross-validation. 
%In contrast, our method optimizes the hyperparameters of each dimension's base kernel, as well as the relative weighting of each order of interaction. 



%A related functional ANOVA GP model\cite{kaufman2010bayesian} decomposes the \emph{mean} function into a weighted sum of GPs.
%However, the effect of a particular degree of interaction cannot be quantified by that approach.
%Also, computationally, the Gibbs sampling approach used in \cite{kaufman2010bayesian} is disadvantageous.

\paragraph{Genetic Searches over Kernel Structures}
A similar method was developed by \cite{diosan2007evolving} using SVMs in place of Gaussian processes, cross-validation error instead of marginal likelihood, and a genetic algorithm instead of a greedy search.  
In \cite{bing2010gp}, a genetic programming construction and optimization method was used to select a kernel for relevance vector machines.

\paragraph{Equation Learning}
A related body of work with similar goals are the various approaches to equation learning.  
For example, \cite{dzeroski1999equation} and \cite{washio1999discovering} attempt to learn parametric forms of equations to describe time series, or relations between quantities.  The advantage of a \gp{} approach is that 

\paragraph{Semiparametric Regression}
Semiparametric regression \TBD{add cites} attempts to combine interpretability with flexibility by building  a composite model out of an interpretable, parametric part (such as linear regression) and a ``catch-all'' nonparametric part (such as a \gp{} with an SE kernel).
Structure search has the same goal, but improves upon semiparametric regression in two ways.
First, construcing a compound model through structure in the kernel allows us to integrate over the parameters of the implicit components of the model.
Second, by automatically searching over components, we remove the limitation that the practictioner must be able to specify in advance the parametric form of any part of the model.


\paragraph{Multiple Kernel Learning}
There is a large body of work attempting to construct a rich kernel through a weighted sum of base kernels, including in \gp{} models.  
\cite{christoudias2009bayesian}
\TBD{cite more} 
However, these methods typically construct all kernels in advance, and simply reweight them.  This sort of approach would be inapplicable in our case, since our search space is exponential in our search depth, and we also must estimate the hyperparameters of each kernel.%One could imagine expanding our search for some fixed depth, instantiating a sum of all
%In the \gp{} framework,  showed how mixtures of kernels can be learnt by gradient descent in the Gaussian process framework.
%They call this \emph{Bayesian localized multiple kernel learning}.

\paragraph{Nonparametric kernel functions}

Since the complexity of our kernel is bounded only by the search depth, we could consider the structure search to be a non-parametric kernel model - in a sense, a doubly nonparametric model.  However, there has been some similar work along these lines in \cite{salakhutdinov2008using}, where a deep neural network was use to learn the kernel of a \gp{}.

%Hyperkernels \cite{ong2002hyperkernels}

\section{Decomposing a function}

Sums of kernels were shown in section~\ref{sec:Structure} to correspond to sums of independent functions.
By splitting kernel expressions found by our search algorithm into additive components, we can also decompose the \gp{} posterior over functions into a joint distribution over functions that are summed togther.
This potentially allows us to probabilistically separate causally the independent processes that give rise to the data.

Moreover, the posterior distribution over additive components is analytically tractable, allowing us to plot the individual components.
%\fTBD{If kept here, make notation consistent through paper}
The conditional distribution of a Gaussian $\function_1$ conditioned on its sum with another Gaussian\fTBD{Make sure notation strictly correct and defined} $\vf = \vf_1 + \vf_2$ where $\vf_1 \sim \GP( \vmu_1, k_1)$ and $\vf_2 \sim \GP( \vmu_2, k_2)$ is given by
\begin{align}
\vf_1 | \vf \sim \mathcal{N} \big( & \vmu_1 + \vk_1\tra (\vK_1 + \vK_2)\inv \left( \vf - \vmu_1 - \vmu_2 \right), \nonumber \\
& \vk_1 - \vk_1\tra (\vK_1 + \vK_2)\inv \vk_1 \big).
\end{align}
%and the covariance between the two components, conditioned on their sum is given by:
%\begin{align}
%\cov(\vf_1^\star, \vf_2^\star) | \vf = \vk_1^{\star\tra} (\vK_1 + \vK_2)\inv \vk_2^\star
%\end{align}
%Derivations can be found in the supplementary material.


We demonstrate the decompositions produced by our method on two time series.
We ran the search to a depth of 10 and used sq-exp, periodic, linear and rational quadratic kernels.

% \input{"tables/kernel_descriptions.tex"}

%The sorts of structure implied by complex kernel might not be immediately interpretable.
%However, there is one property of GPs that allows us to break down the signal into interpretable parts.
%Specifically, a function $f \sim \GP(0, k_1(x, x') + k_2(x,x')$, whose kernel is a sum of kernels, has the same distribution as a sum of two functions $\vf = \vf_1 + \vf_2$ where $\vf_1 \sim \GP( \vmu_1, k_1)$ and $\vf_2 \sim \GP( \vmu_2, k_2)$.
%In addition, the posterior distributions over individual components after conditioning on data are also analytically tractable.

\paragraph{Mauna Loa Atmospheric Carbon Dioxide}

We revisited a dataset explored in \cite{rasmussen38gaussian}, pages 120-126, where a kernel was hand-tailored to fit a GP model to the dataset.

Figures \ref{fig:mauna_grow} show the posterior of the increasingly appropriate model as the search depth increases on this dataset.
Note in particular that, while the data can be smoothly interpolated by a single base kernel model, the extrapolations improve dramatically as the search depth increases.

Figure \ref{fig:mauna_decomp} shows the complete posterior of the maximum depth model together with the posteriors of the additive components.
The plot shows a very plausible extrapolation and highly interpretable components i.e.~a long-term trend, annual periodicity, medium-term deviations and short-term serially correlated `noise'.\fTBD{Remember to update}

\begin{figure*}[h!]
\centering
\newcommand{\wmg}{5.5cm}  % width maunu growth
\newcommand{\hmg}{4cm}  % height maunu growth
\begin{tabular}{c}
 \includegraphics[width=\wmg,height=\hmg]{../figures/decomposition/03-mauna2003-s_max_level_0/03-mauna2003-s_all_small} 
 \includegraphics[width=\wmg,height=\hmg]{../figures/decomposition/03-mauna2003-s_max_level_1/03-mauna2003-s_all_small}
 \includegraphics[width=\wmg,height=\hmg]{../figures/decomposition/03-mauna2003-s_max_level_2/03-mauna2003-s_all_small}
\end{tabular}
\caption{Posterior mean and variance for different depths of kernel search.  In the first row, the function is only modeled as a locally smooth function, and the extrapolation is extemely poor.  In the second, a periodic component is added, and the extrapolation is better.  At depth 3, the kernel can capture most of the relevant structure, and is able to extrapolate reasonably. \TBD{RBG: (1) I think we somehow need to visualize the lengthscales, to make it obvious that the SE kernels really mean different things. (2) Why isn't SE + PE the correct answer?}}
\label{fig:mauna_grow}
\end{figure*}

\begin{figure}[h!]
\newcommand{\wmgd}{9.5cm}  % width mauna decomp
\newcommand{\hmgd}{3.1cm}  % height mauna decomp
\newcommand{\mdrd}{../figures/decomposition/11-Feb-03-mauna2003-s}  % mauna decomp results dir
\begin{tabular}{c}
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_all} \\ = \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_1} \\ + \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_2_zoom} \\ + \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_3} \\ + \\
\hspace{-1cm} \includegraphics[width=\wmgd,height=\hmgd]{\mdrd/03-mauna2003-s_resid}
\end{tabular}
\caption{First row: The posterior on the Mauna Loa dataset after a search of depth 8.  Subsequent rows: automatic decomposition of Mauna Loa data.  The signal has been automatically decomposed into long-term, yearly periodic, medium-term anomaly and short-term noise components, respectively.}
\end{figure}
\label{fig:mauna_decomp}

\paragraph{Airline passenger data}

\TBD{Give an intro to this dataset.}
Figures \ref{fig:airline_grow} and \ref{fig:airline_decomp} show the corresponding plots for the airline data set (\TBD{describe me!}).
\TBD{Make similar observations of good performance and perhaps draw attention to the selection of a linear kernel}.

\begin{figure*}[h!]
\centering
\newcommand{\wag}{4.8cm}  % width airline growth
\newcommand{\hag}{4cm}  % height airline growth
\begin{tabular}{c}
\includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_0/01-airline-s_all_small}
\hspace{-1cm} \includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_1/01-airline-s_all_small}
\hspace{-1cm} \includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_2/01-airline-s_all_small} 
\hspace{-1cm} \includegraphics[width=\wag,height=\hag]{../figures/decomposition/01-airline-s_max_level_3/01-airline-s_all_small}
\end{tabular}
\caption{Posterior mean and variance for different levels of kernel search on the airline dataset.}
\label{fig:airline_grow}
\end{figure*}


\begin{figure}[h!]
\centering
\newcommand{\wagd}{9.5cm}  % width airline decomp
\newcommand{\hagd}{3.5cm}  % height airline decomp
%\newcommand{\ard}{../figures/decomposition/11-Feb-01-airline-s}  % airline results dir
\newcommand{\ard}{../figures/decomposition/31-Jan-v201-airline-months}  % airline results dir
\begin{tabular}{c}
\hspace{-1cm}\includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_all} \\
 = \\ 
%\hspace{-1cm} \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-s_4} \\
% + \\
\hspace{-1cm} \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_1} \\
 + \\
\hspace{-1cm}  \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_2} \\
 + \\
\hspace{-1cm}   \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_3} \\
 + \\
\hspace{-1cm} \includegraphics[width=\wagd,height=\hagd]{\ard/01-airline-months_resid}
\end{tabular}
\caption{First row:  The airline dataset and posterior after a search of depth 8.  Rows 2-4: Additive decomposition of posterior into long-term smooth trend, yearly variation, and short-term noise.  Note that the variance of the noise grows over time, making this a heteroskedastic model.}
\label{fig:airline_decomp}
\end{figure}

%To evaluate the effectiveness of our method we performed two types of experiments.
%First, we examined the ability of our algorithm to discover useful structure in one-dimensional datasets.  Second, we examined the predictive performance of our model on multi-dimensional datasets.
% testing the ability of the algorithm to infer both an accurate regression function estimate and to produce an interpretable decomposition of a regression function.
%The accuracy experiments show that our method consistently matches or beats previous state of the art regression methods (\TBD{not yet statistically significantly within each experiment - more experiments currently running\ldots}).
%The decomposition experiments produce highly interpretable decompositions of time series data and highly plausible extrapolations; a particularly rare property of naively applied (linear) smoothers (\TBD{justify comment - e.g.~local linear is just linear, GP Sq-Exp nose dives towards the mean}).

\paragraph{Solar Irradiance Data} 
\cite{lean1995reconstruction} 

\begin{figure}[h!]
\newcommand{\wsd}{9.5cm}  % width solar decomp
\newcommand{\hsd}{3.1cm}  % height solar decomp
\newcommand{\srd}{../figures/decomposition/11-Feb-02-solar-s}  % solar decomp results dir
\begin{tabular}{c}
\hspace{-1cm} \includegraphics[width=\wsd,height=\hsd]{\srd/02-solar-s_all} \\ = \\
\hspace{-1cm} \includegraphics[width=\wsd,height=\hsd]{\srd/02-solar-s_1} \\ + \\
\hspace{-1cm} \includegraphics[width=\wsd,height=\hsd]{\srd/02-solar-s_2} \\ + \\
\hspace{-1cm} \includegraphics[width=\wsd,height=\hsd]{\srd/02-solar-s_resid}
\end{tabular}
\caption{First row: The posterior on the solar irradiance dataset after a search of depth 10.  Subsequent rows: automatic decomposition.}
\label{fig:solar_decomp}
\end{figure}


Figure \ref{fig:solar_decomp} shows a failure mode of the model. \TBD{Write more here}

\subsection{Multidimensional decomposition}

\begin{figure*}[h]
\centering
\begin{tabular}{ccc}
\includegraphics[width=0.3\textwidth]{../figures/additive_multi_d/interpretable_1st_order1.pdf} &
\includegraphics[width=0.3\textwidth]{../figures/additive_multi_d/interpretable_1st_order2.pdf}& 
\includegraphics[width=0.3\textwidth]{../figures/additive_multi_d/interpretable_2nd_order1.pdf}\\
\end{tabular}
\caption{One-dimensional decompositions of the concrete dataset posterior.  Left, Centre:  Green points indicate the original data, blue points are data after the mean contribution from the other dimensions' first-order terms has been subtracted.  The black line is the posterior mean of a GP with only one term in its kernel.  Right:  The posterior mean of a GP with only one second-order term in its kernel.}
\label{fig:interpretable functions}
\end{figure*}
%
Figure \ref{fig:interpretable functions} demonstrates that an additive posterior can be decomposed into a sum of functions across dimensions.


\section{Empirical evaluation}

In addition to producing highly interpretable decompositions of regression functions, our method also produces state of the art predictions in both extrapolation and interpolation tasks.

\subsection{Extrapolation}

\begin{figure}[h!]
\centering
\begin{tabular}{c}
\hspace{-0.5cm} \includegraphics[width=\columnwidth]{../figures/extrapolation_curves/01-airline-s-ex-curve.pdf}
\end{tabular}
\caption{Extrapolation performance.  Here, we show the average test set MSE as a function of the fraction of the dataset used for training.}
\end{figure}

Figure \ref{extrapolation} shows the learning curves of our method compared against simple linear regression and \gp{} regression with a SE kernel.  The amount of data in the 1D airline timeseries was increased by increments of 10\%. Both linear and SE regression plateau, while unbounded complexity of our method allows the model to continue to capture structure. 

\subsection{Interpolation}

To demonstrate the interpolation accuracy of our method we extend the comparison of \cite{duvenaud2011additive11} to include our method.
5 comparsion algorithms were evaluated on 5 datasets in terms of mean-squared (MSE) error and predictive likelihood using 10 fold cross validation.
Results are presented in tables \ref{tbl:Regression Mean Squared Error} and \ref{tbl:Regression Negative Log Likelihood}.
%
\input{tables/regression_results_mse.tex}
\input{tables/regression_results_nll.tex}
%
Our method outperforms the other methods in all tests.

\TBD{Other versions of tests / learning curves}

\subsubsection{Data sets}


\paragraph{Bach Synthetic Dataset}
\fTBD{Too much text? Move to supplementary}
In addition to standard UCI repository datasets, we generated a synthetic dataset following the same recipe as \cite{DBLP:journals/corr/abs-0909-0844}.
% From a covariance matrix drawn from a Wishart distribution with 1024 degrees of freedom, we select 8 variables.
%We then construct the non-linear function $f(X) = \sum_{i=1}^4 \sum_{j=1+1}^4 X_i X_j + \epsilon$, which sums all 2-way products of the first 4 variables, and adds Gaussian noise $\epsilon$.
This dataset is one which can be predicted well by a kernel which is a sum of two-way interactions over the first 4 variables, ignoring the extra 4 noisy copies.
%
This dataset was designed by \cite{DBLP:journals/corr/abs-0909-0844} to demonstrate the advantages of HKL over a \gp{} with a SE-ARD kernel. 

%If the dataset is large enough, HKL can construct a hull around only those subsets of cross terms that are optimal for predicting the output.
%GP-ARD, in contrast, can only learn to ignore the noisy copy variables, but cannot learn to ignore the higher-term interactions between the predictive variables.
%However, a GP with an additive kernel can learn both to ignore irrelevant variables, and to ignore certain orders of interaction.
%In this example, the additive GP is able to recover the relevant structure.

[TODO: Describe concrete, pumadyn, servo and housing]

\subsection{Details of methods}

Linear regression was included as a baseline, with noise variance set to the empirical variance on the training residuals.
%
Generalized Additive Models (GP GAM) were implemented through a \gp{} whose kernel is a sum of sq-exp kernels across dimensions.
This model was included to demonstrate that the gain in predictive performance was not simply due to the inclusion of additive structure.
%
Hierarchical Kernel Learning (HKL)
%\footnote{Code for HKL available at \texttt{http://www.di.ens.fr/\textasciitilde fbach/hkl/}} 
was run using the all-subsets kernel, which corresponds to the same set of kernels as considered by the additive GP with a squared-exp base kernel.

Additive Gaussian Processes was run using SE as a base kernel, and at most ten degrees of interaction.

\paragraph{source code}All of the experiments in this paper were performed using the standard GPML toolbox\footnote{Available at \texttt{http://www.gaussianprocess.org/gpml/code/}}; code to perform all experiments is available on github.

%$k = 1$, $D = 8$, kernels are SE and RQ (\TBD{currently running other experiments that may be more canonical}).
%We have extended the comparison of \cite{duvenaud2011additive11} to include our method.

%Some points for discussion
%\begin{itemize}
%\item Experiments just using SE kernel can outperform additive kernel surprisingly. This is presumably a regularisation effect of using a finite depth search and/or BIC. We could make this a more Bayesian result (i.e.~more a property of the model) by placing a prior on kernels that depends on the number of components.
%\item Need to discuss design choices e.g.~$k$, depth of search, base kernels.
%\end{itemize}



\section{Discussion}

There is a large and mature literature on automatic structure discovery in unsupervised settings [citations including Josh].  However, relatively little attention has been paid to automatic structure discovery in supervised learning.

Similar searches over large model classes have been succesfully used in machine vision \cite{pinto2009high}
This work respresents a step towards more data-based modeling, as opposed to proposing the model beforehand.


Machine learning can be more data-driven, analogous to the high-thoughput approaches being used in biology. 

[Blessing of abstraction, doing lots with a small amount of data.]

[Extrapolation]

Choosing a kernel family is also a stumbling block for non-experts who wish to use Gaussian Process models.  Automating the search over kernel structures should be especially useful in allowing pracitioners from a wide variety of backgrounds to use appropriate models for their data.

%Indeed, the authors were sometimes surprised by the kernels chosen by the automated search, and the elegant ways in which complex structure was expressed through simple combination of kernels.



\section{Conclusion}

In this paper, we introduced an automated search over composite kernels, defining the search space through a simple grammar consisting of a set of one-dimensional base kernels and operations which combine those kernels.

We examined the composite kernels discovered by this search on several datasets, and demonstrate that they can recover known structure, and result in state-of-the-art predictive performance.  We further demonstrated an automated decomposition of the posterior into a sum of interpretable components.

%\section{Future Work}

%\TBD{RBG: Do we need this section?  Bill hates future work sections since they basically just give the reviewers a list of things to complain about.}

%\paragraph{Support Vector Machines} Structure search could be used to select kernels for support vector machines, using cross-validated prediction error as an objective function.

%However, it would be presumably quite difficult to learn the moderate number of kernel hyperparameters present in composite kernels using only cross-validation.

%\paragraph{Noise Models and Mean Functions}
%The models used in this paper all assumed Gaussian noise, as well a constant mean function for the Gaussian process.
%A more sophisticated approach might also define grammars over composite mean functions and noise models.
%However, in most cases, the covariance function can express structure that might be expressed in the mean function, and often in a more flexible way.
%The same can be said about noise models - for instance, a t-noise can be approximated by a sum of sq-exp kernels with very small lengthscales, but output variances which vary widely.


%\paragraph{Bayes model averaging}
%A natural extension of this work is to average over models in the grammar with a Bayes model averaging approach.
%This would mean augmenting our grammar with a prior over kernel structures, and approximately integrating over both structures and kernels.  

%In the first few levels of the search, where the difference between likelihoods of different kernel structures is usually high, model averaging might not make a substantial difference in terms of prediction.
%However, at deeper levels of the kernel search, integrating over structures and parameters would address overfitting concerns.

%\paragraph{Parsimony}
%Integrating over kernel structures raises the problem of reporting interpretable structures to the user.
%Currently, this is achieved through the use of the BIC complexity penalty, and reporting only the highest-scoring structure found.
%In experiments which used the more sophisticated Laplace approximation as a complexity penalty, the kernels found were more complicated.

%\paragraph{Bayesian Optimization} 
%The greedy search strategy used in this paper is rudimentary.  A natural extension of this work would be a model-based search over composite kernels, where the results of previous searches are used to propose informative directions.  This would be an example of ``meta-learning''.
%\cite{snoek2012practical} [also cite Lizotte, de Frietas]
%
%\paragraph{A Kernel between kernels}
%Placing a \gp{} model on the log-likelihood of different kernels would require constructing a kernel on kernels, as in \cite{ong2002hyperkernels}.

\subsubsection*{Acknowledgements}

\bibliographystyle{format/icml2013}
\bibliography{gpss}
\end{document}



% Stuff that didn't make the cut

%Kernel-based nonparametric regression models have been widely and succesfully used over the last 20 years. [citations needed]
%As part of the development of this model class, a rich set of kernels have been developed for both continuous and structured data.
%Because different kernels reflect different properties of the function being modeled, the choice of kernel can have a large impact on the sorts of structure that can be expressed or discovered by the model.

%When approaching a new dataset, it is often not clear which kernel family is most appropriate, especially for non-experts.
%Often, a standard kernel family (typically squared-exponential) is chosen out of convenience, and kernel parameters are set by cross-validation or maximum marginal likelihood. [citations?]

%Sometimes the argument is made that since most kernels are universal (i.e. consistent in the limit), that the choice of kernel is not important[citation?].  However, our experiments demonstrate that for datasets of moderate size, predictive performance, extrapolation and interpretability all depend heavily on the choice of kernel. \TBD{RBG: Do people actually say the choice of kernel isn't important?  This sounds like a straw man...}

%\paragraph{Main contribution}
%In this paper, we introduce an automated search over composite kernels.
%We define the search space through a simple grammar consisting of a set of one-dimensional base kernels (squared-exp, periodic, linear, etc.) and operations which combine those kernels (add and multiply).
%Our objective function is the marginal likelihood of a Gaussian process model with the given kernel, conditioned on the data.

%We examine the composite kernels discovered by this search on several datasets, and demonstrate that they can recover known structure, and result in state-of-the-art predictive performance.  We further demonstrate that in some cases, the structured posterior can be automatically decomposed into a sum of interpretable components, e.g. in Figure \ref{fig:airline_decomp}. \TBD{RBG: Need to clarify that the contribution isn't that we can decompose it (that was in GPML), but that the decomposition is often interesting.}

%\section{Motivation}
%\TBD{RBG: Do we need a separate motivation section?  The introduction should motivate the work.}

%There is a large and mature literature on automatic structure discovery in unsupervised settings [citations including Josh].
%However, relatively little attention has been paid to automatic structure discovery in supervised learning.

%Similar searches over large model classes have been succesfully used in machine vision [cite Cox + Pinto].
%  This work respresents a step towards more data-based modeling, as opposed to proposing the model beforehand.
%In general, learning the model class from data seems superior proposing the model beforehand.
%In high dimensional problems, it is also hard for a practitioner to propose an appropriate model even after examining a dataset closely.
% Choosing a kernel family is also a stumbling block for non-experts who wish to use Gaussian Process models.

%Machine learning can be more data-driven, analogous to the high-thoughput approaches being used in biology. 

%\subsection{Convenience Priors}
%One of the main advantage of kernel methods, and Gaussian process regression in particular, is that the kernel can be used to specify interesting prior structure to the model.
%However, in practice, GPs are typically used with only the squared-exp kernel.
%[Apparently, this was one of Sholkopf's big dissapointment with GPs... is there a quote from him on this? -David]

%A frequent criticism of the Bayesian framework is that priors are often chosen purely for convenience. 
%Gaussian process priors may be a good examples of a model chosen for its tractability.
%However, specifying a compound kernel typically does not change the order of computational complexity of inference in Gaussian process models.
%The broad use of simple kernels may be simply due to the difficulty of determining which kinds of structure exist in the data.
%As well, even if one understands the structure present in a dataset, it is often non-obvious how to express that through a covariance function. 
%\TBD{RBG: Not sure we need this in the motivation. Is this basically just to answer the objection of why we don't consider the tradeoff between expressiveness and tractability?  Also, the prior can affect computational tractability if we're using approximate algorithms, so this makes it sound like we're only thinking about small datasets.}

%Automating the search over kernel structures should be especially useful in allowing pracitioners from a wide variety of backgrounds to use appropriate models for their data.
%Indeed, the authors were sometimes surprised by the kernels chosen by the automated search, and the elegant ways in which complex structure was expressed through simple combination of kernels.

{'fraction_non_alphanumeric': 0.07116848725587553, 'fraction_numerical': 0.010656586688876786, 'mean_word_length': 5.0612258106271435, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6763577', 'n_tokens_mistral': 25386, 'n_tokens_neox': 23471, 'n_words': 11857}