AbstractGPs.jl

abstract type AbstractGP end

Supertype for various Gaussian process (GP) types. A common interface is provided for interacting with each of these objects. See [1] for an overview of GPs.

[1] - C. E. Rasmussen and C. Williams. "Gaussian processes for machine learning". MIT Press. 2006.

ConstMean{T<:Real} <: MeanFunction

Returns c everywhere.

CustomMean{Tf} <: MeanFunction

A wrapper around whatever unary function you fancy. Must be able to be mapped over an AbstractVector of inputs.

FiniteGP{Tf<:AbstractGP, Tx<:AbstractVector, TΣy}

The finite-dimensional projection of the AbstractGP f at x. Assumed to be observed under Gaussian noise with zero mean and covariance matrix Σ

GP{Tm<:MeanFunction, Tk<:Kernel}

A Gaussian Process (GP) with known mean and kernel. See e.g. [1] for an introduction.

Zero Mean

If only one argument is provided, assume the mean to be zero everywhere:

julia> f = GP(Matern32Kernel());

julia> x = randn(5);

julia> mean(f(x)) == zeros(5)
true

julia> cov(f(x)) == kernelmatrix(Matern32Kernel(), x)
true

Constant Mean

If a Real is provided as the first argument, assume the mean function is constant with that value

julia> f = GP(5.0, Matern32Kernel());

julia> x = randn(5);

julia> mean(f(x)) == 5.0 .* ones(5)
true

julia> cov(f(x)) == kernelmatrix(Matern32Kernel(), x)
true

Custom Mean

Provide an arbitrary function to compute the mean:

julia> f = GP(x -> sin(x) + cos(x / 2), Matern32Kernel());

julia> x = randn(5);

julia> mean(f(x)) == sin.(x) .+ cos.(x ./ 2)
true

julia> cov(f(x)) == kernelmatrix(Matern32Kernel(), x)
true

[1] - C. E. Rasmussen and C. Williams. "Gaussian processes for machine learning". MIT Press. 2006.

LatentFiniteGP(fx<:FiniteGP, lik)

• fx is a FiniteGP.
• lik is the log likelihood function which maps sample from f to corresposing

conditional likelihood distributions.

LatentGP(f<:GP, lik, Σy)

• f is a AbstractGP.
• lik is the log likelihood function which maps sample from f to corresposing

conditional likelihood distributions.

• Σy is the observation noise

ZeroMean{T<:Real} <: MeanFunction

Returns zero(T) everywhere.

approx_posterior(::VFE, fx::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

Compute the optimal approximate posterior [1] over the process f, given observations y of f at x, and inducing points u, where u = f(z) for some inducing inputs z.

[1] - M. K. Titsias. "Variational learning of inducing variables in sparse Gaussian processes". In: Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics. 2009.

cov_diag(f::AbstractGP, x::AbstractVector)

Compute only the diagonal elements of cov(f(x)).

dtc(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

The Deterministic Training Conditional (DTC) [1]. y are observations of f, and u are pseudo-points.

julia> f = GP(Matern52Kernel());

julia> x = randn(1000);

julia> z = range(-5.0, 5.0; length=256);

julia> y = rand(f(x, 0.1));

julia> isapprox(dtc(f(x, 0.1), y, f(z)), logpdf(f(x, 0.1), y); atol=1e-3, rtol=1e-3)
true

[1] - M. Seeger, C. K. I. Williams and N. D. Lawrence. "Fast Forward Selection to Speed Up Sparse Gaussian Process Regression". In: Proceedings of the Ninth International Workshop on Artificial Intelligence and Statistics. 2003

elbo(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

The Titsias Evidence Lower BOund (ELBO) [1]. y are observations of f, and u are pseudo-points, where u = f(z) for some z.

julia> f = GP(Matern52Kernel());

julia> x = randn(1000);

julia> z = range(-5.0, 5.0; length=13);

julia> y = rand(f(x, 0.1));

julia> elbo(f(x, 0.1), y, f(z)) < logpdf(f(x, 0.1), y)
true

[1] - M. K. Titsias. "Variational learning of inducing variables in sparse Gaussian processes". In: Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics. 2009.

marginals(f::FiniteGP)

Compute a vector of Normal distributions representing the marginals of f efficiently. In particular, the off-diagonal elements of cov(f(x)) are never computed.

julia> f = GP(Matern32Kernel());

julia> x = randn(11);

julia> fs = marginals(f(x));

julia> mean.(fs) == mean(f(x))
true

julia> std.(fs) == sqrt.(diag(cov(f(x))))
true

mean_and_cov(f::AbstractGP, x::AbstractVector)

Compute both mean(f(x)) and cov(f(x)). Sometimes more efficient than separately computation, particularly for posteriors.

mean_and_cov(f::FiniteGP)

Equivalent to (mean(f), cov(f)), but sometimes more efficient to compute them jointly than separately.

julia> fx = GP(SqExponentialKernel())(range(-3.0, 3.0; length=10), 0.1);

julia> mean_and_cov(fx) == (mean(fx), cov(fx))
true

mean_and_cov_diag(f::AbstractGP, x::AbstractVector)

Compute both mean(f(x)) and the diagonal elements of cov(f(x)). Sometimes more efficient than separately computation, particularly for posteriors.

posterior(fx::FiniteGP, y::AbstractVector{<:Real})

Constructs the posterior distribution over fx.f given observations y at x made under noise fx.Σy. This is another AbstractGP object. See chapter 2 of [1] for a recap on exact inference in GPs. This posterior process has mean function

m_posterior(x) = m(x) + k(x, fx.x) inv(cov(fx)) (y - mean(fx))

and kernel

k_posterior(x, z) = k(x, z) - k(x, fx.x) inv(cov(fx)) k(fx.x, z)

where m and k are the mean function and kernel of fx.f respectively.

posterior(fx::FiniteGP{<:PosteriorGP}, y::AbstractVector{<:Real})

Constructs the posterior distribution over fx.f when f is itself a PosteriorGP by updating the cholesky factorisation of the covariance matrix and avoiding recomputing it from original covariance matrix. It does this by using update_chol functionality.

Other aspects are similar to a regular posterior.

sampleplot(GP::FiniteGP, samples)

Plot samples from the given FiniteGP. Make sure to run using Plots before using this function.

Example

using Plots
f = GP(SqExponentialKernel())
sampleplot(f(rand(10), 10; markersize=5)

The given example plots 10 samples from the given FiniteGP. The markersize is modified from default of 0.5 to 5.

sampleplot(GP::FiniteGP, samples)

Plot samples from the given FiniteGP. Make sure to run using Plots before using this function.

Example

using Plots
f = GP(SqExponentialKernel())
sampleplot(f(rand(10), 10; markersize=5)

The given example plots 10 samples from the given FiniteGP. The markersize is modified from default of 0.5 to 5.

sampleplot(GP::FiniteGP, samples)

Plot samples from the given FiniteGP. Make sure to run using Plots before using this function.

Example

using Plots
f = GP(SqExponentialKernel())
sampleplot(f(rand(10), 10; markersize=5)

The given example plots 10 samples from the given FiniteGP. The markersize is modified from default of 0.5 to 5.

function update_approx_posterior(
    f_post_approx::ApproxPosteriorGP,
    fx::FiniteGP,
    y::AbstractVector{<:Real}
)

Update the ApproxPosteriorGP given a new set of observations. Here, we retain the same of pseudo-points.

function update_approx_posterior(
    f_post_approx::ApproxPosteriorGP,
    u::FiniteGP,
)

Update the ApproxPosteriorGP given a new set of pseudo-points to append to the existing set of pseudo points.

 update_chol(chol::Cholesky, C12::AbstractMatrix, C22::AbstractMatrix)

Let C be the positive definite matrix comprising blocks

C = [C11 C12;
     C21 C22]

with upper-triangular cholesky factorisation comprising blocks

U = [U11 U12;
     0   U22]

where U11 and U22 are themselves upper-triangular, and U11 = cholesky(C11).U. update_chol computes the updated Cholesky given original chol, C12, and C22.

Arguments

 - chol::Cholesky: The original cholesky decomposition
 - C12::AbstractMatrix: matrix of size (size(chol.U, 1), size(C22, 1))
 - C22::AbstractMatrix: positive-definite matrix

rand(rng::AbstractRNG, f::FiniteGP, N::Int=1)

Obtain N independent samples from the marginals f using rng. Single-sample methods produce a length(f) vector. Multi-sample methods produce a length(f) x NMatrix.

julia> f = GP(Matern32Kernel());

julia> x = randn(11);

julia> rand(f(x)) isa Vector{Float64}
true

julia> rand(MersenneTwister(123456), f(x)) isa Vector{Float64}
true

julia> rand(f(x), 3) isa Matrix{Float64}
true

julia> rand(MersenneTwister(123456), f(x), 3) isa Matrix{Float64}
true

logpdf(f::FiniteGP, y::AbstractVecOrMat{<:Real})

The logpdf of y under f if is y isa AbstractVector. logpdf of each column of y if y isa Matrix.

julia> f = GP(Matern32Kernel());

julia> x = randn(11);

julia> y = rand(f(x));

julia> logpdf(f(x), y) isa Real
true

julia> Y = rand(f(x), 3);

julia> logpdf(f(x), Y) isa AbstractVector{<:Real}
true

logpdf(lfgp::LatentFiniteGP, y::NamedTuple{(:f, :y)})

Returns the joint log density of the gaussian process output f and real output y.

cov(f::AbstractGP, x::AbstractVector, y::AbstractVector)

Compute the length(x) by length(y) cross-covariance matrix between f(x) and f(y).

cov(f::AbstractGP, x::AbstractVector)

Compute the length(x) by length(x) covariance matrix of the multivariate Normal f(x).

cov(fx::FiniteGP, gx::FiniteGP)

Compute the cross-covariance matrix between fx and gx.

julia> f = GP(Matern32Kernel());

julia> x1 = randn(11);

julia> x2 = randn(13);

julia> cov(f(x1), f(x2)) == kernelmatrix(Matern32Kernel(), x1, x2)
true

cov(f::FiniteGP)

Compute the covariance matrix of fx.

Noise-free observations

julia> f = GP(Matern52Kernel());

julia> x = randn(11);

julia> cov(f(x)) == kernelmatrix(Matern52Kernel(), x)
true

Isotropic observation noise

julia> cov(f(x, 0.1)) == kernelmatrix(Matern52Kernel(), x) + 0.1 * I
true

Independent anisotropic observation noise

julia> s = rand(11);

julia> cov(f(x, s)) == kernelmatrix(Matern52Kernel(), x) + Diagonal(s)
true

Correlated observation noise

julia> A = randn(11, 11); S = A'A;

julia> cov(f(x, S)) == kernelmatrix(Matern52Kernel(), x) + S
true

mean(f::AbstractGP, x::AbstractVector)

Computes the mean vector of the multivariate Normal f(x).

mean(fx::FiniteGP)

Compute the mean vector of fx.

julia> f = GP(Matern52Kernel());

julia> x = randn(11);

julia> mean(f(x)) == zeros(11)
true