1.4. Surrogates: polynomial chaos, analytic Sobol’ indices

A surrogate (emulator) is a cheap statistical stand-in for an expensive function of the parameters – typically the log posterior kernel, whose every evaluation costs a model solve plus a filter pass. RISE fits the surrogate once, from a few hundred evaluations, and every subsequent evaluation costs microseconds. Two things come out of one fit:

a fast approximation of the expensive function over the plausible (prior) region of the parameter space, and
the complete variance decomposition of that function: because the polynomial basis is orthonormal under the input measure, Sobol’ sensitivity indices are read directly off the fitted coefficients – no additional sampling.

This page covers the two emulator families – polynomial chaos (PCE) and Gaussian process (Kriging) – and the emulate method. The machinery lives in rise.engine.surrogates (pce, gp, doe, and the abstract surrogate contract they share). A runnable walk-through is rise-modern-tutorials/GSAandUQ/UQ_surrogates/howto.m.

1.4.1. Quick start

With priors and data set on a model (nothing needs to have been estimated):

[s, design] = emulate(m);          % ~3 x nterms posterior evaluations

s.diagnostics                      % R2, Q2 (leave-one-out), N, P, ...

yhat = predict(s, theta);          % instant log-posterior approximation

sb = sobol(s);                     % analytic sensitivity indices

table(sb.names(:), sb.main, sb.total, ...
    'VariableNames', {'param','main','total'})

sb.main is each parameter’s own share of the variance of the log posterior over the prior; sb.total adds its interactions. Parameters with tiny total indices barely move the posterior – prime suspects for weak identification; parameters with large indices are where the data speak.

1.4.2. The prior chooses the polynomials

The defining trick: the polynomial family is derived per parameter from its prior distribution, so the basis is orthonormal under the prior measure (the Wiener-Askey correspondence):

a normal prior gets (probabilists’) Hermite polynomials;
any other RISE prior (beta, gamma, inverse gamma, …) is mapped through its own cdf onto a uniform variable, which gets Legendre polynomials – for a uniform prior this transform is linear, so the classical scheme is reproduced exactly;
parameters without a usable prior (e.g. positions inside a dirichlet block) fall back to Legendre on their bounds.

Two consequences. First, the design of experiments (rise.engine.surrogates.doe, a Latin-hypercube/Sobol/Halton design pushed through the priors’ inverse cdfs) places points where the prior puts mass, not in the corners of the bound box – important when priors are fat-tailed. Second, the Sobol’ decomposition is taken with respect to the prior, which is the measure a Bayesian cares about.

1.4.3. The Gaussian-process emulator: calibrated pointwise uncertainty

The second emulator on the same fit/predict contract:

s = emulate(m, 'method', "gp");

[yhat, ystd] = predict(s, theta);   % ystd varies point by point

Where the PCE’s uncertainty (its leave-one-out RMSE) is one global number, the GP’s ystd is pointwise and calibrated: small near the design points, growing away from them. That is the ingredient for adaptive designs and surrogate-guided optimization. The implementation is detrended (universal) Kriging: a polynomial mean fitted by least squares (trend option; emulate uses "quadratic", the natural shape of a log posterior) plus an anisotropic squared-exponential GP on the residuals, with inputs standardized by the design data and hyperparameters learned by marginal likelihood with restarts. On fs2000 the GP reaches a leave-one-out Q2 of ~0.92 from 135 design points (the order-2 PCE: ~0.81 from 165).

Trade-offs: fitting is \(O(N^3)\) per marginal-likelihood evaluation (designs of a few hundred points are comfortable, a few thousand are not), and there is no analytic Sobol’ decomposition – use the PCE for sensitivity analysis, the GP when pointwise uncertainty matters. Both work as delayed-acceptance screens.

1.4.4. Fitting on your own function

emulate accepts any scalar function of the parameter vector:

s = emulate(m, 'target', @(x) my_statistic_of(x), 'order', 3);

and the lower-level objects can be used standalone, without a model:

p = rise.engine.surrogates.pce(lb, ub, 'order', 2);

X = rise.engine.surrogates.doe(lb, ub, 200);

p = fit(p, X, y);

1.4.5. Options

order (default 2) – maximum total polynomial degree. Order 2 is the workhorse for log posteriors (locally quadratic); raise it if Q2 is poor.
qnorm (default 1) – hyperbolic truncation. Values below 1 drop high-order interaction terms, shrinking the basis in high dimension while keeping the univariate terms.
nsamples / oversampling (default 3) – design size; least-squares PCE wants 2-3 points per basis term.
scheme – latin_hypercube (default), sobol or halton.
use_priors – set false to force the uniform-on-bounds treatment.
use_parallel – evaluate the design under parfor.

1.4.6. Trust, but verify

s.diagnostics.Q2 is the leave-one-out cross-validated \(R^2\) – the honest in-design measure (the in-sample R2 always flatters). For an out-of-sample check, draw a fresh design from the same measure the surrogate was trained on:

Xv = rise.engine.surrogates.doe(lb, ub, 30, 'priors', design.priors);

% evaluate the truth on Xv, then:
holdout_r2(s, Xv, y_true)

Two warnings born of experience:

Validate on the prior, not on the bound box. A uniform design over [lb, ub] concentrates points in regions a fat-tailed prior essentially never visits; the surrogate was (rightly) never trained there and a variance-based score will collapse.
The deep tails are extrapolation. The log posterior can fall by thousands in a corner of the prior; no low-order polynomial tracks that, and it does not need to: the surrogate’s job is to be right where the posterior mass is. Rank (Spearman) correlation on validation points is a robust complement to \(R^2\).

1.4.7. Delayed-acceptance MCMC: exact posterior, fraction of the cost

The payoff of the surrogate layer. Set the da_surrogate property on a Metropolis-Hastings sampler (rwmh or imh) and every proposal is first screened with the surrogate; only survivors pay for a true posterior evaluation, and a second accept/reject step (Christen & Fox, 2005) makes the chain’s stationary distribution exactly the true posterior. The surrogate redistributes computation – it cannot bias the result. A poor surrogate costs acceptance rate, never correctness (this is unit-tested with a deliberately shifted surrogate):

[s, design] = emulate(m);                        % fit once

[ff, lb, ub] = pull_objective(m);

energy = @(x)-ff(x);                             % maximize the log posterior

opts = struct('N', 3000, 'burnin', 500, ...
    'tunedCov', C0, 'da_surrogate', @(x)predict(s, x));

smplr = rise.engine.estimation.rsamplers.rwmh(energy, x0, lb, ub, opts);

results = sample(smplr, struct('do_tuning', true));  % adaptive scale

results{1}.stats.funevals counts the true posterior evaluations actually paid for. On fs2000 the screen settles ~70% of the iterations without touching the true posterior, cutting wall time ~3x while the plain and delayed-acceptance posterior densities lie on top of each other – see rise-modern-tutorials/Estimation/surrogate_accelerated_mcmc/howto.m.

Practical notes:

the surrogate must approximate the same quantity and scale as the sampler’s target (the log posterior); @(x)predict(s, x) from emulate is already right;
proposal covariances built from a prior-wide design are far too large for a concentrated posterior – turn on the sampler’s adaptive tuning (sample(smplr, struct('do_tuning', true))) and let the scale find the standard 0.234 acceptance;
if the surrogate fails at the chain’s initial point, delayed acceptance is disabled for that chain (with a warning) rather than risking a chain that can never move;
the tempered apt and slice samplers do not support delayed acceptance (they error early rather than silently ignoring it).

1.4.8. Adaptive tools on top of the emulators

Three tools close the surrogate lifecycle (fit → refine → optimize → sample):

Refinement – refine(s, f, nadd) adds nadd true evaluations where the emulator is least sure (its pointwise sd for the GP; space-filling for the PCE) and refits, in rounds. Use it when the diagnostics say the fit is not yet trustworthy, or to sharpen the emulator around the region a chain is exploring:

s = refine(s, @(x)my_true_function(x), 30, 'priors', design.priors);

Surrogate-guided mode finding – rise.engine.surrogates.ei_maximize runs the classic EGO loop (Jones, Schonlau & Welch, 1998): fit a GP, evaluate the true function where the expected improvement is largest (exploit the surrogate’s optimum, explore where its uncertainty is big), repeat. Locating the mode of an expensive posterior in a few dozen evaluations:

energy = @(x)-ff(x);   % ff from pull_objective

[xmode, fmode, g] = rise.engine.surrogates.ei_maximize(energy, lb, ub, ...
    'budget', 80, 'priors', design.priors);

The returned GP is reusable – e.g. as a delayed-acceptance screen for the sampling stage that follows.

Active subspaces – rise.engine.surrogates.active_subspace(s) estimates the eigendecomposition of \(E[\nabla f \, \nabla f']\) under the prior, with the gradients taken on the fitted surrogate (microseconds; zero additional true evaluations). Large leading eigenvalues followed by a sharp drop mean the posterior effectively lives in a few linear parameter combinations – the activity field gives a quick “who matters” read that complements the Sobol’ indices.

1.4.9. One call, three views: the sensitivity front-end

sensitivity(m) runs the whole global-sensitivity toolbox on ONE evaluated design (the same one emulate uses – HDMR and Monte-Carlo filtering recycle it rather than sampling again):

res = sensitivity(m);

res.table
%  param    sobol_main  sobol_total  hdmr_first  mcf_ks  mcf_pval

sobol_main / sobol_total: who drives the variance of the log posterior over the prior (analytic, from the PCE);
hdmr_first: the same first-order question answered by an independent metamodel (RISE’s HDMR) fitted on the same design – a built-in cross-check. Only the first-order column is reported: the legacy HDMR engine’s higher-level aggregates show systematic leakage at realistic sample sizes (verified against analytic benchmarks), so read interactions off the PCE columns;
mcf_ks / mcf_pval: who decides whether a draw lands in the top region of the target (default: top 10%), measured by the Kolmogorov-Smirnov separation between the “behave” and “non-behave” conditional distributions – a question variance decompositions do not answer. The full mcf object is returned for its cdf and correlation-pattern plots.

res.surrogate and res.design carry the fitted PCE and the evaluated design for reuse – as a delayed-acceptance screen, in active_subspace, or for refinement.