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.

Warning

CustomMean is generally sufficient for testing purposes, but care should be taken if attempting to differentiate through mean_vector with a CustomMean when using Zygote.jl. In particular, mean_vector(m::CustomMean, x) is implemented as map(m.f, x), which when x is a ColVecs or RowVecs will not differentiate correctly.

In such cases, you should implement mean_vector directly for your custom mean. For example, if f(x) = sum(x), you might implement mean_vector as

mean_vector(::CustomMean{typeof(f)}, x::ColVecs) = vec(sum(x.X; dims=1))
mean_vector(::CustomMean{typeof(f)}, x::RowVecs) = vec(sum(x.X; dims=2))

which avoids ever applying map to a ColVecs or RowVecs.

DTC(fz::FiniteGP)

Similar to VFE, but uses a different objective for approx_log_evidence.

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 Σy

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 likelihood function which maps samples from f to the corresponding conditional likelihood distributions (i.e., lik must return a Distribution compatible with the observations).

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

• f is a AbstractGP.
• lik is the likelihood function which maps samples from f to the corresponding conditional likelihood distributions (i.e., lik must return a Distribution compatible with the observations).
• Σy is the noise under which the latent GP is "observed"; this represents the jitter used to avoid numeric instability and should generally be small.

abstract type MeanFunction end

MeanFunction introduces an API for treating the prior mean function appropriately. On the abstract level, all MeanFunction are functions. However we generally want to evaluate them on a collection of inputs. To this effect, we provide the mean_vector(::MeanFunction, ::AbstractVector) function, which is equivalent to map but with possibilities of optimizations (for ZeroMean and ConstMean for example).

VFE(fz::FiniteGP)

The "Variational Free Energy" sparse approximation [1], used to construct an approximate posterior with inducing inputs fz.x. See posterior(v::VFE, fx::FiniteGP, y::AbstractVector{<:Real}) for a usage example.

[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.

ZeroMean{T<:Real} <: MeanFunction

Returns zero(T) everywhere, T is Float64 by default.

approx_log_evidence(approx::<Approximation>, lfx::LatentFiniteGP, ys)

Compute an approximation to the log of the marginal likelihood (also known as "evidence") under the given approximation to the posterior. The return value of approx_log_evidence can be used to optimise the hyperparameters of lfx.

approx_log_evidence(dtc::DTC, fx::FiniteGP, y::AbstractVector{<:Real})

The Deterministic Training Conditional (DTC) [1]. y are observations of fx, and v.z are inducing points.

julia> f = GP(Matern52Kernel());

julia> x = randn(1000);

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

julia> d = DTC(f(z));

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

julia> isapprox(approx_log_evidence(d, f(x, 0.1), y), logpdf(f(x, 0.1), y); atol=1e-6, rtol=1e-6)
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

approx_log_evidence(vfe::VFE, fx::FiniteGP, y::AbstractVector{<:Real})
elbo(vfe::VFE, fx::FiniteGP, y::AbstractVector{<:Real})

The Titsias Evidence Lower BOund (ELBO) [1]. y are observations of fx, and v.z are inducing points.

julia> f = GP(Matern52Kernel());

julia> x = randn(1000);

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

julia> v = VFE(f(z));

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

julia> elbo(v, f(x, 0.1), y) < 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_vector(m::MeanFunction, x::AbstractVector)::AbstractVector{<:Real}

mean_vector is the function to call to apply a MeanFunction to a collection of inputs.

posterior(fx::FiniteGP, y::AbstractVector{<:Real})
posterior(approx::<Approximation>, fx::FiniteGP, y::AbstractVector{<:Real})
posterior(approx::<Approximation>, lfx::LatentFiniteGP, y::AbstractVector)

Construct the posterior distribution over the latent Gaussian process (fx.f or lfx.fx.f), given the observations y corresponding to the process's finite projection (fx or lfx).

In the two-argument form, this describes exact GP regression with y observed under a Gaussian likelihood, and returns a PosteriorGP.

In the three-argument form, the first argument specifies the approximation to be used (e.g. VFE or defined in other packages such as ApproximateGPs.jl), and returns an ApproxPosteriorGP.

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

Construct the posterior distribution over fx.f given observations y at fx.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})

Construct 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 the original covariance matrix. It does this by using update_chol functionality.

Other aspects are similar to a regular posterior.

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

Compute the optimal approximate posterior [1] over the process f = fx.f, given observations y of f at x, and inducing points vfe.fz.x.

julia> f = GP(Matern52Kernel());

julia> x = randn(1000);

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

julia> vfe = VFE(f(z));

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

julia> post = posterior(vfe, f(x, 0.1), y);

julia> post(z) isa AbstractGPs.FiniteGP
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.

sampleplot([x::AbstractVector=f.x, ]f::FiniteGP; samples=1, kwargs...)

Plot samples from the projection f of a Gaussian process versus x.

When plotting multiple samples, these are treated as a single series (i.e., only a single entry will be added to the legend when providing a label).

Example

using Plots

gp = GP(SqExponentialKernel())
sampleplot(gp(rand(5)); samples=10, linealpha=1.0)

The given example plots 10 samples from the projection of the GP gp. The linealpha is modified from default of 0.35 to 1.

sampleplot(x::AbstractVector, gp::AbstractGP; samples=1, kwargs...)

Plot samples from the finite projection gp(x, 1e-9) versus x.

sampleplot([x::AbstractVector=f.x, ]f::FiniteGP; samples=1, kwargs...)

Plot samples from the projection f of a Gaussian process versus x.

When plotting multiple samples, these are treated as a single series (i.e., only a single entry will be added to the legend when providing a label).

Example

using Plots

gp = GP(SqExponentialKernel())
sampleplot(gp(rand(5)); samples=10, linealpha=1.0)

The given example plots 10 samples from the projection of the GP gp. The linealpha is modified from default of 0.35 to 1.

sampleplot(x::AbstractVector, gp::AbstractGP; samples=1, kwargs...)

Plot samples from the finite projection gp(x, 1e-9) versus x.

sampleplot([x::AbstractVector=f.x, ]f::FiniteGP; samples=1, kwargs...)

Plot samples from the projection f of a Gaussian process versus x.

When plotting multiple samples, these are treated as a single series (i.e., only a single entry will be added to the legend when providing a label).

Example

using Plots

gp = GP(SqExponentialKernel())
sampleplot(gp(rand(5)); samples=10, linealpha=1.0)

The given example plots 10 samples from the projection of the GP gp. The linealpha is modified from default of 0.35 to 1.

sampleplot(x::AbstractVector, gp::AbstractGP; samples=1, kwargs...)

Plot samples from the finite projection gp(x, 1e-9) versus x.

 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

function update_posterior(
    f_post_approx::ApproxPosteriorGP{<:Union{VFE,DTC}},
    fx::FiniteGP,
    y::AbstractVector{<:Real}
)

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

function update_posterior(
    f_post_approx::ApproxPosteriorGP{<:Union{VFE,DTC}},
    z::FiniteGP,
)

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

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) × N Matrix.

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 y isa AbstractVector. The 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)})

\[ log p(y, f; x)\]

The joint log density of the Gaussian process output f and observation y.

rand!(rng::AbstractRNG, f::FiniteGP, y::AbstractVecOrMat{<:Real})

Obtain sample(s) from the marginals f using rng and write them to y.

If y is a matrix, then each column corresponds to an independent sample.

julia> f = GP(Matern32Kernel());

julia> x = randn(11);

julia> y = similar(x);

julia> rand!(f(x), y);

julia> rand!(MersenneTwister(123456), f(x), y);

julia> ys = similar(x, length(x), 3);

julia> rand!(f(x), ys);

julia> rand!(MersenneTwister(123456), f(x), ys);

plot(f::FiniteGP; kwargs...)
plot!([plot, ]f::FiniteGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection f of a Gaussian process versus f.x.

plot(x::AbstractVector, gp::AbstractGP; kwargs...)
plot!([plot, ]x::AbstractVector, gp::AbstractGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection gp(x) of the Gaussian process gp.

plot(x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)
plot!([plot, ]x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)

Plot the predictive mean for the projection f of a Gaussian process and a ribbon of ribbon_scale standard deviations around it versus x.

Examples

Plot the mean and a ribbon of 3 standard deviations:

using Plots

gp = GP(SqExponentialKernel())
plot(gp(rand(5)); ribbon_scale=3)

plot(f::FiniteGP; kwargs...)
plot!([plot, ]f::FiniteGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection f of a Gaussian process versus f.x.

plot(x::AbstractVector, gp::AbstractGP; kwargs...)
plot!([plot, ]x::AbstractVector, gp::AbstractGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection gp(x) of the Gaussian process gp.

plot(x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)
plot!([plot, ]x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)

Plot the predictive mean for the projection f of a Gaussian process and a ribbon of ribbon_scale standard deviations around it versus x.

Examples

Plot the mean and a ribbon of 3 standard deviations:

using Plots

gp = GP(SqExponentialKernel())
plot(gp(rand(5)); ribbon_scale=3)

plot(f::FiniteGP; kwargs...)
plot!([plot, ]f::FiniteGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection f of a Gaussian process versus f.x.

plot(x::AbstractVector, gp::AbstractGP; kwargs...)
plot!([plot, ]x::AbstractVector, gp::AbstractGP; kwargs...)

Plot the predictive mean and a ribbon around it for the projection gp(x) of the Gaussian process gp.

plot(x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)
plot!([plot, ]x::AbstractVector, f::FiniteGP; ribbon_scale=1, kwargs...)

Plot the predictive mean for the projection f of a Gaussian process and a ribbon of ribbon_scale standard deviations around it versus x.

Examples

Plot the mean and a ribbon of 3 standard deviations:

using Plots

gp = GP(SqExponentialKernel())
plot(gp(rand(5)); ribbon_scale=3)

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)

Compute 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

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

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

var(f::FiniteGP)

Compute only the diagonal elements of cov(f).

Examples

julia> fx = GP(Matern52Kernel())(randn(10), 0.1);

julia> var(fx) == diag(cov(fx))
true

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

Compute both mean(f(x)) and cov(f(x)). Sometimes more efficient than computing them separately, 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_var(f::AbstractGP, x::AbstractVector)

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

mean_and_var(f::FiniteGP)

Compute both mean(f) and the diagonal elements of cov(f).

Sometimes more efficient than computing them separately, particularly for posteriors.

Examples

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

julia> mean_and_var(fx) == (mean(fx), var(fx))
true

test_finitegp_primary_and_secondary_public_interface(
    rng::AbstractRNG, fx::FiniteGP; atol::Real=1e-12
)

Basic consistency tests for both the Primary and Secondary Public FiniteGP APIs. Runs test_finitegp_primary_public_interface as part of these tests. You should run these tests if you implement both the primary and secondary public APIs – see API section of the docs for more information about these APIs.

These are consistency checks, not correctness tests in the absolute sense. For example, these tests ensure that samples generated by rand are of the correct size, but does not check that they come from the intended distribution.

test_finitegp_primary_public_interface(
    rng::AbstractRNG, fx::FiniteGP; atol::Real=1e-12
)

Basic consistency tests for the Primary Public FiniteGP API. You should run these tests if you only implement the Primary Public API – see API section of docs for details about this API.

These are consistency checks, not correctness tests in the absolute sense. For example, these tests ensure that samples generated by rand are of the correct size, but does not check that they come from the intended distribution.

test_internal_abstractgps_interface(
    rng::AbstractRNG,
    f::AbstractGP,
    x::AbstractVector,
    z::AbstractVector;
    atol=1e-12,
    σ²::Real=1e-1,
    jitter::Real=AbstractGPs.default_σ²,
)

Basic consistency tests for the Internal AbstractGPs API. Runs test_finitegp_primary_and_secondary_public_interface as part of these tests. Run these tests if you implement the Internal AbstractGPs API – see the API section of the docs for more information about this API.

These are consistency checks, not correctness tests in the absolute sense. For example, these tests ensure that samples generated by rand are of the correct size, but does not check that they come from the intended distribution.

kwargs

• atol: the minimum eigenvalue of cov(f, x) must be greater than this.
• σ²: the observation noise under which things are observed.
• jitter: constant added to diagonal of cov(f, x) and cov(f, z) before inversion.