Higher-dimensional interpolation

This example applies anisotropic Matérn GP interpolation to a test function in dimension greater than two. Since direct contour plots are no longer available, the preview uses prediction scatter plots, leave-one-out (LOO) diagnostics, and one-dimensional cross sections.

What this example does

The script chooses a benchmark function, evaluates it at observation points, selects covariance parameters with REMAP, predicts independent test points, and computes LOO predictions. LOO prediction removes one observation at a time, predicts it from all other observations, and compares the prediction with the removed value.

Mathematical description

The model uses the universal-kriging formulation from the one-dimensional example:

\[Z(x)=p(x)^\top\beta+Z_0(x), \qquad Z_0 \sim \mathcal{GP}(0, k_\theta), \qquad Z_i = Z(x_i),\]

and zi stores realized values \(z_i\). In this script the mean basis is constant. The input dimension is too large for a direct contour plot. The first diagnostic therefore compares realized reference values \(z_t\) with the posterior mean \(\mathbb{E}[Z_t\mid Z_i=z_i]\) on independent test points.

For the LOO diagnostic, each observation is predicted after removing it from the conditioning set:

\[\widehat z_{i,-i} = \mathbb{E}\left[ Z_i \mid Z_j=z_j,\quad j\ne i \right].\]

The cross-section plots fix all coordinates except one and display \(x_j \mapsto \mathbb{E}[Z(x)\mid Z_i=z_i]\) and its pointwise uncertainty.

Outputs

The first figure compares GP posterior means with reference values on independent test points. Points close to the diagonal indicate small test errors.

The second figure compares LOO predictions with observed values. Vertical bars show nominal predictive intervals. Large deviations from the diagonal or many observations outside the intervals indicate poorly selected covariance parameters, missing structure, or predictive intervals that are too narrow.

The last figure shows prediction cross sections. Each panel fixes all coordinates except one and plots the posterior mean and intervals along that coordinate. Black points are projected observations. The red point is the observation used as the cross-section anchor.

Functions used

  • model.loo(xi, zi) returns LOO means, variances, and errors.

  • gp.plot.plot_loo plots LOO diagnostics for high-dimensional problems where spatial plots are impossible.

  • gp.plot.crosssections plots one-dimensional slices through selected observation points.

  • REMAP parameter selection can be used when pure likelihood criteria produce poorly identified lengthscales.

../_images/interpolation_nd_0_0.png ../_images/interpolation_nd_0_1.png ../_images/interpolation_nd_0_2.png

Script: examples/gpmp_example04_nd.py

  1"""Prediction of some classical test functions in dimension > 2
  2
  3An anisotropic Matern covariance function is used for the Gaussian
  4Process (GP) prior. The parameters of this covariance function
  5(variance and ranges) are estimated using the Restricted Maximum
  6A Posteriori (ReMAP).
  7
  8The mean function of the GP prior is assumed to be constant and
  9unknown.
 10
 11The function is sampled on a space-filling Latin Hypercube design, and
 12the data is assumed to be noiseless.
 13
 14----
 15Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
 16Copyright (c) 2022-2026, CentraleSupelec
 17License: GPLv3 (see LICENSE)
 18"""
 19import gpmp.num as gnp
 20import gpmp as gp
 21
 22
 23def choose_test_case(problem):
 24    if problem == 1:
 25        problem_name = "Hartmann4"
 26        f = gp.misc.testfunctions.hartmann4
 27        dim = 4
 28        box = [[0.0] * 4, [1.0] * 4]
 29        ni = 40
 30        xi = gp.misc.designs.ldrandunif(dim, ni, box)
 31        nt = 1000
 32        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 33
 34    elif problem == 2:
 35        problem_name = "Hartmann6"
 36        f = gp.misc.testfunctions.hartmann6
 37        dim = 6
 38        box = [[0.0] * 6, [1.0] * 6]
 39        ni = 200
 40        xi = gp.misc.designs.ldrandunif(dim, ni, box)
 41        nt = 1000
 42        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 43
 44    elif problem == 3:
 45        problem_name = "Borehole"
 46        f = gp.misc.testfunctions.borehole
 47        dim = 8
 48        box = [
 49            [0.05, 100.0, 63070.0, 990.0, 63.1, 700.0, 1120.0, 9855.0],
 50            [0.15, 50000.0, 115600.0, 1110.0, 116.0, 820.0, 1680.0, 12045.0],
 51        ]
 52        ni = 30
 53        xi = gp.misc.designs.maximinldlhs(dim, ni, box)
 54        nt = 1000
 55        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 56
 57    elif problem == 4:
 58        problem_name = "detpep8d"
 59        f = gp.misc.testfunctions.detpep8d
 60        dim = 8
 61        box = [[0.0] * 8, [1.0] * 8]
 62        ni = 60
 63        xi = gp.misc.designs.maximinldlhs(dim, ni, box)
 64        nt = 1000
 65        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 66
 67    elif problem == 5:
 68        problem_name = "Ishigami"
 69        f = gp.misc.testfunctions.ishigami
 70        dim = 3
 71        box = [[-gnp.pi] * 3, [gnp.pi] * 3]
 72        ni = 80
 73        xi = gp.misc.designs.ldrandunif(dim, ni, box)
 74        nt = 1000
 75        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 76
 77    return problem_name, f, dim, box, ni, xi, nt, xt
 78
 79
 80def constant_mean(x, param):
 81    return gnp.ones((x.shape[0], 1))
 82
 83
 84def kernel(x, y, covparam, pairwise=False):
 85    p = 10
 86    return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 87
 88
 89def visualize_predictions(problem_name, zt, zpm):
 90    fig = gp.plot.Figure()
 91    fig.plot(zt, zpm, "ko", markersize=3)
 92    (xmin, xmax), (ymin, ymax) = fig.ax.get_xlim(), fig.ax.get_ylim()
 93    xmin = min(xmin, ymin)
 94    xmax = max(xmax, ymax)
 95    fig.plot([xmin, xmax], [xmin, xmax], "--", linewidth=1)
 96    fig.xylabels("reference values", "posterior mean")
 97    fig.title(f"{problem_name}: test predictions")
 98    fig.grid()
 99    fig.show()
100
101
102def main():
103    problem = 5
104    problem_name, f, dim, box, ni, xi, nt, xt = choose_test_case(problem)
105
106    zi = f(xi)
107    zt = f(xt)
108
109    model = gp.Model(constant_mean, kernel)
110    model, info = gp.kernel.select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior(
111        model, xi, zi, info=True
112    )
113    gp.modeldiagnosis.diag(
114        model, "linear_mean_maternp_anisotropic", info, xi, zi
115    )
116
117    (zpm, zpv) = model.predict(xi, zi, xt)
118
119    visualize_predictions(problem_name, zt, zpm)
120
121    zloom, zloov, eloo = model.loo(xi, zi)
122    gp.plot.plot_loo(zi, zloom, zloov)
123
124    gp.plot.crosssections(
125        model, xi, zi, box, ind_i=[0, 1], ind_dim=list(range(dim))
126    )
127
128    gp.modeldiagnosis.perf(
129        model,
130        xi,
131        zi,
132        loo=True,
133        loo_res=(zloom, zloov, eloo),
134        xtzt=(xt, zt),
135        zpmzpv=(zpm, zpv),
136    )
137
138
139if __name__ == "__main__":
140    main()