Custom covariance

This example shows how to define a covariance function outside GPmp’s built-in Matern helpers and use it in a one-dimensional interpolation example.

What this example does

The script defines a user covariance callable with the same signature expected by gpmp.core.Model: covariance(x, y, covparam, pairwise=False). It then builds a GP model from that covariance, selects covariance parameters, and predicts on a dense one-dimensional grid.

Outputs

The displayed quantities are the reference function, observations, posterior mean, and uncertainty envelope obtained with the custom covariance. For a positive definite covariance and converged parameter selection, the posterior mean interpolates the noise-free observations and the uncertainty envelope widens where observations are sparse.

API points

  • Custom covariance functions should use gpmp.num operations when they must work with both NumPy and torch backends.

  • The covariance callable must accept the pairwise argument used internally by GPmp prediction and likelihood routines.

  • Once the model is built, parameter selection and prediction use the same API as for built-in covariance functions.

Covariance construction

gpmp.core.Model calls the covariance function in three situations:

k(x_i, x_i)

Covariance matrix between observation points. In the example, this is handled by kernel_ii_or_tt.

k(x_i, x_t)

Cross-covariance matrix between observation points and prediction points. In the example, this is handled by kernel_it.

k(x_t, x_t)

Prior covariance at prediction points. GPmp uses only its diagonal when computing posterior variances, unless a full posterior covariance matrix is requested. In the example, this is also handled by kernel_ii_or_tt.

The suffix ii_or_tt means “same-set covariance”: the two arguments are the same point set, either observations x_i or prediction points x_t. The suffix it means “cross covariance” between observations and prediction points.

The pairwise flag controls the returned shape. With pairwise=False, the covariance function returns a full matrix. With pairwise=True, it returns the elementwise covariance vector [k(x_0, y_0), ..., k(x_{n-1}, y_{n-1})]. When y is None, this is the diagonal of k(x, x).

The wrapper kernel dispatches to the same-set or cross-set function. This is the callable passed to gpmp.core.Model. The small nugget added in kernel_ii_or_tt is numerical jitter on same-set covariance matrices. It is not an observation-noise model.

The example uses a one-dimensional Matern covariance with

\[k_\theta(x, x') = \sigma^2 c_2\left(|x-x'|/\rho\right), \qquad \theta = \left(\log(\sigma^2), -\log(\rho)\right).\]

The custom part of the example is not the mathematical covariance itself, but the way the covariance callable is written and dispatched.

../_images/custom_kernel_0_0.png

Script: examples/gpmp_example05_1d_custom_kernel.py

  1"""
  2GP interpolation in 1D, with noiseless data
  3
  4This example demonstrates how to compute GP interpolation with unknown mean
  5(aka ordinary / intrinsic kriging) on a one-dimensional noiseless dataset.
  6
  7A Mat'ern covariance function is used for the Gaussian Process (GP) prior.
  8The parameters of this covariance function are assumed to be known
  9(i.e., no parameter estimation is performed here).
 10
 11The kriging predictor / posterior mean of the GP, interpolates the data.
 12
 13----
 14Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
 15Copyright (c) 2022-2026, CentraleSupelec
 16License: GPLv3 (see LICENSE)
 17----
 18This example is based on the file stk_example_kb01.m from the STK at
 19https://github.com/stk-kriging/stk/
 20by Julien Bect and Emmanuel Vazquez, released under the GPLv3 license.
 21
 22Original copyright notice:
 23
 24   Copyright (c) 2015, 2016, 2018 CentraleSupelec
 25   Copyright (c) 2011-2014 SUPELEC
 26"""
 27
 28import math
 29import numpy as np
 30import gpmp.num as gnp
 31import gpmp as gp
 32
 33
 34def generate_data():
 35    """
 36    Data generation
 37    (xt, zt): target
 38    (xi, zi): input dataset
 39    """
 40    # Build (xt, zt)
 41    dim = 1
 42    nt = 200
 43    box = [[-1], [1]]
 44    xt = gp.misc.designs.regulargrid(dim, nt, box)
 45    zt = gp.misc.testfunctions.twobumps(xt)
 46
 47    shuffle = True
 48    if shuffle:
 49        ni = 5
 50        ind = np.random.choice(nt, ni, replace=False)
 51    else:
 52        ind = [10, 45, 100, 130, 160]
 53    xi = xt[ind]
 54    zi = zt[ind]
 55
 56    return xt, zt, xi, zi
 57
 58
 59def zero_mean(x, param):
 60    return None
 61
 62
 63def constant_mean(x, param):
 64    return gnp.ones((x.shape[0], 1))
 65
 66
 67def kernel_ii_or_tt(x, param, pairwise=False):
 68    """Same-set covariance k(x, x).
 69
 70    This function is used for observation-observation covariance k(x_i, x_i)
 71    and prediction-prediction covariance k(x_t, x_t).  With pairwise=True,
 72    GPmp asks only for the diagonal vector.
 73    """
 74    p = 2
 75    sigma2 = gnp.exp(param[0])
 76    loginvrho = param[1]
 77    # Numerical jitter for same-set covariance matrices, not observation noise.
 78    nugget = 100 * gnp.eps
 79
 80    if pairwise:
 81        # Return diag(k(x, x)) as a vector of shape (n,).
 82        K = sigma2 * gnp.ones((x.shape[0], ))
 83    else:
 84        # Return the full same-set covariance matrix of shape (n, n).
 85        K = gnp.scaled_distance(loginvrho, x, x)
 86        K = sigma2 * gp.kernel.maternp_kernel(p, K) + nugget * gnp.eye(K.shape[0])
 87
 88    return K
 89
 90
 91def kernel_it(x, y, param, pairwise=False):
 92    """Cross-covariance k(x, y) between two point sets."""
 93    p = 2
 94    sigma2 = gnp.exp(param[0])
 95    loginvrho = param[1]
 96
 97    if pairwise:
 98        # Return [k(x_i, y_i)] as a vector of shape (n,).
 99        K = gnp.scaled_distance_elementwise(loginvrho, x, y)
100    else:
101        # Return the full cross-covariance matrix of shape (n_x, n_y).
102        K = gnp.scaled_distance(loginvrho, x, y)
103
104    K = sigma2 * gp.kernel.maternp_kernel(p, K)
105    return K
106
107
108def kernel(x, y, param, pairwise=False):
109    """Dispatch covariance calls made by gpmp.core.Model."""
110    if y is x or y is None:
111        return kernel_ii_or_tt(x, param, pairwise)
112    else:
113        return kernel_it(x, y, param, pairwise)
114
115
116def visualize(xt, zt, xi, zi, zpm, zpv):
117    fig = gp.plot.Figure(isinteractive=True)
118    fig.plot(xt, zt, 'C0', linestyle=(0, (5, 5)), linewidth=1.0)
119    fig.plotdata(xi, zi)
120    fig.plotgp(xt, zpm, zpv)
121    fig.xylabels('x', 'z')
122    fig.show(grid=True, legend=True, legend_fontsize=9)
123
124
125def main():
126    xt, zt, xi, zi = generate_data()
127    mean = constant_mean
128    meanparam = None
129
130    covparam = gnp.array([math.log(0.5**2),    # log(sigma2)
131                          math.log(1 / .7)])   # log(1/rho)
132
133    model = gp.core.Model(mean, kernel, meanparam, covparam)
134
135    zpm, zpv = model.predict(xi, zi, xt)
136    visualize(xt, zt, xi, zi, zpm, zpv)
137
138
139if __name__ == "__main__":
140    main()