Gaussian process posterior mean and variance of the loss after two noise-free runs

Writes the posterior mean and variance at an untried value in closed form when two runs have been made without noise. The weights a_1 and a_2 solve the two-by-two correlation system of the runs; the posterior mean moves the prior mean towards the observed losses by those weights, and the posterior variance is the prior variance less the part the runs explain. Correlations come from the kernel (HE-FM-BOPT-002); for more runs or noisy runs the matrix form is HE-CF-BOPT-001.

Signature

a_1 = (c_1 - c_12 * c_2) / (1 - c_12^2); a_2 = (c_2 - c_12 * c_1) / (1 - c_12^2); mu_n = mu_0 + a_1 * (y_1 - mu_0) + a_2 * (y_2 - mu_0); var_n = sf2 * (1 - (c_1 * a_1 + c_2 * a_2))
Inputs
InputsDefinitionUnit
c_1Kernel correlation between the candidate value and the first runnone
c_12Kernel correlation between the two evaluated parameter valuesnone
c_2Kernel correlation between the candidate value and the second runnone
mu_0Mean of the Gaussian process priorloss units
y_1Loss returned by the model at the first evaluated valueloss units
y_2Loss returned by the model at the second evaluated valueloss units
sf2Variance of the Gaussian process priorsquared loss units
Output
a_1Weight on the deviation of the first observed loss from the prior meannone
a_2Weight on the deviation of the second observed loss from the prior meannone
mu_nPredicted loss at the candidate value given the two runsloss units
var_nVariance of the loss at the candidate value given the two runssquared loss units

Function

Bayesian optimisation of a model calibration loss with a Gaussian process surrogate

Searches a bounded parameter region for the values that minimise a goodness-of-fit loss between model outputs and calibration targets when each model run is expensive. A Gaussian process fitted to the losses already observed gives a normal posterior for the loss at any untried value, and an acquisition function, here expected improvement, chooses the next run. The notation follows the Bayesian Optimisation article and its one-parameter incidence calibration.

Computational function

  • Computational function: Gaussian process posterior and expected improvement over candidate values from any number of runs

    Computes the Gaussian process posterior mean, standard deviation and expected improvement over a set of candidate values from any number of runs, with an optional noise variance, so that the next run can be chosen as the candidate with the highest expected improvement. The inputs differ from the formula's: vectors of evaluated values, losses and candidates instead of three correlations, and the normal distribution function is evaluated in the code.

    Inputs and outputs: X: Evaluated parameter values, a vector or a matrix with one row per run. Unit: parameter units.; y: Observed losses of the runs. Unit: loss units.; cand: Candidate values. Unit: parameter units.; mu0, sf2, ell: Prior mean, prior variance and length-scale. Unit: loss units, squared loss units, parameter units.; noise: Noise variance of a run, 0 for a deterministic model (a jitter of 1e-6 is added for numerical stability). Unit: squared loss units.; mean, sd, ei: Posterior mean, posterior standard deviation and expected improvement at each candidate, in input order. Unit: loss units.

    Assumption: Squared exponential kernel with fixed parameters; the best value in expected improvement is the lowest observed loss without noise and the lowest posterior mean at the evaluated values with noise.

    Worked example (Article example, two runs): With runs at 0.01 and 0.06 (losses 20.4382 and 16.1385), prior mean 20, variance 100 and length-scale 0.02, the candidates 0.020 to 0.050 have posterior means 20.01, 19.62, 19.11, 18.50, 17.84, 17.20 and 16.65, standard deviations 4.60, 6.30, 7.37, 7.74, 7.37, 6.30 and 4.60 and expected improvements 0.51, 1.15, 1.69, 2.05, 2.17, 2.02 and 1.59; the highest is at 0.040, as in the article. X = [0.01, 0.06]; y = [20.4382, 16.1385]; cand = [0.02, 0.025, 0.03, 0.035, 0.04, 0.045, 0.05]; ei = [0.51, 1.15, 1.69, 2.05, 2.17, 2.02, 1.59]

    Worked example (Second iteration, three runs): Adding the run at 0.04 (loss 3.5167), the expected improvement is about 0.4632 at 0.035 and 0.2876 at 0.030, so 0.035 is run next, as in the article. X = [0.01, 0.06, 0.04]; cand = [0.03, 0.035]; ei = [0.2876, 0.4632]

    Worked example (Noisy runs): With a noise variance of 4 on each run, the candidate 0.04 has posterior mean about 17.92, standard deviation 7.48 and expected improvement 2.24 (computed here for illustration). noise = 4; cand = [0.04]; mean = 17.92; sd = 7.48; ei = 2.24

    Excel: For one parameter, with runs in RunPoints and RunLosses, a candidate in Cand and PriorMean, SFVar, LengthScale and NoiseVar named, KMat is =SFVar*EXP(-(RunPoints-TRANSPOSE(RunPoints))^2/(2*LengthScale^2))+NoiseVar*MUNIT(ROWS(RunPoints)), KVec is =SFVar*EXP(-(Cand-RunPoints)^2/(2*LengthScale^2)), PostMean is =PriorMean+SUMPRODUCT(MMULT(MINVERSE(KMat),KVec),RunLosses-PriorMean), PostVar is =SFVar-SUMPRODUCT(KVec,MMULT(MINVERSE(KMat),KVec)), and expected improvement follows HE-FM-BOPT-004 with NORM.S.DIST.

    R: gp_ei <- function(X, y, cand, mu0, sf2, ell, noise = 0) { X <- as.matrix(X); C <- as.matrix(cand); k <- function(A, B) sf2*exp(-pmax(outer(rowSums(A^2), rowSums(B^2), "+")-2*A %*% t(B), 0)/(2*ell^2)); K <- k(X, X)+diag(noise+1e-6, nrow(X)); ks <- k(C, X); m <- as.vector(mu0+ks %*% solve(K, y-mu0)); s <- sqrt(pmax(sf2-rowSums(ks*t(solve(K, t(ks)))), 0)); best <- if (noise > 0) min(mu0+k(X, X) %*% solve(K, y-mu0)) else min(y); z <- (best-m)/pmax(s, 1e-12); data.frame(mean = m, sd = s, ei = ifelse(s > 0, (best-m)*pnorm(z)+s*dnorm(z), pmax(best-m, 0))) } Base R only; gp_ei(c(0.01, 0.06), c(20.4382, 16.1385), seq(0.02, 0.05, by = 0.005), 20, 100, 0.02) returns the first example.

    Python: def gp_ei(X, y, cand, mu0, sf2, ell, noise=0.0): y = np.asarray(y, float); X = np.asarray(X, float).reshape(len(y), -1); C = np.asarray(cand, float).reshape(len(cand), -1); k = lambda A, B: sf2*np.exp(-((A[:, None, :]-B[None, :, :])**2).sum(-1)/(2*ell**2)); K = k(X, X)+(noise+1e-6)*np.eye(len(y)); ks = k(C, X); m = [email protected](K, y-mu0); s = np.sqrt(np.maximum(sf2-(ks*np.linalg.solve(K, ks.T).T).sum(1), 0)); best = (mu0+k(X, X)@np.linalg.solve(K, y-mu0)).min() if noise > 0 else y.min(); z = (best-m)/np.maximum(s, 1e-12); Phi = np.array([0.5*math.erfc(-u/math.sqrt(2)) for u in z]); return {"mean": m, "sd": s, "ei": np.where(s > 0, (best-m)*Phi+s*np.exp(-z**2/2)/math.sqrt(2*math.pi), np.maximum(best-m, 0))} Needs import math and import numpy as np; returns the same values as the R function.

    Test (Posterior passes through a noise-free run): With Cand set to the first run and NoiseVar 0, PostMean equals the first loss and PostVar is close to 0. Expected result: TRUE. FALSE shows a noise term or the wrong length-scale in KMat. Excel check: =AND(ABS(PostMean-INDEX(RunLosses,1))<1E-3,PostVar<1E-3)

    Common error (Searching the acquisition on a coarse grid): The grid's best candidate is 0.040, but the continuous maximum of expected improvement is about 0.0398; a coarse grid can miss narrow peaks when several parameters are calibrated (computed here for illustration).

    Source: Rasmussen CE, Williams CKI. Gaussian Processes for Machine Learning. Cambridge, MA: MIT Press; 2006. Chapter 2, Regression. Section 2.2, equations 2.20 to 2.26 (predictive mean and variance with noise); Frazier PI. A tutorial on Bayesian optimization. arXiv preprint arXiv:1807.02811 [stat.ML], version 1, 8 July 2018. Sections 3, 4.1 and 5 (posterior, expected improvement in closed form, noisy evaluations).

    mu_n = mu_0 + k_n^T (K_n + noise I)^(-1) (y - mu_0); var_n = sf2 - k_n^T (K_n + noise I)^(-1) k_n; EI = (L_best - mu_n) Phi(z) + s_n phi(z)

Try this function

Implementations

  • Excel

    Two-run Gaussian process posterior from named cells

    With Corr1, Corr2, Corr12, PriorMean, Loss1, Loss2 and SFVar named, the formulas return the two weights, the posterior mean and the posterior variance, held in GPWeight1, GPWeight2, PostMean and PostVar.

    =(Corr1-Corr12*Corr2)/(1-Corr12^2); =(Corr2-Corr12*Corr1)/(1-Corr12^2); =PriorMean+GPWeight1*(Loss1-PriorMean)+GPWeight2*(Loss2-PriorMean); =SFVar*(1-(Corr1*GPWeight1+Corr2*GPWeight2))

Assumptions

  • Noise-free model runs

    The losses are observed without simulation noise, as for a deterministic cohort model; a microsimulation with a finite number of simulated people adds a noise variance to the diagonal of the run covariance matrix.

  • Fixed prior mean, variance and length-scale for the posterior

    The prior mean, variance and length-scale are held fixed when the posterior is computed.

Worked examples

  • Posterior at the candidate 0.04 after runs at 0.01 and 0.06

    With correlations of 0.3247 and 0.6065 to the runs and 0.0439 between them, the weights are about 0.2986 and 0.5934; the posterior mean is about 17.8394 and the variance about 54.3143, a standard deviation of about 7.37, as in the article (17.84 and 54.31).

    c_1 = 0.324652; c_2 = 0.606531; c_12 = 0.043937; mu_0 = 20; y_1 = 20.4382; y_2 = 16.1385; sf2 = 100; a_1 = 0.2986; a_2 = 0.5934; mu_n = 17.8394; var_n = 54.3143

  • Posterior at the candidate 0.035, midway between the runs

    Midway the two correlations are equal at about 0.4578 and each weight is about 0.4386; the posterior mean is about 18.4987 and the variance about 59.8421, a standard deviation of about 7.74, as in the article's table.

    c_1 = 0.4578334; c_2 = 0.4578334; c_12 = 0.043937; mu_0 = 20; y_1 = 20.4382; y_2 = 16.1385; sf2 = 100; a_1 = 0.4386; a_2 = 0.4386; mu_n = 18.4987; var_n = 59.8421

  • Posterior at an evaluated value, the run at 0.06

    At a run the correlation with it is 1, so the weights are 0 and 1; the posterior mean equals the observed loss of 16.1385 and the variance is 0 (computed here for illustration).

    c_1 = 0.043937; c_2 = 1; c_12 = 0.043937; mu_0 = 20; y_1 = 20.4382; y_2 = 16.1385; sf2 = 100; a_1 = 0; a_2 = 1; mu_n = 16.1385; var_n = 0

Common errors

  • Ignoring simulation noise in a microsimulation

    Treating noisy losses as exact makes the surrogate pass through the noise and report zero uncertainty at evaluated values; a noise variance on the diagonal (HE-AS-BOPT-005) avoids this.

  • Reading the surrogate as the disease model

    The posterior describes the loss, not the model's outputs, and can take impossible values: after another iteration in the example the posterior mean at 0.030 is about minus 1.09 (computed here for illustration).

Sources

  • Gaussian process posterior mean and variance with noise on the diagonal

    Rasmussen CE, Williams CKI. Gaussian Processes for Machine Learning. Cambridge, MA: MIT Press; 2006. Chapter 2, Regression. Section 2.2, equations 2.20 to 2.26: with additive independent Gaussian noise the covariance of the observations is K(X, X) plus sigma_n squared I; the predictive mean for one test point is k_* transpose (K plus sigma_n squared I) inverse y and the variance k(x_, x_) minus k_* transpose (K plus sigma_n squared I) inverse k_*.

    View source →

  • Posterior with a prior mean and variance reduced by the evaluations

    Frazier PI. A tutorial on Bayesian optimization. arXiv preprint arXiv:1807.02811 [stat.ML], version 1, 8 July 2018. Section 3, equation 3: the posterior mean is a weighted average of the prior mean and an estimate from the data, and the posterior variance equals the prior covariance less a term for the variance removed by observing the evaluations; section 5, noisy evaluations: a diagonal term equal to the noise variance is added, and the maximum of the posterior mean at evaluated points is typically used in place of the best observed value.

    View source →

Canonical Identity

Persistent URI
https://healtheconomics.wiki/formula/gaussian-process-posterior-mean-and-variance-of-the-loss-after-two-noise-free-runs

Stable URI · Machine-readable · Resolvable · CC BY 4.0