Evolution of Covariance Functions for Gaussian
Process Regression using Genetic Programming
Gabriel Kronberger and Michael Kommenda

arXiv:1305.3794v2 [] 22 May 2013

School of Informatics, Communications and Media,
University of Applied Sciences Upper Austria,
Softwarepark 11, 4232, Hagenberg
{gabriel.kronberger,michael.kommenda}@fh-hagenberg.at

Abstract. In this contribution we describe an approach to evolve composite covariance functions for Gaussian processes using genetic programming. A critical aspect of Gaussian processes and similar kernel-based
models such as SVM is, that the covariance function should be adapted to
the modeled data. Frequently, the squared exponential covariance function is used as a default. However, this can lead to a misspecified model,
which does not fit the data well.
In the proposed approach we use a grammar for the composition of covariance functions and genetic programming to search over the space of
sentences that can be derived from the grammar.
We tested the proposed approach on synthetic data from two-dimensional
test functions, and on the Mauna Loa CO2 time series. The results show,
that our approach is feasible, finding covariance functions that perform
much better than a default covariance function. For the CO2 data set a
composite covariance function is found, that matches the performance of
a hand-tuned covariance function.

Keywords: Gaussian Process, Genetic Programming, Structure Identification

1

Introduction

The composition of covariance functions is a non-trivial task and has been described as a black art [2]. On the one hand, it is critical to tune the covariance
function to the data set, that should be modeled, because this is the primary
option to integrate prior knowledge into the learning process [10]; on the other
hand a lot of experience and knowledge about the modeled system is required
to do this correctly. Frequently, and especially for multi-dimensional data sets it
is far from obvious how the covariance function should be structured.
In this work we discuss the composition of covariance functions for Gaussian
processes, that can be used for nonparametric machine learning tasks e.g., for
regression or classification [10]. In this context a Gaussian process is used as a
Bayesian prior over functions, relating the input variables to the target variable.
Gaussian process regression allows modeling of non-linear functional dependencies through different covariance functions, and produces posterior probability
distribution estimates for the target values instead of point estimates only.

2

1.1

Evolution of Covariance Functions for Gaussian Process Regression

Our Contribution

The aim of this paper is to describe the idea of using a grammar for covariance
functions and genetic programming to search for a good covariance function for a
given data set. We also describe our prototype implementation using grammarguided tree-based GP, and finally, present results as a proof-of concept. We
have not yet evaluated the difficulty of this problem for genetic programming,
and in particular, if GP suited well for this kind of problem. The results of
our experiments indicate that the idea is feasible, producing good covariance
functions for low-dimensional data sets.
1.2

Previous Work

In a very recent contribution, the problem of structure identification for covariance functions has been approached using a grammar, or rather a set of rewriting
rules, as a basis for searching over composite covariance functions for Gaussian
processes [2]. This approach is actually very similar to our work; the main difference is, that in our work we use genetic programming to search over the set of
possible structures, while Duvenaud et al. enumerate over composite functions,
starting with standard functions.
Another recent contribution discusses more flexible families of covariance
functions, instead of composing covariance functions from simple terms [12].
Also related is earlier work that describes additive Gaussian processes [3], which
are equivalent to a weighted additive composition of base kernels, but can be
calculated efficiently.
Genetic programming has been used previously to evolve kernel functions for
SVMs with mixed results [4], [6]. The latest contribution found that genetic programming was able to “rediscover multiple standard kernels, but no significant
improvements over standard kernels were obtained” [7]. These results can, however, not be transfered directly to Gaussian processes because of several major
differences between Gaussian processes and SVMs. In particular, in the case of
Gaussian processes hyper-parameters are optimized using a ML-II approach, in
contrast to SVMs, where hyper-parameter values are usually tuned using crossvalidation and grid-search. Additionally, in contrast to all other previous work,
simple embeddings of covariance functions by masking dimensions are supported.

2

Gaussian Processes

A Gaussian process is a non-parametric model that produces predictions solely
from the specified mean and covariance functions and the available training
data [10]. The inference of function values f ∗ for observed input values X ∗ based
on observations of y and X involves the calculation of the covariance matrices K(X, X) and K(X, X ∗ ) and inference from the multi-dimensional Gaussian
shown in Equation 1.
 

 

y
m(X)
K(X, X) + σ 2 I K(X, X ∗)
∼
N
,
(1)
f∗
m(X ∗ )
K(X ∗ , X) K(X ∗ , X ∗ )

Evolution of Covariance Functions for Gaussian Process Regression

3

The term σ 2 I is necessary to account for Gaussian distributed noise with
variance σ 2 . From this definition it follows that the posterior for f ∗ is again
a multi-dimensional Gaussian. For model selection and hyper-parameter learning the marginal likelihood p(y|X) must be calculated. The model is a multidimensional Gaussian so an analytical form of the likelihood can be derived.
Calculation of the marginal likelihood requires a matrix inversion and, thus,
has asymptotic complexity O(n3 ). Usually, the covariance function K(x, x′ ) has
hyper-parameters that must optimized. This is often accomplished in a simple ML-II fashion, optimizing the hyper-parameters w.r.t. the likelihood using
a quasi-Newton method (e.g., BFGS). Since the gradients of the marginal likelihood for the hyper-parameters can be determined with an additional computational complexity of O(n2 ) for each hyper-parameter, it is feasible to use
gradient-based methods. The drawback is that the likelihood is typically multimodal, and especially for covariances with many hyper-parameters (e.g., ARD)
the optimizer can converge to a local optimum. Thus, it is typically suggested
to execute several random restarts. A better solution would be to include priors
on the hyper-parameters and optimizing w.r.t. posterior distribution (MAP).
However, this can only be accomplished using a MCMC approach which is computationally expensive.
Frequently used covariance functions for Gaussian processes include the linear, polynomial, squared exponential (SE), rational quadratic (RQ) and the
Matérn function. Covariance functions can be combined to more complex covariance functions, for instance as products or sums of different covariance functions [10].

3

Genetic Programming

Genetic programming generally refers to the automatic creation of computer
programs using genetic algorithms [8]. The basic principle is to evolve variablelength structures, frequently symbolic expression trees, which represent potential
solutions to the problem. One of the most prominent applications of genetic programming is symbolic regression, the synthesis of regression models without a
predetermined structure. Genetic programming makes it possible to optimize the
structure of solutions in combination with their parameters. Thus, it should also
be possible to synthesize composite covariance functions with genetic programming. In the following, we use a grammar-guided genetic programming system
to make sure that only valid covariance functions are produced. A good survey
of grammar-guided genetic programming is given in [9].

4

Grammar for Covariance Functions

The grammar for covariance functions has been derived from the rules for the
composition of kernels as e.g., discussed in [10]. It should be noted that the grammar shown below is not complete, meaning that several constructions that would

4

Evolution of Covariance Functions for Gaussian Process Regression

lead to a valid covariance function are not possible1 . The following represents
the grammar G(Cov) for covariance functions in EBNF notation2 :
Cov

-> "Prod" "(" Cov { Cov } ")" | "Sum" "(" Cov { Cov } ")" |
"Scale" Cov | "Mask" BitVector Cov | TerminalCov .
TerminalCov -> "SE" | "RQ" | "Matern1" | "Matern3" | "Matern5" |
"Periodic" | "Linear" | "Constant" | "Noise" .
BitVector
-> "[" {"0" | "1" } "]" .

The functions Prod and Sum produce the product and sum of multiple covariance functions, which can again be composite covariance functions. The scale
operator can be used to add a scaling factor to any covariance function. The Mask
operator selects a potentially empty subset of input variables from all possible
input variables. The non-terminal symbol BitVector can be derived to a list
of zeros and ones. The bit vector is used to mask selected dimensions in the
data set, effectively reducing the dimensionality. The length of the bit mask has
to match the total number of dimensions; this is checked when the resulting
covariance function is evaluated.
Finally, the non-terminal symbol TerminalCov can be derived to a range of
default covariance functions. Currently, we only included isometric covariance
functions, but other covariance functions can be added to the grammar easily.
The grammar does not include the hyper-parameters, because they are not optimized by genetic programming. Instead, hyper-parameters are optimized for
each potential solution, using a gradient-descent technique.

5

Experiments

For the experiments we implemented Gaussian processes, a set of commonly
used covariance functions, and the grammar for covariance functions in HeuristicLab3 [11] which already provides an implementation of grammar-guided treebased genetic programming.
The aim of the experiments presented in this contribution is mainly to test
the feasibility of the idea. Two different types of data sets are used for the experiments, and the forecasts of the synthesized covariance functions are compared
to a set of default covariance functions and also to hand-tuned covariance functions. The first data set is the univariate Mauna Loa atmospheric CO2 time
series. This data set has been chosen, because a hand-tuned covariance function
for this data set is presented in [10]. For the second experiment we created several
synthetic data sets sampled randomly from two-dimensional Gaussian process
priors shown in Equation 2. The data generated from these functions are difficult to model with a single isometric covariance function. Multiple covariance
functions have to be combined and the correct dimension masking vectors have
1
2

3

One example is vertical scaling of covariance functions: K ′ (x, x′ ) = a(x)K(x, x′ )a(x′ )
The grammar is largely based on the capabilities of the GPML package by Rasmussen
and Nickisch, http://gaussianprocess.org/gpml/code.
HeuristicLab version 3.3.8 is available from http://dev.heuristiclab.com/

Evolution of Covariance Functions for Gaussian Process Regression

5

to be identified. Each data set contains 882 samples of the function on a regular
two-dimensional grid.
SE+RQ(x, x′ ) = SE(x0 , x′0 ) + RQ(x1 , x′1 )
SE+Matérn(x, x′ ) = SE(x0 , x′0 ) + Matérn1(x1 , x′1 )
SE+Periodic(x, x′ ) = SE(x0 , x′0 ) + RQ(x1 , x′1 )
5.1

(2)

Genetic Programming Parameter Settings

Training of Gaussian processes is computationally expensive, and because it is
necessary to optimize the hyper-parameters for each evaluated covariance function the run time of the genetic programming algorithm grows quickly. Therefore,
we used very restrictive parameter settings, in particular a small population size
of only 50 individuals. All other parameter settings are shown in Table 1.
Table 1. Genetic programming parameter settings for all experiments.
Parameter
Population size
Max. length / height
Initialization
Parent selection
Mutation rate
ML-II iterations
Offspring selection [1]
Max. selection pressure
Max. generations

5.2

Value
50
25 / 7
PTC2
gender-specific (proportional + random )
15%
50
strict (success ratio = 1, comparison factor = 1)
100
20

Results on Mauna Loa CO2 Data Set

The results for the CO2 time series are positive. The algorithm was able to
consistently find covariance functions that fit well in the training period (1958 –
2004), accurate forecasts over the testing period (2004 – 2012). The structures
of two exemplary solutions are shown in Equation 3. The first solution (K1)
is actually very similar to the hand-tuned covariance solution proposed in [10].
The second covariance function is more complex and has only a slightly better
likelihood. Unfortunately, genetic programming often leads to overly complex
solutions which is a critical drawback of our approach. Both solutions have been
found after only 800 evaluated solution candidates and achieve a negative loglikelihood of 129.8 and 116, respectively. The correlation coefficients for the
forecasts in the test partition are above 0.99. Figure 1 shows the output of the
first model.

6

Evolution of Covariance Functions for Gaussian Process Regression

K1(x, x′ ) = SE(x, x′ ) + Periodic(x, x′ ) + Matérn1(x, x′ )+
SE(x, x′ ) + Matérn5(x, x′ ) + Const
K2(x, x′ ) = Matérn3(x, x′ ) ∗ Perioric(x, x′ ) ∗ RQ(x, x′ ) ∗
(Matérn1(x, x′ ) + Matérn3(x, x′ ) + Matérn5(x, x′ ) +

(3)

Perioric(x, x′ ) + Linear(x, x′ )) ∗
(Matérn1(x, x′ ) + Matérn3(x, x′ ) + RQ(x, x′ ))

Estimated Values (training)

Estimated Values (test)

80
60
40
20
0
-20
-40
0

149

298

447

596

745

Fig. 1. The output and forecast for the Mauna Loa CO2 time series of a Gaussian
process using the first evolved covariance function (K1) shown in Equation 3.

5.3

Results on Synthetic Data Sets

The results for the synthetic two-dimensional data set are shown in Table 2.
In this experiment we trained multiple Gaussian process models using several
frequently used covariance functions. We trained many models using random
restarts for each data set and covariance function, and report the best negative log-likelihood for each pair. As expected, the models with the isometric
covariance functions do not fit well. In contrast, the composite covariance functions produced by genetic programming fit much better. For comparison, we also
report the negative log-likelihood, that can be achieved with the optimal covariance function for each data set. In these experiments, the exact structure of the
covariance could not be rediscovered, thus, the evolved functions are worse than
the optimal solution.

6

Summary and Discussion

In this contribution we described an approach for the synthesis of composite covariance function for Gaussian processes using grammar-guided genetic programming. In the proposed approach a set of commonly used covariance functions is

Evolution of Covariance Functions for Gaussian Process Regression

7

Covariance

Table 2. Best negative log-likelihood achieved for the three synthetic two-dimensional
test functions, with default covariance functions and with evolved composite covariance
functions.

SE
RQ
Periodic
Matérn
Evolved
Optimal

Problem instance
SE+RQ SE+Matérn SE+Periodic
-204
-492
440
-272
-492
103
-221
-492
479
-27
31
304
-803
-760
-640
-2180
-2187
-2131

used to compose more complex covariance functions, using sums or products of
several covariance functions. The set of valid covariance functions is defined via a
grammar and genetic programming is used to search the space of possible derivations from this grammar. The hyper-parameters of covariance functions are not
subject to the evolutionary search, but are optimized w.r.t. the likelihood using
a standard gradient-descent optimizer (i.e., LBFGS).
The proposed approach was tested on two types of low-dimensional problems
as a proof of concept. We found, that for the univariate Mauna Loa CO2 time
series it is possible to consistently find good covariance functions with genetic
programming. The identified solutions perform as well as a hand-tuned covariance function for this problem. The results for our two-dimensional synthetic
functions show that it is possible to find composite covariance functions, which
perform much better than default covariance functions on these data sets.
In contrast to previous work by the genetic programming community [7],
which focused mainly on kernel synthesis for SVMs, this contribution discusses
kernel synthesis for Gaussian processes, which are non-parametric fully Bayesian
models. For Gaussian process models the hyper-parameters can be optimized
with a standard gradient-descent approach, and it is not strictly necessary to
execute cross-validation [10]. Previous work either used grid-search and crossvalidation to tune hyper-parameters, which is very computationally expensive, or
did not consider hyper-parameter optimization at all. Additionally, we are using
a grammar to compose covariance functions from simple covariance functions
instead of evolving the full function.
In the statistics community, a very recent contribution has also discussed
the usage of grammars for the composition of covariance functions [2]. The main
difference to this work is that here genetic programming is used to search over the
derivations of the grammar. Another relevant difference is that the grammar in
this contribution also supports simple embeddings through the masking function.
It should be noted, that we have not yet analyzed if genetic programming is well
suited for this task, and in particular we did not compare the approach to simple
enumeration or random search.
One question that remains for future work is whether composed covariance
functions also work well for data sets with more variables. We have observed that

8

Evolution of Covariance Functions for Gaussian Process Regression

simple covariance functions often work very well, and tuned covariance functions
do not have a strong beneficial effect for these data sets.
Another interesting topic for future research is to look at alternative ways for
searching over the space of covariance functions defined by a grammar. Recently,
an interesting approach has been described that uses variational methods for
Bayesian learning of probabilistic context free grammars for this task [5]. This
idea could be especially useful for Bayesian models such as Gaussian processes.
Acknowledgments The authors would like to thank Jeffrey Emanuel for the
initial idea leading to this contribution. This work has been supported by the
Austrian Research Promotion Agency (FFG), on behalf of the Austrian Federal
Ministry of Economy, Family and Youth (BMWFJ), within the program “Josef
Ressel-Centers”.

References
1. Affenzeller, M., Winkler, S., Wagner, S., Beham, A.: Genetic Algorithms and Genetic Programming: Modern Concepts and Practical Applications. Numerical Insights, CRC Press, Singapore (2009)
2. Duvenaud, D., Lloyd, J.R., Grosse, R., Tenenbaum, J.B., Ghahramani, Z.: Structure Discovery in Nonparametric Regression through Compositional Kernel Search.
ArXiv e-prints (Feb 2013)
3. Duvenaud, D., Nickisch, H., Rasmussen, C.E.: Additive Gaussian processes. arXiv
preprint arXiv:1112.4394 (2011)
4. Gagne, C., Schoenauer, M., Sebag, M., Tomassini, M.: Genetic programming for
kernel-based learning with co-evolving subsets selection. In: Proc. of PPSN IX.
LNCS, vol. 4193, pp. 1008–1017. Springer-Verlag (2006)
5. Hasegawa, Y., Iba, H.: Latent variable model for estimation of distribution algorithm based on a probabilistic context-free grammar. IEEE Transactions on
Evolutionary Computation 13(4), 858–878 (Aug 2009)
6. Howley, T., Madden, M.G.: An evolutionary approach to automatic kernel construction. In: Proc. of ICANN 2006. LNCS, vol. 4132, pp. 417–426. Springer-Verlag
(2006)
7. Koch, P., Bischl, B., Flasch, O., Beielstein, T., Weihs, C., Konen, W.: Tuning and
evolution of support vector kernels. Evolutionary Intelligence 5, 153–170 (2012)
8. Koza, J.R.: Genetic programming: on the programming of computers by means of
natural selection. MIT Press, Cambridge, MA, USA (1992)
9. McKay, R.I., Hoai, N.X., Whigham, P.A., Shan, Y., O’Neill, M.: Grammar-based
genetic programming: a survey. Genetic Programming and Evolvable Machines
11(3/4), 365–396 (Sep 2010)
10. Rasmussen, C.E., Williams, C.K.: Gaussian Processes for Machine Learning. MIT
Press (2006)
11. Wagner, S.: Heuristic optimization software systems – Modeling of heuristic optimization algorithms in the HeuristicLab software environment. Ph.D. thesis, Institute for Formal Models and Verification, Johannes Kepler University, Linz (2009)
12. Wilson, A.G., Prescott Adams, R.: Gaussian Process Covariance Kernels for Pattern Discovery and Extrapolation. ArXiv e-prints (Feb 2013)

