class gpr::GaussianProcessRegression

Overview

Gaussian Process Regression algorithm. More…

#include <GaussianProcessRegression.h>
 
class GaussianProcessRegression {
public:
    // structs
 
    struct LaplacePosteriorFactor;
    struct LaplacePosteriorSignature;
    struct LaplacePosteriorStamp;
    struct LaplacePosteriorState;
    struct LaplaceTrainingCovarianceGeometry;
    struct SliceLaplaceState;
    struct farthest_point_sampling_t;
 
    // fields
 
    bool failedOptimizer;
    std::unique_ptr<SexpatCF> sexpat_cov_func;
 
    // construction
 
    GaussianProcessRegression();
    virtual ~GaussianProcessRegression();
 
    // methods
 
    void initialize(const InputParameters& parameters, const AtomsConfiguration& conf_info);
    void setHyperparameters(const Observation& all_obs, const AtomsConfiguration& conf_info, const bool update_sexpat_cf_param = true, const bool update_const_cf_param = true, const bool update_sqrt_prior_param = true);
    void setUpDefault();
    void evaluateTrainingCovarianceMatrix(const EigenMatrix& x, const Eigen::VectorXd& x_ind, EigenMatrix& cov_matrix, bool serial_per_rank = false);
    void evaluateCovarianceMatrix(const EigenMatrix& x1, const EigenMatrix& x2, const Eigen::VectorXd& x1_ind, const Eigen::VectorXd& x2_ind, EigenMatrix& C, bool serial_per_rank = false) const;
    void extractUniqueIndices(const Eigen::VectorXd& input_ind, Eigen::VectorXd& output_ind) const;
    void evaluateEnergyAndGradient(const Eigen::VectorXd& w, const EigenMatrix& x, const Eigen::VectorXd& x_ind, const Eigen::VectorXd& y, EnergyAndGradient& energy_and_gradient);
    void optimize(const Observation& observation);
    void optimizeFactorized(const Observation& observation);
    void calculatePotential(Observation& image1);
    void calculatePotentialLaplace(Observation& image1);
    void calculatePotentialLaplaceSlice(Observation& image1);
    void calculateVarianceLaplaceSlice(Observation& image1);
    void setUncertaintyMode(UncertaintyMode mode);
    UncertaintyMode getUncertaintyMode() const;
    void setUsePsisLaplace(bool use);
    bool getUsePsisLaplace() const;
    void setUsePsisDeterministicGrid(bool use);
    bool getUsePsisDeterministicGrid() const;
    void setUseDeltaMethod(bool use);
    bool getUseDeltaMethod() const;
    void buildRFFModel(const Observation& training_data, int D_rff);
    void calculatePotentialRFF(Observation& image);
    void calculateVarianceRFF(Observation& image);
    int getRFFFeatures() const;
    bool isRFFBuilt() const;
    void setRFFFeaturesRequested(int D_rff);
    int getRFFFeaturesRequested() const;
    void setHyperparameterOptimizationMode(HyperparameterOptimizationMode mode);
    HyperparameterOptimizationMode getHyperparameterOptimizationMode() const;
    void setScgMaxIterationsForTesting(Index_t max_iter);
    void resetLaplacePosteriorSignatureBuildsForTesting() const;
    std::size_t laplacePosteriorSignatureBuildsForTesting() const;
    void setForceFactorOwnerLaplaceMeanForTesting(bool value);
    void resetLaplaceMeanBlockPartitionCallsForTesting() const;
    std::size_t laplaceMeanBlockPartitionCallsForTesting() const;
    void calculatePotentialDispatched(Observation& image1);
    void calculateVarianceLaplace(Observation& image1);
    void calculateVarianceDispatched(Observation& image1);
    bool hasBeenOptimized() const;
    void setJitterSigma2(const Index_t value);
    void setParameters(const GPRSetup& parameters);
    LikGaussian& getLikGaussian();
    const LikGaussian& getLikGaussian() const;
    SexpatCF& getSexpAtCovarianceFunction();
    const SexpatCF& getSexpAtCovarianceFunction() const;
    ConstantCF& getConstantCovarianceFunction();
    const ConstantCF& getConstantCovarianceFunction() const;
    Index_t getNumberOfPotentialCalls() const;
    const EigenMatrix& getL() const;
    const Eigen::VectorXd& getAlpha() const;
    const EigenMatrix& getRMatrix() const;
    const Eigen::VectorXd& getRIndices() const;
    bool hasTrainingPermutation() const;
    void assembleDirectionalCovarianceDerivatives(const EigenMatrix& x, const Eigen::VectorXd& ind_Ddim, const Eigen::VectorXd& uDdim, const Eigen::VectorXd& direction, EigenMatrix& first, EigenMatrix& second, bool serial_per_rank = false) const;
    void extractCoordinatesByIndex(const EigenMatrix& x, const Eigen::VectorXd& ind_Ddim, const Index_t ind, FieldMatrixd& x_loc) const;
    void assignBlockToMatrix(const Eigen::VectorXd& ind_Ddim1, const Eigen::VectorXd& ind_Ddim2, const Index_t row_val, const Index_t col_val, const FieldMatrixd& field, EigenMatrix& matrix, bool transpose_field = false) const;
    void accumulateBlockContribution(const Eigen::VectorXd& ind_Ddim1, const Eigen::VectorXd& ind_Ddim2, const Index_t row_val, const Index_t col_val, const FieldMatrixd& field, const EigenMatrix& invC, const Eigen::VectorXd& b, double& trace_accum, double& quad_accum, bool transpose_field = false);
    void removeColumn(const Index_t column_ind, EigenMatrix& matrix);
    Index_t getNumberOfRepetitiveIndices(const Eigen::VectorXd& ind_Ddim) const;
    virtual void calculateVariance(Observation& image1);
    virtual void refreshVariancePrefactor(const Eigen::VectorXd&);
    virtual double localVariancePrefactor(const Eigen::VectorXd&, const Eigen::VectorXd&) const;
    virtual double localVariancePrefactorFromQuadratic(const Eigen::VectorXd&, double) const;
    bool evaluateDirectionalCurvature(const Eigen::VectorXd& x_ind, const EigenMatrix& x, const Eigen::VectorXd& y, const Eigen::VectorXd& direction, double& curvature);
    bool assembleExactHessian(const Eigen::VectorXd& x_ind, const EigenMatrix& x, const Eigen::VectorXd& y, int n_hyperparameters, Eigen::MatrixXd& H);
    bool assembleExactHessianLocalDense(const Eigen::VectorXd& x_ind, const EigenMatrix& x, const Eigen::VectorXd& y, int n_hyperparameters, Eigen::MatrixXd& H);
    void updateModelWithFullData(const Observation& full_observation);
    void updateConfiguration(const AtomsConfiguration& conf_info);
    const AtomsConfiguration& getAtomsConfig() const;
    void calculatePotentialSVGD(Observation& image1);
    void calculateVarianceSVGD(Observation& image1);
    void calculatePotentialBMA(Observation& image1);
    void calculateVarianceBMA(Observation& image1);
    void evaluateEnergyAndGradientLocal(const Eigen::VectorXd& w, const EigenMatrix& x, const Eigen::VectorXd& x_ind, const Eigen::VectorXd& y, EnergyAndGradient& energy_and_gradient);
    void setSVGDConfig(const laplace::SteinVariationalConfig& cfg);
    const laplace::SteinVariationalConfig& getSVGDConfig() const;
};

Detailed Documentation

Gaussian Process Regression algorithm.

Methods

void initialize(const InputParameters& parameters, const AtomsConfiguration& conf_info)

Initialize GPR class.

Parameters:

conf_info	Information about the configurations necessary for the GPR model
parameters	Structure of parameters

void setHyperparameters(const Observation& all_obs, const AtomsConfiguration& conf_info, const bool update_sexpat_cf_param = true, const bool update_const_cf_param = true, const bool update_sqrt_prior_param = true)

Set up hyperparameters for GPR.

Parameters:

all_obs	Structure of all observations
conf_info	Structure of atoms configuration
update_sexpat_cf_param	Set to “true” to update parameters of the SexpAt covariance function
update_const_cf_param	Set to “true” to update parameters of the Constant covariance function
update_sqrt_prior_param	Set to “true” to update parameters of the SQRT prior

void setUpDefault()

Set up default parameters.

This function allocates memory for the private pointers and calls for setUpDefault() method on each pointer.

void evaluateTrainingCovarianceMatrix(const EigenMatrix& x, const Eigen::VectorXd& x_ind, EigenMatrix& cov_matrix, bool serial_per_rank = false)

Evaluate training covariance matrix (gp_cov + noise covariance).

The function takes in matrix x that contains training input vectors to GPR. Returns the noisy covariance matrix cov_matrix for observations y, which is sum of noisless covariance matrix and diagonal term, for example, from Gaussian noise.

Every element ij of noisless covariance matrix contains covariance between inputs i and j in x. The noisless covariance matrix is formed with all covariance functions.

This is essentially the application of the kernel function to the training points

Parameters:

x	Training input vectors to GP
cov_matrix	Noisy covariance matrix for observations y

void evaluateCovarianceMatrix(const EigenMatrix& x1, const EigenMatrix& x2, const Eigen::VectorXd& x1_ind, const Eigen::VectorXd& x2_ind, EigenMatrix& C, bool serial_per_rank = false) const

Evaluate covariance matrix between two input vectors.

Function takes in Gaussian process GP and two matrixes TX and X that contain input vectors to GP. Returns covariance matrix C. Every element ij of C contains covariance between inputs i in TX and j in X. PREDCF is an optional array specifying the indexes of covariance functions, which are used for forming the matrix. If empty or not given, the matrix is formed with all functions.

Parameters:

x1	First input vector
x2	Second input vector
C	Covariance matrix

void extractUniqueIndices(const Eigen::VectorXd& input_ind, Eigen::VectorXd& output_ind) const

Extract unique indices from the input vector.

void evaluateEnergyAndGradient(const Eigen::VectorXd& w, const EigenMatrix& x, const Eigen::VectorXd& x_ind, const Eigen::VectorXd& y, EnergyAndGradient& energy_and_gradient)

Evaluate the energy function (un-normalized negative marginal log posterior) and its gradient.

This function takes a Gaussian process structure GP together with a matrix X of input vectors and a matrix Y of targets, and evaluates the energy function E and its gradient G. Each row of X corresponds to one input vector and each row of Y corresponds to one target vector.

The energy is minus log posterior cost function: E = EDATA + EPRIOR = - log p(Y|X, th) - log p(th), where th represents the parameters (lengthScale, magnSigma2…), X is inputs and Y is observations (regression) or latent values (non-Gaussian likelihood).

Parameters:

w	Vector of combined parameters from covariance functions and likelyhood

void optimize(const Observation& observation)

Optimize paramaters of a Gaussian process.

Main function to optimize the Gaussian Process hyperparameters.

This function optimizes the parameters of a GP structure given matrix X of training inputs and vector Y of training targets.

This function orchestrates the optimization process. It runs in a loop to ensure that the resulting hyperparameters are stable.

1. Selects a representative subset of data.

2. Performs numerical optimization of hyperparameters on this subset.

3. Monitors the new parameters for instability.

4. If instability is detected, the history size for the subset is increased, and the entire process is repeated.

5. If the parameters are stable, the loop terminates, and the final model is updated with the full dataset.

Parameters:

observation

The full set of observed data (coordinates, energies, gradients).

void optimizeFactorized(const Observation& observation)

Alternative optimization via eigendecomposition + noise grid search.

Assembles the noiseless covariance matrix, eigendecomposes it, then uses FactorizedLikelihood + NoiseGridSearch to find optimal noise. The analytical prefactor gives the signal variance (magnSigma2).

Only works for isotropic lengthscale (single lengthscale for all pairtypes). For ARD, falls back to optimize().

Parameters:

observation

The full set of observed data.

void calculatePotential(Observation& image1)

Calculate potential within GPR model.

This auxiliary function calculates energy and gradient predictions for given images R according to a GPR model.

Parameters:

image1	Observation at image 1
num_of_calls	Number of calls for GPR potential

void calculatePotentialLaplace(Observation& image1)

Laplace-mixture-mean potential prediction (INLA, Rue et al. 2009).

Equivalent to calculatePotential in API writes into image1.E and image1.G but the predictive mean is the weighted-mixture mean over the CCD grid centred at the current MAP hyperparameters. The posterior factor state is keyed by the training-state stamp; query coordinates do not enter that stamp.

Falls back to the standard calculatePotential if any precondition fails (Hessian non-PD, directional oracle refuses, training data absent). Always safe to call.

void calculatePotentialLaplaceSlice(Observation& image1)

Slice-sampling Laplace marginalisation (UncertaintyMode :: LAPLACE_SLICE). Caches M = slice_laplace_opts_.n_samples posterior samples (theta_k, alpha_k = K(theta_k)^-1 y) on the first call after a training-data or theta_* change, then averages per-sample mean predictions across the cache. No new Cholesky on hot calls; per-call cost is O(M * n_train * n_query). Falls back to MAP if the slice burst fails (signature-invalid, NaN initial log-target, kernel state inconsistent). See Murray & Adams 2010 NIPS, arxiv 1006.0868.

void calculateVarianceLaplaceSlice(Observation& image1)

Variance counterpart to calculatePotentialLaplaceSlice. Reports the between-sample variance of the per-sample mean predictions (Var_theta[mu(theta; x*)]); the within-sample variance is left for a follow-up since it requires K^-1 k_cross per sample per query (one extra triangular solve each, doubling per-predict cost). The between-sample variance is the dominant component when the hyperparameter posterior is broad (which is the only regime where LAPLACE differs meaningfully from MAP), so the approximation is fit-for-purpose for the variance-gated acquisition path.

void setUncertaintyMode(UncertaintyMode mode)

Select which predictive path consumers should use. Default MAP matches the pre-Laplace behaviour bit-for-bit; LAPLACE routes through calculatePotentialLaplace (CCD mixture), BMA is reserved for the reference path (not yet wired).

void setUsePsisLaplace(bool use)

Switch the LAPLACE_SLICE arm between the doubling slice sampler (default, false) and the PSIS-Laplace sampler (Vehtari 2024, arXiv:1507.02646 v9, vault note PSIS_Laplace_Design_2026_05.org). PSIS-Laplace replaces the slice’s M*~5 Cholesky burst with K Cholesky decomps + a Pareto k_hat reliability check.

void setUsePsisDeterministicGrid(bool use)

Switch the PSIS-Laplace sampler between random IS draws (default, false) and a deterministic axial CCD-style grid (true) of K = 1 + 2 p nodes placed at theta_MAP +/- sigma_j e_j on the principal axes of H^{-1}. Deterministic grid is np-reproducible; random IS retains MC variance per rank. Standard INLA Gaussian-mode quadrature (Rue 2009 sec 3.2, GPstuff INLA mode 0). Vault note INLA_Delta_Design_2026_05.org.

void setUseDeltaMethod(bool use)

Switch the LAPLACE_SLICE arm to the analytical INLA-Delta variance estimator (Rue 2009 sec 3.2; GPstuff INLA mode 0 ‘gaussian’). Replaces the per-call sampling/quadrature loop with the closed-form formula: Var[mu(x*) | y] = Var_MAP(x*) + g(x*)^T H^{-1} g(x*), g(x*)_j ~= k_cross^T dalpha/dtheta_j. SymPy proof D1-D7 in python/tests/test_sympy_inla_delta.py. Per-call cost O(p N n_q) (no per-call Cholesky); the slice/ PSIS sampling burst is skipped entirely. Vault note INLA_Delta_Design_2026_05.org.

void buildRFFModel(const Observation& training_data, int D_rff)

Random Fourier Features (RFF) layered predict path. Approximates the SE kernel with D_rff random cosine features, then solves Bayesian linear regression in the feature space. Per-call predict cost drops from O(N (D+1)) (cached MAP cross-cov + matvec) to O(D_rff D) (constant in N). Layered ON TOP of the LAPLACE_SLICE + INLA-Delta variance: the RFF mean replaces the cheap-mean MAP path; INLA-Delta still provides the calibrated variance for the gate. D_rff = 0 disables RFF (default; bit-identical to main).

void setRFFFeaturesRequested(int D_rff)

Configure the layered RFF predict path. D_rff > 0 enables; optimize() rebuilds the RFF model on the full observation after each successful SCG fit (matches the LaplacePosteriorStamp invalidation policy: theta_MAP changes -> rebuild). 0 disables (default; bit-identical to main).

void calculatePotentialDispatched(Observation& image1)

Dispatched predict: routes to calculatePotential (MAP) or calculatePotentialLaplace based on uncertainty_mode_. This is what AtomicDimer / GPNEB should call; the explicit non-dispatched methods stay available for direct use (e.g. tests).

void calculateVarianceLaplace(Observation& image1)

Laplace-mixture predictive variance (law of total variance over the CCD grid centred at theta*). Same shape contract as calculateVariance: writes per-test-point variance into image1.E (N_im scalars) and per-dim gradient variance into image1.G (N_im x D). Falls back to calculateVariance (MAP) if no Laplace fit can be built.

void calculateVarianceDispatched(Observation& image1)

Dispatched variance: calculateVariance (MAP) or calculateVarianceLaplace based on uncertainty_mode_. GPNEB and AcquisitionStrategy consumers call this to respect the CLI flag.

bool hasBeenOptimized() const

Has optimize() been called at least once in this lifetime? Some callers (tests, telemetry) gate on whether SCG has produced a defensible MAP estimate.

void setJitterSigma2(const Index_t value)

Set sigma2 jitter value.

void setParameters(const GPRSetup& parameters)

Set parameters of the GPR model.

const EigenMatrix& getL() const

Read-only accessors for training-state members. Used by downstream Laplace / TCOV-INLA helpers that previously had to friend class GaussianProcessRegression to reach the protected L / a / R_matrix / R_indices fields.

void assembleDirectionalCovarianceDerivatives(const EigenMatrix& x, const Eigen::VectorXd& ind_Ddim, const Eigen::VectorXd& uDdim, const Eigen::VectorXd& direction, EigenMatrix& first, EigenMatrix& second, bool serial_per_rank = false) const

Assemble the directional first and second derivatives of the training covariance matrix for the optimized Sexpat parameters.

Public + const so that out-of-class callers (e.g. laplace::TCovSecondOrder) can request directional kernel derivatives without needing friend access; promoted from private for that purpose.

void extractCoordinatesByIndex(const EigenMatrix& x, const Eigen::VectorXd& ind_Ddim, const Index_t ind, FieldMatrixd& x_loc) const

Extract coordinates from x using value ind from ind_Ddim.

Note

Exploits the contiguous block structure of ind_Ddim (blocks of N_im identical indices) to perform a single block copy from row ind*N_im instead of scanning all rows. x_loc must be pre-sized to the correct number of rows.

void assignBlockToMatrix(const Eigen::VectorXd& ind_Ddim1, const Eigen::VectorXd& ind_Ddim2, const Index_t row_val, const Index_t col_val, const FieldMatrixd& field, EigenMatrix& matrix, bool transpose_field = false) const

Assign field to matrix according to indices form ind_Ddim1 and ind_Ddim2.

This function performs the following operation (in MATLAB form):

\[matrix(ind_Ddim == row_val, ind_Ddim == col_val) = field;\]

Note

matrix should be pre-allocated

Parameters:

ind_Ddim1	Vector of dimensions used in rows
ind_Ddim2	Vector of dimensions used in columns
row_val	Index of row
col_val	Index of column
field	Reference field
matrix	Reference matrix
transpose_field	True if the field should be transposed

void accumulateBlockContribution(const Eigen::VectorXd& ind_Ddim1, const Eigen::VectorXd& ind_Ddim2, const Index_t row_val, const Index_t col_val, const FieldMatrixd& field, const EigenMatrix& invC, const Eigen::VectorXd& b, double& trace_accum, double& quad_accum, bool transpose_field = false)

Accumulate trace and quadratic-form contributions from one block.

Adds the contribution of the logical block matrix(ind_Ddim1 == row_val, ind_Ddim2 == col_val) = field to trace_accum += tr(invC * dK_block) quad_accum += b^T dK_block b without materializing the full dK matrix.

void removeColumn(const Index_t column_ind, EigenMatrix& matrix)

Remove column from the dense matrix.

Parameters:

column_ind	Column index
matrix	Dense matrix

Index_t getNumberOfRepetitiveIndices(const Eigen::VectorXd& ind_Ddim) const

Count number of the first repetitive indices in ind_Ddim.

virtual void calculateVariance(Observation& image1)

Compute predictive variance for energy and gradient. Corresponds to MATLAB get_variance.m.

virtual void refreshVariancePrefactor(const Eigen::VectorXd&)

Hook called inside calculateVarianceLaplace at each CCD grid point, right after calculateMeanPrediction (so b := L^{-1} y is valid) and before calculateVariance(scratch_var). Purpose: let subclasses refresh any hyperparameter-dependent cached scalar that their override of calculateVariance consumes.

Default: no-op. StudentTProcessRegression overrides this to recompute tp_prefactor_ = (2 hp_b + S) / (2 hp_a + n - 2) per grid point so the mixture variance reflects the TP’s full hyperparameter posterior (the MAP prefactor is stale at every off-centre CCD point).

bool evaluateDirectionalCurvature(const Eigen::VectorXd& x_ind, const EigenMatrix& x, const Eigen::VectorXd& y, const Eigen::VectorXd& direction, double& curvature)

Evaluate the exact directional curvature p^T H p for the current GP objective state.

Returns false when the current kernel parametrization cannot supply a backend-local exact directional derivative path, in which case callers should fall back to the finite-difference probe.

bool assembleExactHessian(const Eigen::VectorXd& x_ind, const EigenMatrix& x, const Eigen::VectorXd& y, int n_hyperparameters, Eigen::MatrixXd& H)

Assemble the full p x p observed-information Hessian of the negative log marginal likelihood at the current hyperparameters.

Uses the evaluateDirectionalCurvature oracle plus the polarisation identity validated symbolically in python/tests/test_sympy_laplace.py::test_polarisation_identity_recovers_off_diagonal:

H(i,i) = e_i^T H e_i (p direct calls) H(i,j) = ( (e_i+e_j)^T H (e_i+e_j)

• e_i^T H e_i

• e_j^T H e_j ) / 2 (p(p-1)/2 extra calls)

Total: p(p+1)/2 directional oracle calls, each of which reuses the cached Cholesky factor L and alpha. Symmetrised on the way out.

Parameters:

x_ind	unique-index vector (same semantics as evaluateDirectionalCurvature).
x	training input matrix.
y	training targets.
n_hyperparameters	p length of the hyperparameter vector the directional oracle is parameterised by (the direction argument to evaluateDirectionalCurvature has this length).
H	[out] p x p observed Hessian. Resized to (p, p) on successful return.

Returns:

true if the oracle returned a finite curvature at every probe direction. false on any failure; H is left NaN.

void updateModelWithFullData(const Observation& full_observation)

Updates the GPR model’s internal state using the full dataset.

After hyperparameters are optimized, this function re-evaluates the covariance matrix and its Cholesky decomposition using all available data points.

Parameters:

full_observation

The complete set of observations.

void updateConfiguration(const AtomsConfiguration& conf_info)

Update the internal atoms configuration. Use this when the set of active/frozen atoms changes dynamically.

void calculatePotentialSVGD(Observation& image1)

SVGD particle posterior predict + variance. Routed via calculatePotentialDispatched / calculateVarianceDispatched when uncertainty-mode=svgd.

void calculatePotentialBMA(Observation& image1)

Bayesian Model Averaging gold-standard reference predict and variance. Marginalises the predictive distribution over a 2D log10 grid (lengthScale x noise_var) weighted by the posterior p(theta_k | D) ~ exp(-NLML(theta_k)). Routed via calculatePotentialDispatched / calculateVarianceDispatched when uncertainty-mode=bma. The grid is fitted from scratch on each call (no cache); per-call cost is n_grid noiseless Cholesky decompositions, intentionally heavier than LAPLACE since BMA is the validation reference for the LAPLACE mixture.

void evaluateEnergyAndGradientLocal(const Eigen::VectorXd& w, const EigenMatrix& x, const Eigen::VectorXd& x_ind, const Eigen::VectorXd& y, EnergyAndGradient& energy_and_gradient)

Local-only counterpart to evaluateEnergyAndGradient: every linear algebra step uses Eigen LLT on the rank’s local heap rather than the ScaLAPACK collective Cholesky. Required by callers that want to evaluate the GP NLML at many hyperparameter samples in parallel across MPI ranks (each rank holds a different theta); a collective inside that loop would either deadlock or feed different matrices into the same global call.

Mutates rank-local kernel state via setParameters(w) the same way the distributed evaluator does; pair with ScopedKernelState at the call-site if the caller needs to restore (magn, ls, sigma2) after a sweep.

Does NOT touch member L, b, a, or C (the distributed Cholesky state); writes only into the caller-supplied energy_and_gradient out parameter.

void setSVGDConfig(const laplace::SteinVariationalConfig& cfg)

SVGD config getters/setters. Default-constructed SteinVariationalConfig matches the prior hardcoded literals (T=30, eps=0.05, conv_tol=1e-3) so existing callers stay bit-identical until they opt in.