Gaussian process class

class lsqfitgp.GP(covfun=None, *, solver='chol', checkpos=True, checksym=True, checkfinite=True, checklin=True, posepsfac=1, halfmatrix=False, **kw)

Object that represents a Gaussian process.

A GP is structured like a pair of dictionaries, one for “processes”, and one for “elements”. The processes represent independent Gaussian processes, i.e., infinite-dimensional Normally distributed variables. The elements represent finite-dimensional Normal variables, typically finite subsets of the processes.

The methods to define processes start with “def”, while those to define elements starts with “add”. The basic methods are defproc and addx.

A GP object is immutable. Methods that modify the Gaussian process return a new object which differs only in the requested modification, leaving the original untouched.

Parameters:

covfunKernel, optional

An instance of Kernel representing the covariance kernel of the default process of the GP object. It can be left unspecified.

solverstr, default ‘chol’

The algorithm used to decompose the prior covariance matrix. See decompose for the available solvers.

checkposbool, default True

Raise a LinAlgError if the prior covariance matrix turns out non positive within numerical error.

checksymbool, default True

Check that the prior covariance matrix is symmetric.

checkfinitebool, default True

Check that the prior covariance matrix does not contain infs or nans.

checklinbool, default True

The method addlintransf will check that the given transformation is linear on a random input tensor.

posepsfacnumber, default 1

The threshold used to check if the prior covariance matrix is positive definite is multiplied by this factor.

halfmatrixbool, default False

If checksym=False, compute only half of the covariance matrices by unrolling their lower triangular part as flat arrays. This may actually be a large performance hit if the input arrays have large item size or if the implementation of the kernel takes advantage of non-broadcasted inputs.

**kw

Additional keyword arguments are passed to the solver, see decompose.

Attributes:

DefaultProcess

Key of the default process.

Methods

addx(x[, key, deriv, proc])	Add points where the Gaussian process is evaluated.
addlintransf(transf, keys, key, *[, checklin])	Define a finite linear transformation of the evaluated process.
addtransf(tensors, key, *[, axes])	Apply a linear transformation to already specified process points.
addcov(covblocks[, key, decomps])	Add user-defined prior covariance matrix blocks.
defproc(key[, kernel, deriv])	Add an independent process.
deflintransf(key, transf, procs, *[, deriv, ...])	Define a new process as a linear combination of other processes.
deftransf(key, ops, *[, deriv])	Define a new process as a linear combination of other processes.
deflinop(key, transfname, arg, proc)	Define a new process as the transformation of an existing one.
defderiv(key, deriv, proc)	Define a new process as the derivative of an existing one.
defxtransf(key, transf, proc)	Define a new process by transforming the inputs of another one.
defrescale(key, scalefun, proc)	Define a new process as a rescaling of an existing one.
prior([key, raw])	Return an array or a dictionary of arrays of gvars representing the prior for the Gaussian process.
pred(given[, key, givencov, fromdata, raw, ...])	Compute the posterior.
predfromfit(*args, **kw)	Like pred with fromdata=False.
predfromdata(*args, **kw)	Like pred with fromdata=True.
marginal_likelihood(given[, givencov])	Compute the logarithm of the probability of the data.
decompose(posdefmatrix[, solver])	Decompose a nonnegative definite matrix.

class DefaultProcess

Key of the default process.

addcov(covblocks, key=None, *, decomps=None)

Add user-defined prior covariance matrix blocks.

Covariance matrices defined with addcov represent arbitrary finite-dimensional zero-mean Gaussian variables, assumed independent from all other variables in the GP object.

Parameters:

covblocksarray or dictionary of arrays

If an array: a covariance matrix (or tensor) to be added under key key. If a dictionary: a mapping from pairs of keys to the corresponding covariance matrix blocks. A missing off-diagonal block in the dictionary is interpreted as a matrix of zeros, unless the corresponding transposed block is specified.

keyhashable

If covblocks is an array, the dictionary key under which covblocks is added. Can not be specified if covblocks is a dictionary.

decompsDecomposition or dict of Decompositions

Pre-computed decompositions of (not necessarily all) diagonal blocks, as produced by decompose. The keys are single GP keys and not pairs like in covblocks.

Returns:

gpGP

A new GP object with the applied modifications.

Raises:

KeyError

A key is already used in the GP.

ValueError

covblocks and/or key and decomps are malformed or inconsistent.

TypeError

Wrong type of covblocks or decomps.

addlintransf(transf, keys, key, *, checklin=None)

Define a finite linear transformation of the evaluated process.

Parameters:

transfcallable

A function with signature f(array1, array2, ...) -> array which computes the linear transformation. The function must be jax-traceable, i.e., use jax.numpy instead of numpy.

keyssequence

Keys of parts of the process to be passed as inputs to the transformation.

keyhashable

The key of the newly defined points.

checklinbool

If True (default), check that the given function is linear in its inputs. The default can be overridden at initialization of the GP object. Note that an affine function (x -> a + bx) is not linear.

Returns:

gpGP

A new GP object with the applied modifications.

Raises:

RuntimeError

The transformation seems not to be linear. To disable the linearity check, initialize the GP with checklin=False.

addtransf(tensors, key, *, axes=1)

Apply a linear transformation to already specified process points. The result of the transformation is represented by a new key.

Parameters:

tensorsdict

Dictionary mapping keys of the GP to arrays/scalars. Each array is matrix-multiplied with the process array represented by its key, while scalars are just multiplied. Finally, the keys are summed over.

keyhashable

A new key under which the transformation is placed.

axesint

Number of axes to be summed over for matrix multiplication, referring to trailing axes for tensors in ` tensors``, and to heading axes for process points. Default 1.

Returns:

gpGP

A new GP object with the applied modifications.

Notes

The multiplication between the tensors and the process is done with np.tensordot with, by default, 1-axis contraction. For >2d arrays this is different from numpy’s matrix multiplication, which would act on the second-to-last dimension of the second array.

addx(x, key=None, *, deriv=0, proc=_base.GPBase.DefaultProcess)

Add points where the Gaussian process is evaluated.

The GP object keeps the various x arrays in a dictionary. If x is an array, you have to specify its dictionary key with the key parameter. Otherwise, you can directly pass a dictionary for x.

To specify that on the given x a derivative of the process instead of the process itself should be evaluated, use the parameter deriv.

addx may or may not copy the input arrays.

Parameters:

xarray or dictionary of arrays

The points to be added.

keyhashable

If x is an array, the dictionary key under which x is added. Can not be specified if x is a dictionary.

derivDeriv-like

Derivative specification. A Deriv object or something that can be converted to Deriv.

prochashable

The process to be evaluated on the points. If not specified, use the default process.

Returns:

gpGP

A new GP object with the applied modifications.

classmethod decompose(posdefmatrix, solver='chol', **kw)

Decompose a nonnegative definite matrix.

The decomposition can be used to calculate linear algebra expressions where the (pseudo)inverse of the matrix appears.

Parameters:

posdefmatrixarray

A nonnegative definite nonempty symmetric square matrix. If the array is not square, it must have a shape of the kind (k, n, m, …, k, n, m, …) and is reshaped to (k * n * m * …, k * n * m * …).

solverstr

Algorithm used to decompose the matrix.

‘chol’

Cholesky decomposition after regularizing the matrix with a Gershgorin estimate of the maximum eigenvalue.

**kw

Additional options.

epsrel, epsabspositive float or ‘auto’

Specify the threshold for considering small the eigenvalues:

eps = epsrel * maximum_eigenvalue + epsabs

epsrel=’auto’ sets epsrel = matrix_size * float_epsilon, while epsabs=’auto’ sets epsabs = float_epsilon. Default is epsrel=’auto’, epsabs=0.

Returns:

decompDecomposition

An object representing the decomposition of the matrix. The available methods and properties are (K being the matrix):

matrix():

Return K.

ginv():

Compute K⁻.

ginv_linear(X):

Compute K⁻X.

pinv_bilinear(A, r)

Compute A’K⁺r.

pinv_bilinear_robj(A, r)

Compute A’K⁺r, and r can be an array of arbitrary objects.

ginv_quad(A)

Compute A’K⁻A.

ginv_diagquad(A)

Compute diag(A’K⁻A).

correlate(x)

Compute Zx such that K = ZZ’, Z can be rectangular.

back_correlate(X)

Compute Z’X.

pinv_correlate(x):

Compute Z⁺x.

minus_log_normal_density(r, …)

Compute a Normal density and its derivatives.

eps

The threshold below which eigenvalues are not calculable.

n

Number of rows/columns of K.

m

Number of columns of Z.

Notes

The decomposition operations are JAX-traceable, but they are not meant to be differentiated. The method minus_log_normal_density provides required derivatives with a custom implementation, given the derivatives of the inputs.

defderiv(key, deriv, proc)

Define a new process as the derivative of an existing one.

\[g(x) = \frac{\partial^n}{\partial x^n} f(x)\]

Parameters:

keyhashable

The key of the new process.

derivDeriv-like

Derivation order.

prochashable

The key of the process to be derived.

Returns:

gpGP

A new GP object with the applied modifications.

deflinop(key, transfname, arg, proc)

Define a new process as the transformation of an existing one.

Parameters:

keyhashable

Key for the new process.

transfnamehashable

A transformation recognized by the transf method of the kernel.

arg

A valid argument to the transformation.

prochashable

Key of the process to be transformed.

Returns:

gpGP

A new GP object with the applied modifications.

deflintransf(key, transf, procs, *, deriv=0, checklin=False)

Define a new process as a linear combination of other processes.

Let f_i(x), i = 1, 2, … be already defined processes, and T a linear map from processes to a single process. The new process is

h(x) = T(f_1, f_2, …)(x).

Parameters:

keyhashable

The name that identifies the new process in the GP object.

transfcallable

A function with signature transf(callable, callable, ...) -> callable.

procssequence

The keys of the processes to be passed to the transformation.

derivDeriv-like

The linear combination is derived as specified by this parameter.

checklinbool

If True, check if the transformation is linear. Default False.

Returns:

gpGP

A new GP object with the applied modifications.

Notes

The linearity check may fail if the transformation does nontrivial operations with the inner function input.

defproc(key, kernel=None, *, deriv=0)

Add an independent process.

Parameters:

keyhashable

The name that identifies the process in the GP object.

kernelKernel

A kernel for the process. If None, use the default kernel. The difference between the default process and a process defined with the default kernel is that, although they have the same kernel, they are independent.

derivDeriv-like

Derivatives to take on the process defined by the kernel.

Returns:

gpGP

A new GP object with the applied modifications.

defrescale(key, scalefun, proc)

Define a new process as a rescaling of an existing one.

\[g(x) = s(x)f(x)\]

Parameters:

keyhashable

The key of the new process.

scalefuncallable

A function from the domain of the process to a scalar.

prochashable

The key of the process to be transformed.

Returns:

gpGP

A new GP object with the applied modifications.

deftransf(key, ops, *, deriv=0)

Define a new process as a linear combination of other processes.

Let f_i(x), i = 1, 2, … be already defined processes, and g_i(x) be deterministic functions. The new process is defined as

h(x) = g_1(x) f_1(x) + g_2(x) f_2(x) + …

Parameters:

keyhashable

The name that identifies the new process in the GP object.

opsdict

A dictionary mapping process keys to scalars or scalar functions. The functions must take an argument of the same kind of the domain of the process.

derivDeriv-like

The linear combination is derived as specified by this parameter.

Returns:

gpGP

A new GP object with the applied modifications.

defxtransf(key, transf, proc)

Define a new process by transforming the inputs of another one.

\[g(x) = f(T(x))\]

Parameters:

keyhashable

The key of the new process.

transfcallable

A function mapping the new kind input to the input expected by the transformed process.

prochashable

The key of the process to be transformed.

Returns:

gpGP

A new GP object with the applied modifications.

marginal_likelihood(given, givencov=None, **kw)

Compute the logarithm of the probability of the data.

The probability is computed under the Gaussian prior and Gaussian error model. It is also called marginal likelihood. If \(y\) is the data and \(g\) is the Gaussian process, this is

\[\log \int p(y|g) p(g) \mathrm{d} g.\]

Unlike pred, you can’t compute this with a fit result instead of data. If you used the Gaussian process as latent variable in a fit, use the whole fit to compute the marginal likelihood. E.g. lsqfit always computes the logGBF (it’s the same thing).

The input is an array or dictionary of arrays, given. The contents of given represent the input data.

Parameters:

givendictionary of arrays

The data for some/all of the points in the GP. The arrays can contain either gvars or normal numbers, the latter being equivalent to zero-uncertainty gvars.

givencovdictionary of arrays, optional

Covariance matrix of given. If not specified, the covariance is extracted from given with gvar.evalcov(given).

**kw

Additional keyword arguments are passed to the matrix decomposition.

Returns:

logpscalar

The logarithm of the marginal likelihood.

pred(given, key=None, givencov=None, *, fromdata=None, raw=False, keepcorr=None)

Compute the posterior.

The posterior can be computed either for all points or for a subset, and either directly from data or from a posterior obtained with a fit. The latter case is for when the Gaussian process was used in a fit with other parameters.

The output is a collection of gvars, either an array or a dictionary of arrays. They are properly correlated with gvars returned by prior and with the input data/fit.

The input is a dictionary of arrays, given, with keys corresponding to the keys in the GP as added by addx or addtransf.

Parameters:

givendictionary of arrays

The data or fit result for some/all of the points in the GP. The arrays can contain either gvars or normal numbers, the latter being equivalent to zero-uncertainty gvars.

keyNone, key or list of keys, optional

If None, compute the posterior for all points in the GP (also those used in given). Otherwise only those specified by key.

givencovdictionary of arrays, optional

Covariance matrix of given. If not specified, the covariance is extracted from given with gvar.evalcov(given).

fromdatabool

Mandatory. Specify if the contents of given are data or already a posterior.

rawbool, optional

If True, instead of returning a collection of gvars, return the mean and the covariance. When the mean is a dictionary, the covariance is a dictionary whose keys are pairs of keys of the mean (the same format used by gvar.evalcov). Default False.

keepcorrbool, optional

If True (default), the returned gvars are correlated with the prior and the data/fit. If False, they have the correct covariance between themselves, but are independent from all other preexisting gvars.

Returns:

If raw=False (default):

posteriorarray or dictionary of arrays

A collections of gvars representing the posterior.

If raw=True:

pmeanarray or dictionary of arrays

The mean of the posterior. Equivalent to gvar.mean(posterior).

pcov2D array or dictionary of 2D arrays

The covariance matrix of the posterior. If pmean is a dictionary, the keys of pcov are pairs of keys of pmean. Equivalent to gvar.evalcov(posterior).

predfromdata(*args, **kw)

Like pred with fromdata=True.

predfromfit(*args, **kw)

Like pred with fromdata=False.

prior(key=None, *, raw=False)

Return an array or a dictionary of arrays of gvars representing the prior for the Gaussian process. The returned object is not unique but the gvars stored inside are, so all the correlations are kept between objects returned by different calls to prior.

Calling without arguments returns the complete prior as a dictionary. If you specify key, only the array for the requested key is returned.

Parameters:

keyNone, key or list of keys

Key(s) corresponding to one passed to addx or addtransf. None for all keys.

rawbool

If True, instead of returning a collection of gvars return their covariance matrix as would be returned by gvar.evalcov. Default False.

Returns:

If raw=False (default):

priornp.ndarray or dict

A collection of gvars representing the prior.

If raw=True:

covnp.ndarray or dict

The covariance matrix of the prior.