GoGP is a library for probabilistic programming around Gaussian processes. It uses Infergo for automatic differentiation and inference.

GoGP is built around a dual view on the Gaussian process

• as a stochastic process,
• as a probabilistic model with respect to kernel.

Gaussian process instance

GP, the Gaussian process type, encapsulates similarity and noise kernels, their parameters, and observation inputs and outputs:

// Type GP is the barebone implementation of GP.
type GP struct {
	NDim                   int       // number of dimensions
	Simil, Noise           Kernel    // kernels
	ThetaSimil, ThetaNoise []float64 // kernel parameters

	X [][]float64 // inputs
	Y []float64   // outputs

	Parallel bool // when true, covariances are computed in parallel
}

Public methods defined on GP fall into two groups: Gaussian process fitting and prediction, on one hand, and probabilistic model interface, on the other hand.

Gaussian process methods

Absorb updates GP with observations, Produce returns predicted outputs for unseen inputs.

func (*GP) Absorb(x [][]float64, y []float64)
func (*GP) Produce(x [][]float64) (mu, sigma []float64)

When kernel parameters are fixed (known or learned), Absorb and Produce are all that is needed for posterior Gaussian process inference.

Probabilistic model methods

Observe and Gradient turn a GP instance into an elemental Infergo model (that is, a model with supplied gradient). The model can be passed to inference algorithms or used within another model.

func (*GP) Observe(x []float64) float64
func (*GP) Gradient() []float64 // `elemental' model

If the length of x is equal to the number of kernel hyperparameters (Simil.NTheta() + Noise.NTheta()) then the gradient of marginal likelihood is computed with respect to kernel hyperparameters only. Observations must be provided in the fields of the GP instance. Otherwise, if the length of x is greater than the number of parameters, the rest of x is interpreted as observations. In the latter case, the gradient is computed with respect to both kernel hyperparameters and observations.

Hyperparameter priors

Type Model is a wrapper model combining a GP instance and priors on hyperparameters into an elemental Infergo model.

type Model struct {
	*GP
	Priors model.Model
}

Kernels

There are two kernel kinds:

• similarity kernel;
• noise kernel.

Both kinds must satisfy the Kernel interface:

type Kernel interface {
	Observe([]float64) float64
	NTheta() int
}

The Observe method computes the variance or covariance, the NTheta method returns the number of kernel parameters.

A similarity (or covariance) kernel receives concanetation of kernel parameters and coordinates of two points. A noise kernel receives concatenation of kernel parameters and coordinates of a single point. Here is an example implementation of the RBF (or normal) kernel:

type RBF struct{}

func (RBF) Observe(x []float64) float64 {
	l, xa, xb := x[0], x[1], x[2]
	d := (xa - xb) / l
	return math.Exp(-d * d / 2)
}

func (RBF) NTheta() int { return 1 }

Frequently used primitive kernels are included in the library.

Case studies

GoGP includes case studies, illustrating, on simple examples, common patterns of GoGP use. We briefly summarize here some of the case studies.

Basic usage

In the basic case, similar to that supported by many Gaussian process libraries, a GP instance directly serves as the model for inference on hyperparameters (or the hyperparameters can be just fixed).

The library user specifies the kernel:

type Basic struct{}
func (Basic) Observe(x []float64) float64 {
    return x[0] * kernel.Normal.Observe(x[1:])
}
func (Basic) NTheta() int { return 2 }

gp := &gp.GP{
    NDim:  1,
    Simil: Basic{},
}

MLE inference on hyperparameters and prediction can then be performed through library functions.

Priors on hyperparameters

If priors on hyperparameters are to be specified, the library user provides both the kernel and the prior beliefs:

// Similarity kernel
type Simil struct{}
func (Simil) Observe(x []float64) float64 {
	const (
		c = iota // trend scale
		l        // trend length scale
		xa        // first point
		xb        // second point
	)

	return x[c]*kernel.Matern52.Cov(x[l], x[xa], x[xb])
}
func (Simil) NTheta() int { return 2 }

// Noise kernel
type Noise struct{}
func (Noise) Observe(x []float64) float64 {
	return 0.01 * kernel.UniformNoise.Observe(x)
}
func (Noise) NTheta() int { return 1 }

// Hyperparameter priors
type Priors struct{}
func (*Priors) Observe(x []float64) float64 {
	return Normal.Logps(1, 1, x...)
}

m := &gp.Model{
	GP: &gp.GP{
		NDim:     1,
		Simil:    Simil,
		Noise:    Noise,
	},
	Priors: &Priors{},
}

Uncertain observation inputs

When observation inputs are uncertain, beliefs about inputs can be specified, and the log-likelihood gradient can be computed with respect to both hyperparameters and observation inputs. For example, in a time series, one can assume that observation inputs come from a renewal process and let the inputs move relative to each other. Then, forecasting can be performed relative to posterior observation inputs.

Non-Gaussian noise

In basic usage, the observation noise is assumed to be Gaussian. This usage is supported by initializing GP with a noise kernel, along with a similarity kernel. When the noise is not Gaussian, an analytical solution for posterior Gaussian process inference does not always exist. However, non-Gaussian noise is straightforwardly supported by GoGP through supplying a model for beliefs on observation outputs:

type Noise struct {
	Y []float64 // noisy outputs
}
func (m *Noise) Observe(x []float64) (lp float64){
	// Laplacian noise
	for i := range m.Y {
		lp += Expon.Logp(1/math.Exp(x[0]),
			math.Abs(m.Y[i]-x[1+i]))
	}
	return lp
}