class gpr::LaplaceINLA

Overview

Laplace / INLA uncertainty quantifier over GP/TP hyperparameters. More…

#include <LaplaceINLA.h>
 
class LaplaceINLA {
public:
    // typedefs
 
    typedef std::function<std::tuple<EigenMatrix, EigenMatrix, Eigen::VectorXd>(const Eigen::VectorXd&theta)> KernelBuilder;
    typedef std::function<double(const Eigen::VectorXd&theta)> LogPrior;
 
    // enums
 
    enum Likelihood;
 
    // methods
 
    void setTheta0(const Eigen::VectorXd& theta_star);
    void setHessian(const Eigen::MatrixXd& H);
    void setLikelihood(Likelihood lik);
    Likelihood likelihood() const;
    void setTPHyperprior(double a, double b);
    void setObservationLikelihood(std::shared_ptr<ObservationLikelihood> ol);
    void buildGrid(double f0 = -1.0, bool include_corners = true);
    void fit(const Eigen::VectorXd& y, KernelBuilder kernel_builder, LogPrior log_prior);
    INLAPrediction predict(const Eigen::VectorXd& y, KernelBuilder kernel_builder) const;
    const std::vector<Eigen::VectorXd>& gridPoints() const;
    const Eigen::VectorXd& gridWeights() const;
    const Eigen::VectorXd& finalWeights() const;
    int p() const;
    int nGridPoints() const;
    const Eigen::VectorXd& theta0() const;
    const Eigen::MatrixXd& lScale() const;
    Eigen::MatrixXd hessianInverse() const;
    Eigen::MatrixXd precisionMatrix() const;
};

Detailed Documentation

Laplace / INLA uncertainty quantifier over GP/TP hyperparameters.

Usage: LaplaceINLA inla; inla.setTheta0(theta_star); // MAP from SCG inla.setHessian(H); // p x p, from assembleExactHessian inla.buildGrid(); // 1 + 2p + 2^p CCD points inla.fit(y, kernel_builder, log_prior); INLAPrediction p = inla.predict(y, kernel_builder);

The KernelBuilder callback signature takes the RAW hyperparameter vector theta directly (caller owns the log-base convention). This differs from BayesianModelAveraging::KernelBuilder, which takes (lengthscale, noise) in linear space; callers who want to share a kernel closure across LaplaceINLA and BMA should write two adapters rather than rely on implicit coercion inside either class.

LOAD-BEARING CONTRACT: the LAST element of theta is always log-noise (natural log). LaplaceINLA extracts it internally for FactorizedLikelihood::evaluateGP / evaluateTP.

Typedefs

typedef std::function<std::tuple<EigenMatrix, EigenMatrix, Eigen::VectorXd>(const Eigen::VectorXd&theta)> KernelBuilder

Builds kernel matrices from the RAW hyperparameter vector theta. CONVENTION: the last element of theta is always log-noise (natural log). Callers own the log-base for the other components. Returns (K_XX, K_test_X, k_test_test_diag).

typedef std::function<double(const Eigen::VectorXd&theta)> LogPrior

Log prior over hyperparameters. Called at every grid point. Return 0.0 for a flat improper prior; any log density otherwise.

Methods

void setTheta0(const Eigen::VectorXd& theta_star)

MAP hyperparameters (log space if using the project’s conventions).

void setHessian(const Eigen::MatrixXd& H)

p x p observed-information Hessian at theta*. Must be positive-definite.

void setTPHyperprior(double a, double b)

TP hyperprior (only used when lik_ == T_CLOSED_FORM).

void setObservationLikelihood(std::shared_ptr<ObservationLikelihood> ol)

Attach a non-Gaussian observation likelihood. Only consulted when lik_ == NON_GAUSSIAN; the branch throws if this was never attached (fail fast rather than silently falling back to Gaussian).

void buildGrid(double f0 = -1.0, bool include_corners = true)

Build the CCD grid in theta-space from (theta*, H).

Places 1 center, 2p axial, and 2^p factorial-corner points at z = 0, +/- f0 * e_i, and (+/- f0/sqrt(p), …, +/- f0/sqrt(p)) in standardised z-space, then maps back to theta-space via theta_k = theta* + U Lambda^{-1/2} z_k.

Grid-rule weights w_k^grid are chosen so the rule integrates z (to zero) and z z^T (to identity) exactly against N(0, I). These match Rue et al. 2009 eq. 6.5 and the moments validated by test_ccd_grid_moments_p2 in python/tests/test_sympy_laplace.py.

Parameters:

f0

Radius in standardised space. Pass a negative value (default) to use Rue’s default sqrt(p + 2).

include_corners

Include the 2^p factorial corner points (default true). When false, only the 1 + 2p axial subset is built; the corner Cholesky factorisations that dominate per-call BMA / LAPLACE-CCD wall are skipped. Trade-off: corners contribute ~25-40% of the total grid weight at p = 4 with f0 = sqrt(6), so axial-only loses some calibrated-uncertainty accuracy but cuts per-call cost from K = 25 to K = 9 (vault note Cheap_Mean_Bounded_INLA_Adaptive_ Quadrature_2026_05.org). Use true for paper-quality variance, false for the production-fast variance gate.

void fit(const Eigen::VectorXd& y, KernelBuilder kernel_builder, LogPrior log_prior)

Evaluate the per-grid-point correction rho_k = exp(-Delta phi_k) and build the final normalised weights.

Needs the kernel builder so it can evaluate the true log p(y|theta_k)

• log p(theta_k) at each grid point. The Rue correction ratio is

  rho_k = [ p(y|theta_k) * p(theta_k) * pi_G(theta*) ] / [ p(y|theta*) * p(theta*) * pi_G(theta_k) ]

where pi_G is the Laplace Gaussian surrogate N(theta*, H^-1). In the linear-Gaussian limit phi is exactly quadratic and rho_k = 1 everywhere (see test_rho_k_is_one_on_quadratic_phi).

Parameters:

y	Training targets.
kernel_builder	Same callback as BMA uses.
log_prior	Log prior over hyperparameters (pass a zero-returning lambda for a flat improper prior).

INLAPrediction predict(const Eigen::VectorXd& y, KernelBuilder kernel_builder) const

Predictive moments at the test points embedded in the kernel builder’s K_test_X / k_test_test outputs.

Returns:

INLAPrediction with law-of-total-variance applied to the CCD mixture over theta_k’s posterior predictive.

const std::vector<Eigen::VectorXd>& gridPoints() const

Grid accessors (for testing).

const Eigen::VectorXd& theta0() const

MAP hyperparameter point set by setTheta0. Used by SVGD initFromCCD.

const Eigen::MatrixXd& lScale() const

Posterior-covariance square root L_scale = U * Lambda^{-1/2} in z-space. The posterior covariance is H^{-1} = L_scale * L_scale^T. Used by the 2nd-order TCOV correction in GaussianProcessRegression::calculatePotentialLaplace.

Eigen::MatrixXd hessianInverse() const

p x p inverse of the observed-information Hessian at theta*. Assembled as L_scale * L_scale^T from the cached eigendecomposition of H; cheaper than calling H_.inverse() and numerically stable (eigenvalues are strictly positive by the PD guard in setHessian). Returns an empty matrix if setHessian was never called.

Eigen::MatrixXd precisionMatrix() const

Effective observed-information precision used to map the CCD grid. This includes any precision-form regularization applied by setHessian().