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``. 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. 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. 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 :math:`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. 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); 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``. Trust, but verify ----------------- ``s.diagnostics.Q2`` is the leave-one-out cross-validated :math:`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 :math:`R^2`. 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). 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 :math:`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. 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.