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='SquaredExponential', 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]]))

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, the product of the inverse with the difference between the targets and the mean are cached for later use, a second inverted matrix needed to compute the mean function, and the negative marginal log-likelihood are all also cached. This method has no return value, but it does modify the state of the object.

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

get_K_matrix()

Returns current value of the covariance matrix

Computes the covariance matrix (covariance of the inputs with respect to themselves) as a numpy array. Does not include the nugget parameter, as this is dependent on how the nugget is fit.

Returns:Covariance matrix conditioned on the emulator inputs. Will be a numpy array with shape (n,n). Return type:ndarray

get_cov_matrix(other_inputs)

Computes the covariance matrix for a set of inputs

Compute the covariance matrix for the emulator. Assumes the second set of inputs is the inputs on which the GP is conditioned. Thus, calling with the inputs returns the covariance matrix, while calling with a different set of values gives the information needed to make predictions. Note that this does not include the nugget, which (if relevant) is computed separately.

Parameters:otherinputs (ndarray) – Input values for which the covariance is desired. Must be a 2D numpy array with the second dimension matching D for this emulator. Returns:Covariance matrix for the provided inputs relative to the emulator inputs. If the other_inputs array has first dimension of length M, then the returned array will have shape (n,M) Return type:

get_design_matrix(inputs)

Returns the design matrix for a set of inputs

For a given set of inputs, compute the design matrix based on the GP mean function.

Parameters:inputs – 2D numpy array for input values to be used in computing the design matrix. Second dimension must match the number of dimensions of the input data (D).

has_nugget

Boolean indicating if the GP has a nugget parameter

Returns:Boolean indicating if GP has a nugget Return type:bool

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_data,) Returns:partial derivatives of the negative log-posterior with respect to the hyperparameters (array with shape (n_data,)) Return type:ndarray

logpost_hessian(theta)

Calculate the Hessian of the negative log-posterior

NOTE: NOT CURRENTLY SUPPORTED

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_data,) 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_corr

Returns number of correlation length parameters

Returns:Number of correlation length parameters Return type:int

n_mean

Returns number of mean parameters

Returns:Number of mean parameters Return type:int

n_params

Returns number of hyperparameters

Returns the number of fitting hyperparameters for the emulator. The number depends on the choice of covariance function, nugget strategy, and possibly the number of inputs depending on the covariance function. Note that this does not include the mean function, which is fit analytically.

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

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:Nugget fitting method Return type:str

predict(testing, unc=True, deriv=False, include_nugget=True, full_cov=False)

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. If the uncertainties are computed, they are returned as the second output from the method, either as an (n_predict,) shaped numpy array if full_cov=False or an array with (n_predict, n_predict) if the full covariance is computed (default is only to compute the variance, which is the diagonal of the full covariance).

Derivatives have been deprecated due to changes in how the mean function is computed, so setting deriv=True will have no effect and will raise a DeprecationWarning.

The include_nugget kwarg 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.
• include_nugget (bool) – (optional) Flag indicating if the nugget should be included in the predictive variance. Only relevant if unc = True. Default is True.
• full_cov (bool) – (optional) Flag indicating if the full covariance should be computed for the uncertainty. Only relevant if unc = True. Default is False.

Returns:

PredictResult object holding numpy arrays with the predictions and uncertainties. Predictions and uncertainties have shape (n_predict,). If the unc or flag is set to False, then the unc array is replaced by None.

Return type:

PredictResult

priors

Specifies the priors using a GPPriors object. Property returns GPPriors object. Can be set with either a GPPriors object, a dictionary holding the arguments needed to construct a GPPriors object, or None to use the default priors.

If None is provided, use default priors, which are weak prior information for the mean function, Inverse Gamma distributions for the correlation lengths and nugget (if the nugget is fit), and weak prior information for the covariance. The Inverse Gamma distributions are chosen to put 99% of the distribution mass between the minimum and maximum spacing of the inputs. More fine-grained control of the default priors can be obtained by using the class method of the GPPriors object.

Note that the default priors are not weak priors – to use weak prior information, the user must explicitly construct a GPPriors object with the appropriate number of parameters and nugget type.

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 GPParams object. If no parameters have been fit, the data for that object will be set to None. Note that the number of parameters depends on the mean function and whether or not a nugget is fit, so the length of this array will vary across instances.

Returns:Current parameter values as a GPParams object. If the emulator has not been fit, the parameter values will not have been set. Return type:GPParams

When parameter values have been 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 (GPParams) – Values of the hyperparameters to use in fitting. Must be a GPParams object, a numpy array with length n_params, or None to specify no parameters.

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)