pymc.gp.HSGP#

class pymc.gp.HSGP(m, L=None, c=None, drop_first=False, parametrization='noncentered', *, boundary='dirichlet', mean_func=<pymc.gp.mean.Zero object>, cov_func)[source]#

Hilbert Space Gaussian process approximation.

The gp.HSGP class is an implementation of the Hilbert Space Gaussian process. It is a reduced rank GP approximation that uses a fixed set of basis vectors whose coefficients are random functions of a stationary covariance function’s power spectral density. Its usage is largely similar to gp.Latent. Like gp.Latent, it does not assume a Gaussian noise model and can be used with any likelihood, or as a component anywhere within a model. Also like gp.Latent, it has prior and conditional methods. It supports any sum of covariance functions that implement a power_spectral_density method. (Note, this excludes the Periodic covariance function, which uses a different set of basis functions for a low rank approximation, as described in HSGPPeriodic.).

For information on choosing appropriate m, L, and c, refer to Ruitort-Mayol et al. or to the PyMC examples that use HSGP.

To work with the HSGP in its “linearized” form, as a matrix of basis vectors and a vector of coefficients, see the method prior_linearized.

Parameters:
m: list

The number of basis vectors to use for each active dimension (covariance parameter active_dim).

L: list

The boundary of the space for each active_dim. Choose L such that the domain [-L, L] contains all points in the column of X given by the active_dim.

c: float

The proportion extension factor. Used to construct L from X. Defined as S = max|X| such that X is in [-S, S]. L is calculated as c * S. One of c or L must be provided. Further information can be found in Ruitort-Mayol et al.

drop_first: bool

Default False. Sometimes the first basis vector is quite “flat” and very similar to the intercept term. When there is an intercept in the model, ignoring the first basis vector may improve sampling. With boundary="neumann" the first basis vector is exactly constant (its prior variance is the spectral density at zero frequency divided by the volume of the box); with the mixed conditions it is not constant. This argument will be deprecated in future versions.

parametrization: str

Whether to use the centered or noncentered parametrization when multiplying the basis by the coefficients.

boundary: str or sequence of str, default “dirichlet”

Boundary condition of the Laplace eigenbasis at the two ends of the approximation box. A single string applies to every active dimension; a sequence with one entry per active dimension sets the condition dimension by dimension, e.g. ("neumann", "dirichlet") for zero slope at both ends of the first dimension and zero value at both ends of the second. Each name reads as "<at lower end>-<at upper end>":

  • "dirichlet": sine basis, the approximate GP is pinned to zero at both ends and its prior variance shrinks to zero there.

  • "neumann": cosine basis (including the constant), the approximate GP has zero slope at both ends and its prior variance is doubled there.

  • "dirichlet-neumann" / "neumann-dirichlet": zero value at one end and zero slope at the other.

For ExpQuad and Matérn kernels all conditions deviate from the exact stationary GP only within roughly 1.5 lengthscales of an end, and by the same amount, so the guidance for m and c is the same for all of them (heavy-tailed kernels such as RatQuad reach further into the box, for every condition). A boundary condition is a modeling assumption about the edge, not a better approximation. Use a Neumann end where the function is known to have zero slope (symmetry axis, zero flux, plateau) and a Dirichlet end where it is known to vanish (e.g. a radial profile decaying to background); place the end of the box on that point (see the note on the box below). A mismatched condition is worse than the default. When forecasting a short horizon just past the data with a small c, a Neumann end keeps the forecast uncertainty where Dirichlet collapses it to zero; over horizons of a lengthscale or more, or to approximate the unconstrained GP, increase c instead (approx_hsgp_hyperparams with the prediction range included in x_range).

"neumann" evaluates the spectral density at frequency zero. For RatQuad this is finite only for alpha > input_dim / 2 and becomes very large as alpha approaches that value; at or below it the constant basis vector gets a nan or inf coefficient (also mid-chain, if alpha is a free parameter). Use drop_first=True or another boundary condition in that case.

cov_func: Covariance function, must be an instance of `Stationary` and implement a

power_spectral_density method.

mean_func: None, instance of Mean

The mean function. Defaults to zero.

Notes

The approximation box is [center - L, center + L] per active dimension, where center is the midpoint of the X first passed to prior or prior_linearized and L is either given or c times the half-range of that X. To place the ends of the box on physical boundaries (required for the mixed conditions to mean what you intend), make X span the physical domain: build the basis on a grid covering the domain with prior_linearized and index the observed rows in the likelihood, e.g.

r_grid = np.linspace(0.0, R, 200)[:, None]  # whole physical domain
with pm.Model():
    gp = pm.gp.HSGP(m=[30], L=[R / 2], boundary="neumann-dirichlet", cov_func=cov_func)
    phi, sqrt_psd = gp.prior_linearized(r_grid)
    beta = pm.Normal("beta", size=gp.n_basis_vectors)
    f = pm.Deterministic("f", phi @ (beta * sqrt_psd))
    pm.Normal("y", mu=f[observed_idx], sigma=sigma, observed=y)

Inputs outside the box, for prior or conditional, are not checked: the basis simply reflects, and the result is not a GP prediction.

References

  • Ruitort-Mayol, G., and Anderson, M., and Solin, A., and Vehtari, A. (2022). Practical Hilbert Space Approximate Bayesian Gaussian Processes for Probabilistic Programming

  • Solin, A., Sarkka, S. (2019) Hilbert Space Methods for Reduced-Rank Gaussian Process Regression.

Examples

# A three dimensional column vector of inputs.
X = np.random.rand(100, 3)

with pm.Model() as model:
    # Specify the covariance function.
    # Three input dimensions, but we only want to use the last two.
    cov_func = pm.gp.cov.ExpQuad(3, ls=0.1, active_dims=[1, 2])

    # Specify the HSGP.
    # Use 25 basis vectors across each active dimension for a total of 25 * 25 = 625.
    # The value `c = 4` means the boundary of the approximation
    # lies at four times the half width of the data.
    # In this example the data lie between zero and one,
    # so the boundaries occur at -1.5 and 2.5.  The data, both for
    # training and prediction should reside well within that boundary..
    gp = pm.gp.HSGP(m=[25, 25], c=4.0, cov_func=cov_func)

    # If the function is known to have zero slope at the edges of its domain
    # (symmetry, zero flux), encode that instead of extending the box:
    # gp = pm.gp.HSGP(m=[25, 25], c=1.0, boundary="neumann", cov_func=cov_func)

    # Place a GP prior over the function f.
    f = gp.prior("f", X=X)

...

# After fitting or sampling, specify the distribution
# at new points with .conditional
Xnew = np.linspace(-1, 2, 50)[:, None]

with model:
    fcond = gp.conditional("fcond", Xnew=Xnew)

Methods

HSGP.__init__(m[, L, c, drop_first, ...])

HSGP.conditional(name, Xnew[, dims])

Return the (approximate) conditional distribution evaluated over new input locations Xnew.

HSGP.marginal_likelihood(name, X, *args, ...)

HSGP.predict(Xnew[, point, given, diag, model])

HSGP.prior(name, X[, dims, hsgp_coeffs_dims])

Return the (approximate) GP prior distribution evaluated over the input locations X.

HSGP.prior_linearized(X)

Linearized version of the HSGP.

Attributes

L