The CMP Method

Complete Maximum a Posteriori method for Bayesian calibration of computer models under model discrepancy.

When computational models are calibrated against experimental observations, model discrepancy (or model-form error) is virtually inevitable due to physical simplifications, numerical discretization, or un-modeled phenomena. Standard Bayesian calibration following the classical Kennedy & O’Hagan (KOH) framework introduces a Gaussian Process (GP) discrepancy term \(z(\mathbf{x}) \sim \mathcal{GP}(0, k_{\boldsymbol{\psi}}(\cdot, \cdot))\) parameterized by hyperparameters \(\boldsymbol{\psi}\).

However, existing inference paradigms suffer from severe trade-offs:

  • Full Bayes (FB): Solves the joint problem over parameters \(\boldsymbol{\theta}\) and hyperparameters \(\boldsymbol{\psi}\). While statistically rigorous, it requires sampling high-dimensional, highly correlated joint spaces via MCMC, which is computationally prohibitive for expensive engineering codes and prone to identifiability issues.
  • Fixed/Sequential Modular (KOH): Estimates a single plug-in value \(\hat{\boldsymbol{\psi}}_{\text{KOH}}\) and conditions parameter inference on it. This dramatically underestimates parameter posterior uncertainty (the false certitude effect) and fails completely on multimodal posteriors by collapsing onto a single mode.
  • Fast Maximum a Posteriori (FMP): Allows \(\boldsymbol{\psi}\) to vary as a function of \(\boldsymbol{\theta}\), but ignores hyperparameter priors and the changing volume of the conditional distribution, leading to skewed posterior mode weights.

The Complete Maximum a Posteriori (CMP) method solves these limitations by establishing a modular, closed-form approximation of the marginal parameter posterior \(p(\boldsymbol{\theta} \mid \mathbf{y})\) that retains the full epistemic uncertainty of the hyperparameters without ever requiring high-dimensional joint sampling.

The CMP Formulation

The intuition begins with the exact decomposition of the joint posterior:

\[p(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \mathbf{y}) = p(\boldsymbol{\theta} \mid \mathbf{y}) \, p(\boldsymbol{\psi} \mid \boldsymbol{\theta}, \mathbf{y})\]

Rather than evaluating the joint distribution at a single fixed hyperparameter vector \(\hat{\boldsymbol{\psi}}\), the CMP method evaluates it along an optimal hyperparameter curve (or ridge) \(c(\boldsymbol{\theta}) = \left(\boldsymbol{\theta}, \boldsymbol{\psi}^*(\boldsymbol{\theta})\right)\), where:

\[\boldsymbol{\psi}^*(\boldsymbol{\theta}) = \arg\max_{\boldsymbol{\psi} \in \Psi} p(\boldsymbol{\psi} \mid \boldsymbol{\theta}, \mathbf{y}) = \arg\max_{\boldsymbol{\psi} \in \Psi} p(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \mathbf{y})\]

The Scale-Shift Assumption & Marginal Posterior

To integrate out the hyperparameters analytically, CMP assumes that the conditional distribution \(p(\boldsymbol{\psi} \mid \boldsymbol{\theta}, \mathbf{y})\) changes across parameter space \(\boldsymbol{\theta}\) primarily through affine transformations (shifting its mode and scaling its spread), governed by a fixed underlying shape function \(s(\boldsymbol{\phi})\):

\[\boldsymbol{\phi} = \mathbf{S}_{\boldsymbol{\theta}} (\boldsymbol{\psi} - \boldsymbol{\psi}^*(\boldsymbol{\theta})) \implies p(\boldsymbol{\psi} \mid \boldsymbol{\theta}, \mathbf{y}) \approx s(\boldsymbol{\phi}) \, \lvert\det \mathbf{S}_{\boldsymbol{\theta}}\rvert\]

Evaluating at \(\boldsymbol{\psi} = \boldsymbol{\psi}^*(\boldsymbol{\theta})\) (where \(\boldsymbol{\phi} = \mathbf{0}\)) yields \(p(\boldsymbol{\psi}^* \mid \boldsymbol{\theta}, \mathbf{y}) \, \lvert\det \mathbf{S}_{\boldsymbol{\theta}}\rvert^{-1} \approx s(\mathbf{0}) = \text{const}\). Multiplying by \(p(\boldsymbol{\theta} \mid \mathbf{y})\), the marginal parameter posterior is obtained in closed form:

\[p_{\text{CMP}}(\boldsymbol{\theta} \mid \mathbf{y}) \propto p(\mathbf{y} \mid \boldsymbol{\theta}, \boldsymbol{\psi}^*(\boldsymbol{\theta})) \, p(\boldsymbol{\theta}) \, p(\boldsymbol{\psi}^*(\boldsymbol{\theta})) \, \frac{1}{\lvert\det \mathbf{S}_{\boldsymbol{\theta}}\rvert}\]

This expression demonstrates that the marginal posterior can be accurately approximated along the ridge \(c(\boldsymbol{\theta})\), provided one divides by the volume correction factor \(\lvert\det \mathbf{S}_{\boldsymbol{\theta}}\rvert\) and accounts for the hyperparameter prior \(p(\boldsymbol{\psi}^*(\boldsymbol{\theta}))\).

Evaluation of the Hessian Correction Factor

The volume factor \(\lvert\det \mathbf{S}_{\boldsymbol{\theta}}\rvert\) reflects the local curvature (information volume) of the hyperparameter log-posterior at its peak:

\[\lvert\det \mathbf{S}_{\boldsymbol{\theta}}\rvert \propto \sqrt{\left\lvert\det \boldsymbol{\nabla}^2_{\boldsymbol{\psi}} \log p(\boldsymbol{\theta}, \boldsymbol{\psi}^* \mid \mathbf{y})\right\rvert}\]

For Gaussian likelihoods with covariance \(\mathbf{K}_{\boldsymbol{\psi}} = k_{\boldsymbol{\psi}}(\mathbf{X}, \mathbf{X}) + \boldsymbol{\Sigma}_e\) and residuals \(\mathbf{r} = \mathbf{y} - f(\mathbf{X}; \boldsymbol{\theta})\), this Hessian can be computed analytically:

\[\partial_{\psi_l} \partial_{\psi_k} \log p(\mathbf{y} \mid \boldsymbol{\theta}, \boldsymbol{\psi}) = \frac{1}{2} \text{Tr}\left( (\boldsymbol{\alpha}\boldsymbol{\alpha}^\top - \mathbf{K}_{\boldsymbol{\psi}}^{-1}) \frac{\partial^2 \mathbf{K}_{\boldsymbol{\psi}}}{\partial \psi_l \partial \psi_k} \right) - \text{Tr}\left( \text{Sym}(\mathbf{A}_l \boldsymbol{\alpha}\boldsymbol{\alpha}^\top) \frac{\partial \mathbf{K}_{\boldsymbol{\psi}}}{\partial \psi_k} \right) + \frac{1}{2} \text{Tr}\left( \mathbf{A}_l \mathbf{A}_k \right)\]

where \(\boldsymbol{\alpha} = \mathbf{K}_{\boldsymbol{\psi}}^{-1} \mathbf{r}\) and \(\mathbf{A}_i = \mathbf{K}_{\boldsymbol{\psi}}^{-1} \frac{\partial \mathbf{K}_{\boldsymbol{\psi}}}{\partial \psi_i}\). CMP is therefore a generalized, modular extension of the Laplace approximation applied along the parameter manifold.

Theoretical & Numerical Validation

1. Multimodal Gaussian Mixture

In a multimodal setting where the joint posterior is a mixture of Gaussians with well-separated modes:

  • KOH evaluates along a single average line (blue line), entirely missing the second mode.
  • FMP captures both modes but distorts their relative weights.
  • CMP (red ridge) analytically recovers the exact marginal posterior \(p(\boldsymbol{\theta} \mid \mathbf{y})\), preserving both modes and their precise mixture proportions.

2. Inadequate Model with Bimodal Posterior

When calibrating an inadequate model \(f(x; \theta) = (1-\theta)(x+0.15) + x \sin(2\theta x)\) against linear synthetic data \(y(x)=x\):

The true posterior has two distinct solutions (\(\theta \approx -0.1\) and \(\theta \approx 0.9\)). KOH captures only one mode with extreme false certainty. FMP overweights the minor mode. The CMP posterior matches the Full Bayes MCMC benchmark across the entire support, with the Hessian correction factor correctly penalizing low-curvature regions.

3. Sphere Drag Law (\(C_D - \text{Re}\))

Calibrating the empirical 3-parameter modified Stokes law \(C_D = \frac{A}{\text{Re}^B} + C\) against historical experimental data across \(\text{Re} \in [10^{-1}, 10^5]\):

Standard modular calibration (KOH) underestimates parameter variance by more than 50% due to the false certitude effect. CMP restores the full posterior variance (e.g., parameter \(A\) std: \(5.65\) vs Full Bayes \(6.87\), while KOH collapses to \(2.09\)), while keeping the dimension of the MCMC sampling space strictly restricted to 3 instead of 6.

Publications & Code

  • Journal Article: O. Kahol, P. M. Congedo, O. P. Le Maître, E. Denimal Goy, “Efficient treatment of the model error in the calibration of computer codes: the complete maximum a posteriori method”, International Journal for Uncertainty Quantification, 15(5): 57–76, 2025. DOI: 10.1615/int.j.uncertaintyquantification.2025056317.
  • C++ Implementation: Integrated within the open-source library CMP++.