The GaussianProcess Class

Implementation of a Gaussian Process Emulator.

This class provides a representation of a Gaussian Process Emulator. It contains methods for fitting the GP to a given set of hyperparameters, computing the negative log marginal likelihood plus prior (so negative log posterior) and its derivatives, and making predictions on unseen data. Note that routines to estimate hyperparameters are not included in the class definition, and are instead provided externally to facilitate implementation of high performance versions.

The required arguments to initialize a GP is a set of training data consisting of the inputs and targets for each of those inputs. These must be numpy arrays whose first axis is of the same length. Targets must be a 1D array, while inputs can be 2D or 1D. If inputs is 1D, then it is assumed that the length of the second axis is unity.

Optional arguments are the particular mean function to use (default is zero mean), the covariance kernel to use (default is the squared exponential covariance), a list of prior distributions for each hyperparameter (default is no prior information on any hyperparameters) and the method for handling the nugget parameter. The nugget is additional “noise” that is added to the diagonal of the covariance kernel as a variance. This nugget can represent uncertainty in the target values themselves, or simply be used to stabilize the numerical inversion of the covariance matrix. The nugget can be fixed (a non-negative float), can be found adaptively (by passing the string "adaptive" to make the noise only as large as necessary to successfully invert the covariance matrix), can be fit as a hyperparameter (by passing the string "fit"), or pivoting can be used to ignore any collinear matrix rows and ensure a zero nugget is used (by passing the string "pivot").

The internal emulator structure involves arrays for the inputs, targets, and hyperparameters. Other useful information are the number of training examples n, the number of input parameters D, and the number of hyperparameters n_params. These parameters can be obtained externally by accessing these attributes.

Example:

>>> import numpy as np
>>> from mogp_emulator import GaussianProcess
>>> x = np.array([[1., 2., 3.], [4., 5., 6.]])
>>> y = np.array([4., 6.])
>>> gp = GaussianProcess(x, y)
>>> print(gp)
Gaussian Process with 2 training examples and 3 input variables
>>> gp.n
2
>>> gp.D
3
>>> gp.n_params
5
>>> gp.fit(np.zeros(gp.n_params))
>>> x_predict = np.array([[2., 3., 4.], [7., 8., 9.]])
>>> gp.predict(x_predict)
(array([4.74687618, 6.84934016]), array([0.01639298, 1.05374973]),
array([[8.91363045e-05, 7.18827798e-01, 3.74439445e-16],
       [4.64005897e-06, 3.74191346e-02, 1.94917337e-17]]))
class mogp_emulator.GaussianProcess.GaussianProcess(inputs, targets, mean=None, kernel=<mogp_emulator.Kernel.SquaredExponential object>, priors=None, nugget='adaptive', inputdict={}, use_patsy=True)

Implementation of a Gaussian Process Emulator.

This class provides a representation of a Gaussian Process Emulator. It contains methods for fitting the GP to a given set of hyperparameters, computing the negative log marginal likelihood plus prior (so negative log posterior) and its derivatives, and making predictions on unseen data. Note that routines to estimate hyperparameters are not included in the class definition, and are instead provided externally to facilitate implementation of high performance versions.

The required arguments to initialize a GP is a set of training data consisting of the inputs and targets for each of those inputs. These must be numpy arrays whose first axis is of the same length. Targets must be a 1D array, while inputs can be 2D or 1D. If inputs is 1D, then it is assumed that the length of the second axis is unity.

Optional arguments are the particular mean function to use (default is zero mean), the covariance kernel to use (default is the squared exponential covariance), a list of prior distributions for each hyperparameter (default is no prior information on any hyperparameters) and the method for handling the nugget parameter. The nugget is additional “noise” that is added to the diagonal of the covariance kernel as a variance. This nugget can represent uncertainty in the target values themselves, or simply be used to stabilize the numerical inversion of the covariance matrix. The nugget can be fixed (a non-negative float), can be found adaptively (by passing the string "adaptive" to make the noise only as large as necessary to successfully invert the covariance matrix), can be fit as a hyperparameter (by passing the string "fit"), or pivoting can be used to ignore any collinear matrix rows and ensure a zero nugget is used (by passing the string "pivot").

The internal emulator structure involves arrays for the inputs, targets, and hyperparameters. Other useful information are the number of training examples n, the number of input parameters D, and the number of hyperparameters n_params. These parameters can be obtained externally by accessing these attributes.

Example:

>>> import numpy as np
>>> from mogp_emulator import GaussianProcess
>>> x = np.array([[1., 2., 3.], [4., 5., 6.]])
>>> y = np.array([4., 6.])
>>> gp = GaussianProcess(x, y)
>>> print(gp)
Gaussian Process with 2 training examples and 3 input variables
>>> gp.n
2
>>> gp.D
3
>>> gp.n_params
5
>>> gp.fit(np.zeros(gp.n_params))
>>> x_predict = np.array([[2., 3., 4.], [7., 8., 9.]])
>>> gp.predict(x_predict)
(array([4.74687618, 6.84934016]), array([0.01639298, 1.05374973]),
array([[8.91363045e-05, 7.18827798e-01, 3.74439445e-16],
       [4.64005897e-06, 3.74191346e-02, 1.94917337e-17]]))
__init__(inputs, targets, mean=None, kernel=<mogp_emulator.Kernel.SquaredExponential object>, priors=None, nugget='adaptive', inputdict={}, use_patsy=True)

Create a new GaussianProcess Emulator

Creates a new GaussianProcess Emulator from either the input data and targets to be fit and optionally a mean function, covariance kernel, and nugget parameter/method.

Required arguments are numpy arrays inputs and targets, described below. Additional arguments are mean to specify the mean function (default is None for zero mean), kernel to specify the covariance kernel (default is the squared exponential kernel), priors to indicate prior distributions on the hyperparameters, and nugget to specify how the nugget parameter is handled (see below; default is to fit the nugget adaptively).

inputs is a 2D array-like object holding the input data, whose shape is n by D, where n is the number of training examples to be fit and D is the number of input variables to each simulation. If inputs is a 1D array, it is assumed that D = 1.

targets is the target data to be fit by the emulator, also held in an array-like object. This must be a 1D array of length n.

prior must be a list of length n_params whose elements are either Prior-derived objects or None. Each element is used as the prior for the corresponding parameter (with None indicating an uninformative prior). Passing the empty list or None as this argument (in its entirety) may be used as an abbreviation for a list of n_params where all list elements are None.

nugget controls how additional noise is added to the emulator targets when fitting. This can be specified in several ways. If a string is provided, it can take the values of "adaptive" or "fit", which indicate that the nugget will be chosen in the fitting process. The nugget can also be chosen as part of the fitting process, with three options for nugget handling: "adaptive" means that the nugget will be made only as large as necessary to invert the covariance matrix, "fit" means that the nugget will be treated as a hyperparameter to be optimized, and "pivot" indicates that pivoting will be used to ignore any collinear rows and ensure use of a nugget of zero. Alternatively, a non-negative float can be used to specify a fixed noise level. If no value is specified for the nugget parameter, "adaptive" is the default.

Parameters:
  • inputs – Numpy array holding emulator input parameters. Must be 1D with length n or 2D with shape n by D, where n is the number of training examples and D is the number of input parameters for each output. :type inputs: ndarray
  • targets (ndarray) – Numpy array holding emulator targets. Must be 1D with length n
  • mean (None or MeanFunction) – Mean function to be used (optional, default is None for a zero mean)
  • kernel (Kernel or str) – Covariance kernel to be used (optional, default is Squared Exponential) Can provide either a Kernel object or a string matching the kernel type to be used.
  • priors (list or None) – List of priors to be used. Must be None (default) or an empty list (indicates uninformative priors) or list of length n_params. Any parameter for which you wish to specify an uninformative prior, pass None. Number of parameters is the number of parameters in the mean function plus D + 2 (one correlation length per input plus a covariance scale and a nugget). If the nugget will not be fit, the list can have length n_params - 1.
  • nugget (float or str) – Noise to be added to the diagonal, specified as a string or a float. A non-negative float specifies the noise level explicitly, while a string indicates that the nugget will be found via fitting, either as "adaptive", "fit", or "pivot" (see above for a description). Default is "adaptive".
Returns:

New GaussianProcess instance

Return type:

GaussianProcess

D

Returns number of inputs for the emulator

Returns:Number of inputs for the emulator object
Return type:int
fit(theta)

Fits the emulator and sets the parameters

Pre-calculates the matrices needed to compute the log-likelihood and its derivatives and make subsequent predictions. This is called any time the hyperparameter values are changed in order to ensure that all the information is needed to evaluate the log-likelihood and its derivatives, which are needed when fitting the optimal hyperparameters.

The method computes the mean function and covariance matrix and inverts the covariance matrix using the method specified by the value of nugget_type. The factorized matrix and the product of the inverse with the difference between the targets and the mean are cached for later use, and the negative marginal log-likelihood is also cached. This method has no return value, but it does modify the state of the object.

Parameters:theta (ndarray) – Values of the hyperparameters to use in fitting. Must be a numpy array with length n_params
Returns:None
get_K_matrix()

Returns current value of the covariance matrix as a numpy array. Does not include the nugget parameter, as this is dependent on how the nugget is fit.

inputs

Returns inputs for the emulator as a numpy array

Returns:Emulator inputs, 2D array with shape (n, D)
Return type:ndarray
logpost_deriv(theta)

Calculate the partial derivatives of the negative log-posterior

Calculate the partial derivatives of the negative log-posterior with respect to the hyperparameters. Note that this function is normally used only when fitting the hyperparameters, and it is not needed to make predictions.

During normal use, the logpost_deriv method is called after evaluating the logposterior method. The implementation takes advantage of this by reusing cached results, as the factorized covariance matrix is expensive to compute and is used by the logposterior, logpost_deriv, and logpost_hessian methods. If the function is evaluated with a different set of parameters than was previously used to set the log-posterior, the method calls fit (and subsequently resets the cached information).

Parameters:theta (ndarray) – Value of the hyperparameters. Must be array-like with shape (n_params,)
Returns:partial derivatives of the negative log-posterior with respect to the hyperparameters (array with shape (n_params,))
Return type:ndarray
logpost_hessian(theta)

Calculate the Hessian of the negative log-posterior

Calculate the Hessian of the negative log-posterior with respect to the hyperparameters. Note that this function is normally used only when fitting the hyperparameters, and it is not needed to make predictions. It is also used to estimate an appropriate step size when fitting hyperparameters using the lognormal approximation or MCMC sampling.

When used in an optimization routine, the logpost_hessian method is called after evaluating the logposterior method. The implementation takes advantage of this by storing the inverse of the covariance matrix, which is expensive to compute and is used by the logposterior and logpost_deriv methods as well. If the function is evaluated with a different set of parameters than was previously used to set the log-posterior, the method calls fit to compute the needed information and changes the cached values.

Parameters:theta (ndarray) – Value of the hyperparameters. Must be array-like with shape (n_params,)
Returns:Hessian of the negative log-posterior (array with shape (n_params, n_params))
Return type:ndarray
logposterior(theta)

Calculate the negative log-posterior at a particular value of the hyperparameters

Calculate the negative log-posterior for the given set of parameters. Calling this method sets the parameter values and computes the needed inverse matrices in order to evaluate the log-posterior and its derivatives. In addition to returning the log-posterior value, it stores the current value of the hyperparameters and log-posterior in attributes of the object.

Parameters:theta (ndarray) – Value of the hyperparameters. Must be array-like with shape (n_params,)
Returns:negative log-posterior
Return type:float
n

Returns number of training examples for the emulator

Returns:Number of training examples for the emulator object
Return type:int
n_params

Returns number of hyperparameters

Returns the number of hyperparameters for the emulator. The number depends on the choice of mean function, covariance function, and nugget strategy, and possibly the number of inputs for certain choices of the mean function.

Returns:Number of hyperparameters
Return type:int
nugget

Returns emulator nugget parameter

Returns current value of the nugget parameter. If the nugget is to be selected adaptively, by pivoting, or by fitting the emulator and the nugget has not been fit, returns None.

Returns:Current nugget value, either a float or None
Return type:float or None

The nugget parameter controls how noise is added to the covariance matrix in order to stabilize the inversion or smooth the emulator predictions. If nugget is a non-negative float, then that particular value is used for the nugget. Note that setting this parameter to be zero enforces that the emulator strictly interpolates between points. Alternatively, a string can be provided. A value of "fit" means that the nugget is treated as a hyperparameter, and is the last entry in the theta array. Alternatively, if nugget is set to be "adaptive", the fitting routine will adaptively make the noise parameter as large as is needed to ensure that the emulator can be fit. Finally, pivoting can be selected by setting to "pivot", which will ignore any collinear rows in the covariance matrix.

Internally, this modifies both the way the nugget is chosen (which can be determined via the nugget_type property) and the value itself (the nugget property)

Parameters:nugget (float or str) – Noise to be added to the diagonal, specified as a string or a float. A non-negative float specifies the noise level explicitly, while a string indicates that the nugget will be found via fitting, either as "adaptive", "pivot", or "fit" (see above for a description).
Returns:None
Return type:None
nugget_type

Returns method used to select nugget parameter

Returns a string indicating how the nugget parameter is treated, either "adaptive", "pivot", "fit", or "fixed". This is automatically set when changing the nugget property.

Returns:Current nugget fitting method
Return type:str
predict(testing, unc=True, deriv=True, include_nugget=True)

Make a prediction for a set of input vectors for a single set of hyperparameters

Makes predictions for the emulator on a given set of input vectors. The input vectors must be passed as a (n_predict, D), (n_predict,) or (D,) shaped array-like object, where n_predict is the number of different prediction points under consideration and D is the number of inputs to the emulator. If the prediction inputs array is 1D and D == 1 for the GP instance, then the 1D array must have shape (n_predict,). Otherwise, if the array is 1D it must have shape (D,), and the method assumes n_predict == 1. The prediction is returned as an (n_predict, ) shaped numpy array as the first return value from the method.

Optionally, the emulator can also calculate the variances in the predictions and the derivatives with respect to each input parameter. If the uncertainties are computed, they are returned as the second output from the method as an (n_predict,) shaped numpy array. If the derivatives are computed, they are returned as the third output from the method as an (n_predict, D) shaped numpy array.

The final input to the method determines if the predictive variance should include the nugget or not. For situations where the nugget represents observational error and predictions are estimating the true underlying function, this should be set to False. However, most other cases should probably include the nugget, as the emulator is using it to represent some of the uncertainty in the underlying simulator, so the default value is True.

Parameters:
  • testing (ndarray) – Array-like object holding the points where predictions will be made. Must have shape (n_predict, D) or (D,) (for a single prediction)
  • unc (bool) – (optional) Flag indicating if the uncertainties are to be computed. If False the method returns None in place of the uncertainty array. Default value is True.
  • deriv (bool) – (optional) Flag indicating if the derivatives are to be computed. If False the method returns None in place of the derivative array. Default value is True.
  • include_nugget (bool) – (optional) Flag indicating if the nugget should be included in the predictive variance. Only relevant if unc = True. Default is True.
Returns:

Tuple of numpy arrays holding the predictions, uncertainties, and derivatives, respectively. Predictions and uncertainties have shape (n_predict,) while the derivatives have shape (n_predict, D). If the unc or deriv flags are set to False, then those arrays are replaced by None.

Return type:

tuple

priors

The current list priors used in computing the log posterior

To set the priors, must be a list or None. Entries can be None or a subclass of Prior. None indicates weak prior information. An empty list or None means all uninformative priors. Otherwise list should have the same length as the number of hyperparameters, or alternatively can be one shorter than the number of hyperparameters if nugget_type is "adaptive", "pivot" or "fixed" meaning that the nugget hyperparameter is not fit but is instead fixed or found adaptively. If the nugget hyperparameter is not fit, the prior for the nugget will automatically be set to None even if a distribution is provided.

targets

Returns targets for the emulator as a numpy array

Returns:Emulator targets, 1D array with shape (n,)
Return type:ndarray
theta

Returns emulator hyperparameters

Returns current hyperparameters for the emulator as a numpy array if they have been fit. If no parameters have been fit, returns None. Note that the number of parameters depends on the mean function, so the length of this array will vary across instances.

Returns:Current parameter values (numpy array of length n_params), or None if the parameters have not been fit.
Return type:ndarray or None

When set, pre-calculates the matrices needed to compute the log-likelihood and its derivatives and make subsequent predictions. This is called any time the hyperparameter values are changed in order to ensure that all the information is needed to evaluate the log-likelihood and its derivatives, which are needed when fitting the optimal hyperparameters.

The method computes the mean function and covariance matrix and inverts the covariance matrix using the method specified by the value of nugget. The factorized matrix and the product of the inverse with the difference between the targets and the mean are cached for later use, and the negative marginal log-likelihood is also cached. This method has no return value, but it does modify the state of the object.

Parameters:theta (ndarray) – Values of the hyperparameters to use in fitting. Must be a numpy array with length n_params

The PredictResult Class

class mogp_emulator.GaussianProcess.PredictResult

Prediction results object

Dictionary-like object containing mean, uncertainty (variance), and derivatives with respect to the inputs of an emulator prediction. Values can be accessed like a dictionary with keys 'mean', 'unc', and 'deriv' (or indices 0, 1, and 2 for the mean, uncertainty, and derivative for backwards compatability), or using attributes (p.mean if p is an instance of PredictResult). Also supports iteration and unpacking with the ordering (mean, unc, deriv) to be consistent with indexing behavior.

Code is mostly based on scipy’s OptimizeResult class, with some additional code to support iteration and integer indexing.

Variables:
  • mean – Predicted mean for each input point. Numpy array with shape (n_predict,)
  • unc – Predicted variance for each input point. Numpy array with shape (n_predict,)
  • deriv – Predicted derivative with respect to the inputs for each input point. Numpy array with shape (n_predict, D)