ML / REML / REMAP parameter selection

This example compares the main covariance-parameter selection criteria used in GPmp: maximum likelihood (ML), restricted maximum likelihood (REML), and restricted maximum a posteriori (REMAP). The scripts use the same one-dimensional interpolation setup so that only the selection criterion changes.

What this example does

The rendered preview uses REMAP, selects covariance parameters, and plots the posterior prediction. The full scripts below expose the corresponding REMAP, REML, and ML procedures. All three procedures optimize a scalar criterion over covparam but differ in how they handle mean parameters and prior regularization.

The method names encode the parameter layout, the criterion, and the prior terms. In select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior, sigma2_rho means [log(sigma2), -log(rho_0), ..., -log(rho_{d-1})]. remap means restricted maximum a posteriori selection. logsigma2_logrho_prior states which prior terms are added to the restricted likelihood.

Mathematical criteria

Let \(P_i\) denote the matrix of mean basis functions evaluated at the observation points, and let \(K_{ii}(\theta)\) be the covariance matrix. For a linear mean \(m_\beta(x)=h(x)^\top\beta\), ML profiles out \(\beta\). The profiled estimator is

\[\widehat\beta_\theta = \left(P_i^\top K_{ii}^{-1}P_i\right)^{-1} P_i^\top K_{ii}^{-1}z_i.\]

Up to constants independent of \(\theta\), the ML criterion is

\[J_{\mathrm{ML}}(\theta) = \frac12 \log |K_{ii}| + \frac12 (z_i-P_i\widehat\beta_\theta)^\top K_{ii}^{-1} (z_i-P_i\widehat\beta_\theta).\]

REML is computed from contrasts that remove the linear mean. Let \(W_i\in\mathbb{R}^{n\times(n-q)}\) have full column rank and satisfy \(W_i^\top P_i=0\). GPmp builds \(W_i\) from a complete QR factorization of \(P_i\). With

\[G_\theta = W_i^\top K_{ii}(\theta) W_i,\]

the REML criterion implemented in GPmp is, up to constants independent of \(\theta\),

\[J_{\mathrm{REML}}(\theta) = \frac12 \log |G_\theta| + \frac12 (W_i^\top z_i)^\top G_\theta^{-1} (W_i^\top z_i).\]

REMAP adds prior regularization on covariance parameters:

\[J_{\mathrm{REMAP}}(\theta) = J_{\mathrm{REML}}(\theta) - \log \pi(\theta).\]

How to interpret the comparison

ML optimizes the likelihood after estimating mean parameters. REML accounts for the degrees of freedom used by the mean and is often preferable when the mean is unknown. See [13]. REMAP adds prior terms to the restricted likelihood. This may stabilize poorly identified variance or lengthscale parameters. Different criteria can lead to different posterior uncertainty, even when the posterior mean looks similar.

Functions used

  • select_parameters_sigma2_rho_with_reml selects parameters by restricted likelihood for the sigma2_rho covariance-parameter convention.

  • select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior adds prior regularization on log(sigma2) and logrho for the same convention.

  • The returned info object records the selected covparam, optimization status, objective history, and criterion callables.

../_images/parameter_selection_0_0.png

REMAP

Script: examples/gpmp_example20_1d_interpolation_variation_remap.py

  1"""
  2Plot and optimize the restricted negative log-likelihood
  3
  4Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
  5Copyright (c) 2022-2026, CentraleSupelec
  6License: GPLv3 (see LICENSE)
  7"""
  8
  9import gpmp.num as gnp
 10import gpmp as gp
 11import matplotlib.pyplot as plt
 12
 13
 14def generate_data():
 15    """
 16    Data generation.
 17
 18    Returns
 19    -------
 20    tuple
 21        (xt, zt): target data
 22        (xi, zi): input dataset
 23    """
 24    dim = 1
 25    nt = 200
 26    box = [[-1], [1]]
 27    xt = gp.misc.designs.regulargrid(dim, nt, box)
 28    zt = gp.misc.testfunctions.twobumps(xt)
 29
 30    ni = 6
 31    xi = gp.misc.designs.ldrandunif(dim, ni, box)
 32    zi = gp.misc.testfunctions.twobumps(xi)
 33
 34    return xt, zt, xi, zi
 35
 36
 37def constant_mean(x, param):
 38    return gnp.ones((x.shape[0], 1))
 39
 40
 41def kernel(x, y, covparam, pairwise=False):
 42    p = 3
 43    return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 44
 45
 46def visualize_results(xt, zt, xi, zi, zpm, zpv):
 47    """
 48    Visualize the results using gp.plot.plotutils (a matplotlib wrapper).
 49
 50    Parameters
 51    ----------
 52    xt : numpy.ndarray
 53        Target x values
 54    zt : numpy.ndarray
 55        Target z values
 56    xi : numpy.ndarray
 57        Input x values
 58    zi : numpy.ndarray
 59        Input z values
 60    zpm : numpy.ndarray
 61        Posterior mean
 62    zpv : numpy.ndarray
 63        Posterior variance
 64    """
 65    fig = gp.plot.Figure(isinteractive=True)
 66    fig.plot(xt, zt, "k", linewidth=1, linestyle=(0, (5, 5)))
 67    fig.plotdata(xi, zi)
 68    fig.plotgp(xt, zpm, zpv, colorscheme="simple")
 69    fig.xylabels("$x$", "$z$")
 70    fig.title("Posterior GP with parameters selected by ReMAP")
 71    fig.show(grid=True, xlim=[-1.0, 1.0], legend=True, legend_fontsize=9)
 72
 73
 74def main():
 75    xt, zt, xi, zi = generate_data()
 76
 77    model = gp.core.Model(constant_mean, kernel)
 78
 79    # Automatic selection of parameters using REMAP
 80    model, info = gp.kernel.select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior(
 81        model, xi, zi, info=True
 82    )
 83    gp.modeldiagnosis.diag(
 84        model, "linear_mean_maternp_anisotropic", info, xi, zi
 85    )
 86
 87    # Prediction
 88    zpm, zpv = model.predict(xi, zi, xt)
 89
 90    # Visualization
 91    print("\nVisualization")
 92    print("-------------")
 93    plot_likelihood = True
 94    if plot_likelihood:
 95        gp.modeldiagnosis.plot_selection_criterion_sigma_rho(
 96            model, info, criterion_name="restricted maximum a posteriori"
 97        )
 98
 99    visualize_results(xt, zt, xi, zi, zpm, zpv)
100
101
102if __name__ == "__main__":
103    main()

REML

Script: examples/gpmp_example21_1d_interpolation_variation_reml.py

  1"""
  2Plot and optimize the restricted negative log-likelihood
  3
  4Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
  5Copyright (c) 2022-2026, CentraleSupelec
  6License: GPLv3 (see LICENSE)
  7"""
  8
  9import gpmp.num as gnp
 10import gpmp as gp
 11import matplotlib.pyplot as plt
 12
 13
 14def generate_data():
 15    """
 16    Data generation.
 17
 18    Returns
 19    -------
 20    tuple
 21        (xt, zt): target data
 22        (xi, zi): input dataset
 23    """
 24    dim = 1
 25    nt = 200
 26    box = [[-1], [1]]
 27    xt = gp.misc.designs.regulargrid(dim, nt, box)
 28    zt = gp.misc.testfunctions.twobumps(xt)
 29
 30    ni = 8
 31    xi = gp.misc.designs.ldrandunif(dim, ni, box)
 32    zi = gp.misc.testfunctions.twobumps(xi)
 33
 34    return xt, zt, xi, zi
 35
 36
 37def constant_mean(x, param):
 38    return gnp.ones((x.shape[0], 1))
 39
 40
 41def kernel(x, y, covparam, pairwise=False):
 42    p = 3
 43    return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 44
 45
 46def visualize_results(xt, zt, xi, zi, zpm, zpv):
 47    """
 48    Visualize the results using gp.plot.plotutils (a matplotlib wrapper).
 49
 50    Parameters
 51    ----------
 52    xt : numpy.ndarray
 53        Target x values
 54    zt : numpy.ndarray
 55        Target z values
 56    xi : numpy.ndarray
 57        Input x values
 58    zi : numpy.ndarray
 59        Input z values
 60    zpm : numpy.ndarray
 61        Posterior mean
 62    zpv : numpy.ndarray
 63        Posterior variance
 64    """
 65    fig = gp.plot.Figure(isinteractive=True)
 66    fig.plot(xt, zt, "k", linewidth=1, linestyle=(0, (5, 5)))
 67    fig.plotdata(xi, zi)
 68    fig.plotgp(xt, zpm, zpv, colorscheme="simple")
 69    fig.xylabels("$x$", "$z$")
 70    fig.title("Posterior GP with parameters selected by ReML")
 71    fig.show(xlim=[-1.0, 1.0])
 72
 73
 74def main():
 75    xt, zt, xi, zi = generate_data()
 76
 77    meanparam0 = None
 78    covparam0 = None
 79    model = gp.core.Model(constant_mean, kernel, meanparam0, covparam0)
 80
 81    # Parameter initial guess
 82    covparam0 = gp.kernel.anisotropic_parameters_initial_guess(model, xi, zi)
 83
 84    # selection criterion
 85    selection_criterion = gp.kernel.negative_log_restricted_likelihood
 86    nlrl, nlrl_pregrad, nlrl_nograd, dnlrl = gp.kernel.make_selection_criterion_with_gradient(
 87        model, selection_criterion, xi, zi
 88    )
 89
 90    # optimize parameters
 91    covparam_reml, info = gp.kernel.autoselect_parameters(
 92        covparam0, nlrl_pregrad, dnlrl, silent=False, info=True
 93    )
 94
 95    model.covparam = gnp.asarray(covparam_reml)
 96    info["covparam0"] = covparam0
 97    info["covparam"] = covparam_reml
 98    info["selection_criterion"] = nlrl
 99
100    gp.modeldiagnosis.diag(
101        model, "linear_mean_maternp_anisotropic", info, xi, zi
102    )
103
104    # Prediction
105    zpm, zpv = model.predict(xi, zi, xt)
106
107    # Visualization
108    print("\nVisualization")
109    print("-------------")
110    visualize_results(xt, zt, xi, zi, zpm, zpv)
111
112    zloom, zloov, eloo = model.loo(xi, zi)
113    gp.plot.plot_loo(zi, zloom, zloov)
114
115
116if __name__ == "__main__":
117    main()

ML

Script: examples/gpmp_example22_1d_interpolation_variation_ml.py

  1"""
  2Plot and optimize the restricted negative log-likelihood
  3
  4Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
  5Copyright (c) 2022-2026, CentraleSupelec
  6License: GPLv3 (see LICENSE)
  7"""
  8
  9import gpmp.num as gnp
 10import gpmp as gp
 11import matplotlib.pyplot as plt
 12
 13
 14def generate_data():
 15    """
 16    Data generation.
 17
 18    Returns
 19    -------
 20    tuple
 21        (xt, zt): target data
 22        (xi, zi): input dataset
 23    """
 24    c = 1.0
 25    dim = 1
 26    nt = 200
 27    box = [[-1], [1]]
 28    xt = gp.misc.designs.regulargrid(dim, nt, box)
 29    zt = gp.misc.testfunctions.twobumps(xt) + c
 30
 31    ni = 8
 32    xi = gp.misc.designs.ldrandunif(dim, ni, box)
 33    zi = gp.misc.testfunctions.twobumps(xi) + c
 34
 35    return xt, zt, xi, zi
 36
 37
 38def constant_mean(x, param):
 39    return param * gnp.ones((x.shape[0], 1))
 40
 41
 42def kernel(x, y, covparam, pairwise=False):
 43    p = 3
 44    return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 45
 46
 47def visualize_results(xt, zt, xi, zi, zpm, zpv):
 48    """
 49    Visualize the results using gp.plot.plotutils (a matplotlib wrapper).
 50
 51    Parameters
 52    ----------
 53    xt : numpy.ndarray
 54        Target x values
 55    zt : numpy.ndarray
 56        Target z values
 57    xi : numpy.ndarray
 58        Input x values
 59    zi : numpy.ndarray
 60        Input z values
 61    zpm : numpy.ndarray
 62        Posterior mean
 63    zpv : numpy.ndarray
 64        Posterior variance
 65    """
 66    fig = gp.plot.Figure(isinteractive=True)
 67    fig.plot(xt, zt, "k", linewidth=1, linestyle=(0, (5, 5)))
 68    fig.plotdata(xi, zi)
 69    fig.plotgp(xt, zpm, zpv, colorscheme="bw")
 70    fig.xylabels("$x$", "$z$")
 71    fig.title("Posterior GP with parameters selected by ML")
 72    fig.show(xlim=[-1.0, 1.0])
 73
 74
 75def main():
 76    xt, zt, xi, zi = generate_data()
 77
 78    meanparam0 = None
 79    covparam0 = None
 80    model = gp.core.Model(
 81        constant_mean, kernel, meanparam0, covparam0, meantype="parameterized"
 82    )
 83
 84    # Parameter initial guess
 85    (
 86        meanparam0,
 87        covparam0,
 88    ) = gp.kernel.anisotropic_parameters_initial_guess_constant_mean(model, xi, zi)
 89
 90    param0 = gnp.concatenate((meanparam0, covparam0))
 91
 92    # selection criterion
 93    nll, nll_pregrad, nll_nograd, dnll = gp.kernel.make_selection_criterion_with_gradient(
 94        model,
 95        gp.kernel.negative_log_likelihood,
 96        xi,
 97        zi,
 98        parameterized_mean=True,
 99        meanparam_len=1,
100    )
101
102    param_ml, info = gp.kernel.autoselect_parameters(
103        param0, nll_pregrad, dnll, silent=False, info=True
104    )
105
106    model.meanparam = gnp.asarray(param_ml[0])
107    model.covparam = gnp.asarray(param_ml[1:])
108
109    info["covparam0"] = param0[1:]
110    info["covparam"] = param_ml[1:]
111    info["selection_criterion"] = nll
112
113    gp.modeldiagnosis.diag(
114        model, "linear_mean_maternp_anisotropic", info, xi, zi
115    )
116
117    # Prediction
118    zpm, zpv = model.predict(xi, zi, xt)
119
120    # Visualization
121    print("\nVisualization")
122    print("-------------")
123    visualize_results(xt, zt, xi, zi, zpm, zpv)
124
125    zloom, zloov, eloo = model.loo(xi, zi)
126    gp.plot.plot_loo(zi, zloom, zloov)
127
128    return model
129
130
131if __name__ == "__main__":
132    model = main()