Marginal Gaussian Process (MGP)

Marginal Gaussian Processes (MGP) are Gaussian Processes taking into account the uncertainty of the hyperparameters defined as a density probability function. Especially we suppose that the function to model \(f : \Omega \mapsto \mathbb{R}\), where \(\Omega \subset \mathbb{R}^d\) and \(d\) is the number of design variables, lies in a linear embedding \(\mathcal{A}\) such as \(\mathcal{A} = \{ u = Ax, x\in\Omega\}\), \(A \in \mathbb{R}^{d \times d_e}\) and \(f(x)=f_{\mathcal{A}}(Ax)\) with \(f(x)=f_{\mathcal{A}} : \mathcal{A} \mapsto \mathbb{R}\) et \(d_e \ll d\).

Then, we must use a kernel \(k(x,x')=k_{\mathcal{A}}(Ax,Ax')\) whose each component of the transfert matrix \(A\) is an hyperparameter. Thus we have \(d_e \times d\) hyperparameters to find. (Note that \(d_e\) is defined as \(n\_comp\) in the code).

Moreover, we suppose that \(A\) follows a normal distribution \(\mathcal{N}(vect(\hat A),\Sigma)\), with a mean vector \(vect(\hat A)\) and a covariance matrix \(\Sigma\), which are taken into account in the Gaussian Process building. Using a Gaussian Process with \(\hat A\) as hyperparameters and some approximations, we obtain the prediction and variance of the MGP considering the normal distribution of \(A\). For more theory on the MGP, refer to [1].

Limitations

This implementation only considers a Gaussian kernel. Moreover, the kernel \(k_{\mathcal{A}}\) is a uniform Gaussian kernel. The training of the MGP can be very time consuming.

Usage

import matplotlib.pyplot as plt
import numpy as np

from smt.sampling_methods import LHS
from smt.surrogate_models import MGP

# Construction of the DOE
dim = 3

def fun(x):
    import numpy as np

    res = (
        np.sum(x, axis=1) ** 2
        - np.sum(x, axis=1)
        + 0.2 * (np.sum(x, axis=1) * 1.2) ** 3
    )
    return res

sampling = LHS(
    xlimits=np.asarray([(-1, 1)] * dim), criterion="m", random_state=42
)
xt = sampling(8)
yt = np.atleast_2d(fun(xt)).T

# Build the MGP model
sm = MGP(
    theta0=[1e-2],
    print_prediction=False,
    n_comp=1,
)
sm.set_training_values(xt, yt[:, 0])
sm.train()

# Get the transfert matrix A
emb = sm.embedding["C"]

# Compute the smallest box containing all points of A
upper = np.sum(np.abs(emb), axis=0).item()
lower = -upper

# Test the model
u_plot = np.atleast_2d(np.arange(lower, upper, 0.01)).T
x_plot = sm.get_x_from_u(u_plot)  # Get corresponding points in Omega
y_plot_true = fun(x_plot)
y_plot_pred = sm.predict_values(u_plot)
sigma_MGP = sm.predict_variances(u_plot)
sigma_KRG = sm.predict_variances_no_uq(u_plot)

u_train = sm.get_u_from_x(xt)  # Get corresponding points in A

# Plots
fig, ax = plt.subplots()
ax.plot(u_plot, y_plot_pred, label="Predicted")
ax.plot(u_plot, y_plot_true, "k--", label="True")
ax.plot(u_train, yt, "k+", mew=3, ms=10, label="Train")
ax.fill_between(
    u_plot[:, 0],
    y_plot_pred[:, 0] - 3 * sigma_MGP[:, 0],
    y_plot_pred[:, 0] + 3 * sigma_MGP[:, 0],
    color="r",
    alpha=0.5,
    label="Variance with hyperparameters uncertainty",
)
ax.fill_between(
    u_plot[:, 0],
    y_plot_pred[:, 0] - 3 * sigma_KRG[:, 0],
    y_plot_pred[:, 0] + 3 * sigma_KRG[:, 0],
    color="b",
    alpha=0.5,
    label="Variance without hyperparameters uncertainty",
)

ax.set(xlabel="x", ylabel="y", title="MGP")
fig.legend(loc="upper center", ncol=2)
fig.tight_layout()
fig.subplots_adjust(top=0.74)
plt.show()
___________________________________________________________________________

                                    MGP
___________________________________________________________________________

 Problem size

      # training points.        : 8

___________________________________________________________________________

 Training

   Training ...
   Training - done. Time (sec):  0.5436754
../../_images/mgp_Test_test_mgp.png

Options

List of options

Option

Default

Acceptable values

Acceptable types

Description

print_global

True

None

[‘bool’]

Global print toggle. If False, all printing is suppressed

print_training

True

None

[‘bool’]

Whether to print training information

print_prediction

True

None

[‘bool’]

Whether to print prediction information

print_problem

True

None

[‘bool’]

Whether to print problem information

print_solver

True

None

[‘bool’]

Whether to print solver information

poly

constant

[‘constant’, ‘linear’, ‘quadratic’]

[‘str’]

Regression function type

corr

squar_exp

[‘pow_exp’, ‘abs_exp’, ‘squar_exp’, ‘act_exp’, ‘matern52’, ‘matern32’]

None

Correlation function type

pow_exp_power

1.9

None

[‘float’]

Power for the pow_exp kernel function (valid values in (0.0, 2.0]). This option is set automatically when corr option is squar, abs, or matern.

categorical_kernel

MixIntKernelType.CONT_RELAX

[<MixIntKernelType.CONT_RELAX: ‘CONT_RELAX’>, <MixIntKernelType.GOWER: ‘GOWER’>, <MixIntKernelType.EXP_HOMO_HSPHERE: ‘EXP_HOMO_HSPHERE’>, <MixIntKernelType.HOMO_HSPHERE: ‘HOMO_HSPHERE’>, <MixIntKernelType.COMPOUND_SYMMETRY: ‘COMPOUND_SYMMETRY’>]

None

The kernel to use for categorical inputs. Only for non continuous Kriging

hierarchical_kernel

MixHrcKernelType.ALG_KERNEL

[<MixHrcKernelType.ALG_KERNEL: ‘ALG_KERNEL’>, <MixHrcKernelType.ARC_KERNEL: ‘ARC_KERNEL’>]

None

The kernel to use for mixed hierarchical inputs. Only for non continuous Kriging

nugget

2.220446049250313e-14

None

[‘float’]

a jitter for numerical stability

theta0

[0.01]

None

[‘list’, ‘ndarray’]

Initial hyperparameters

theta_bounds

[1e-06, 20.0]

None

[‘list’, ‘ndarray’]

bounds for hyperparameters

hyper_opt

TNC

[‘Cobyla’, ‘TNC’]

[‘str’]

Optimiser for hyperparameters optimisation

eval_noise

False

[True, False]

[‘bool’]

noise evaluation flag

noise0

[0.0]

None

[‘list’, ‘ndarray’]

Initial noise hyperparameters

noise_bounds

[2.220446049250313e-14, 10000000000.0]

None

[‘list’, ‘ndarray’]

bounds for noise hyperparameters

use_het_noise

False

[True, False]

[‘bool’]

heteroscedastic noise evaluation flag

n_start

10

None

[‘int’]

number of optimizer runs (multistart method)

xlimits

None

None

[‘list’, ‘ndarray’]

definition of a design space of float (continuous) variables: array-like of size nx x 2 (lower, upper bounds)

design_space

None

None

[‘BaseDesignSpace’, ‘list’, ‘ndarray’]

definition of the (hierarchical) design space: use smt.utils.design_space.DesignSpace as the main API. Also accepts list of float variable bounds

random_state

41

None

[‘NoneType’, ‘int’, ‘RandomState’]

Numpy RandomState object or seed number which controls random draws for internal optim (set by default to get reproductibility)

n_comp

1

None

[‘int’]

Number of active dimensions

prior

{‘mean’: [0.0], ‘var’: 1.25}

None

[‘dict’]

Parameters for Gaussian prior of the Hyperparameters