Overcoming model simplifications when quantifying predictive uncertainty

TL;DR

This paper establishes conditions under which simplified models can reliably quantify predictive uncertainty, and proposes calibration strategies to overcome model simplifications for improved uncertainty estimation.

Contribution

It derives theoretical conditions for useful model simplifications and introduces calibration methods to mitigate simplification effects in uncertainty quantification.

Findings

Defined optimal and suboptimal model simplifications.

Proposed calibration strategies to compensate for simplifications.

Related conditions to common regularization schemes.

Abstract

It is generally accepted that all models are wrong -- the difficulty is determining which are useful. Here, a useful model is considered as one that is capable of combining data and expert knowledge, through an inversion or calibration process, to adequately characterize the uncertainty in predictions of interest. This paper derives conditions that specify which simplified models are useful and how they should be calibrated. To start, the notion of an optimal simplification is defined. This relates the model simplifications to the nature of the data and predictions, and determines when a standard probabilistic calibration scheme is capable of accurately characterizing uncertainty. Furthermore, two additional conditions are defined for suboptimal models that determine when the simplifications can be safely ignored. The first allows a suboptimally simplified model to be used in a way that replicates the performance of an optimal model. This is achieved through the judicial selection of a prior term for the calibration process that explicitly includes the nature of the data, predictions and modelling simplifications. The second considers the dependency structure between the predictions and the available data to gain insights into when the simplifications can be overcome by using the right calibration data. Furthermore, the derived conditions are related to the commonly used calibration schemes based on Tikhonov and subspace regularization. To allow concrete insights to be obtained, the analysis is performed under a linear expansion of the model equations and where the predictive uncertainty is characterized via second order moments only.

Figures

Equations

$$\mathbf{d}$$

$$\mathbf{p}$$

$$p(\mathbf{x}),p(\mathbf{d}|\mathbf{x}),p(\mathbf{p}|\mathbf{x}).$$

$$\hat{p}(\mathbf{v}),\hat{p}(\mathbf{d}|\mathbf{v}),\hat{p}(\mathbf{p}|\mathbf{v}).$$

$$p(\mathbf{x})=N(\mathbf{x};0,{\mathbf{\Sigma }}_{\mathbf{x}}).$$

$$p(\mathbf{d}|\mathbf{x})=N(\mathbf{d};\mathbf{Gx},{\mathbf{\Sigma }}_{\mathbf{\delta }}).$$

$$p(\mathbf{p}|\mathbf{x})=N(\mathbf{p};\mathbf{Yx},{\mathbf{\Sigma }}_{\mathbf{\rho }}).$$

$$p(\mathbf{x}|\mathbf{d})=\frac{p(\mathbf{x})p(\mathbf{d}|\mathbf{x})}{\int p(\mathbf{x})p(\mathbf{d}|\mathbf{x})d\mathbf{x}}.$$

$$\mathbf{\mu }{}_{\mathbf{x}|\mathbf{d}}$$

$${\mathbf{\Sigma }}_{\mathbf{x}|\mathbf{d}}$$

$$p(\mathbf{p}|\mathbf{d})=\int p(\mathbf{p}|\mathbf{x})p(\mathbf{x}|\mathbf{d})d\mathbf{x}.$$

$$\mathbf{\mu }{}_{\mathbf{p}|\mathbf{d}}$$

$$=\mathbf{Y}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }(\mathbf{G}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }+{\mathbf{\Sigma }}_{\mathbf{\delta }}{)}^{-1}\mathbf{d},$$

$${\mathbf{\Sigma }}_{\mathbf{p}|\mathbf{d}}$$

$$=\mathbf{Y}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{Y}}^{\top }+{\mathbf{\Sigma }}_{\mathbf{\rho }}-\mathbf{Y}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }(\mathbf{G}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }+{\mathbf{\Sigma }}_{\mathbf{\delta }}{)}^{-1}\mathbf{G}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{Y}}^{\top }.$$

$$\mathbf{\mu }{}_{\mathbf{p}|\mathbf{d}}$$

$$=\mathbf{Y}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }(\mathbf{G}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }+{\mathbf{\Sigma }}_{\mathbf{\delta }}{)}^{-1}\mathbf{d},$$

$${\mathbf{\Sigma }}_{\mathbf{p}|\mathbf{d}}$$

$$=\mathbf{Y}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{Y}}^{\top }+{\mathbf{\Sigma }}_{\mathbf{\rho }}-\mathbf{Y}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }(\mathbf{G}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{G}}^{\top }+{\mathbf{\Sigma }}_{\mathbf{\delta }}{)}^{-1}\mathbf{G}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{Y}}^{\top }.$$

$$\mathbf{d}$$

$$\mathbf{p}$$

$${\mathbf{x}}^{\prime }=\mathbf{C}{\mathbf{v}}^{\prime }.$$

$$\tilde{\mathbf{G}}$$

$$\tilde{\mathbf{Y}}$$

$$\mathbf{C}=\mathbf{US}{\mathbf{V}}^{\top }=\left[\begin{array}{cc}{\mathbf{U}}_{\mathbf{C}}& {\mathbf{U}}_{\bar{\mathbf{C}}}\end{array}\right]\left[\begin{array}{c}{\mathbf{S}}_{\mathbf{C}}\\ 0\end{array}\right]{\mathbf{V}}_{\mathbf{C}}^{\top }={\mathbf{U}}_{\mathbf{C}}{\mathbf{S}}_{\mathbf{C}}{\mathbf{V}}_{\mathbf{C}}^{\top },$$

$$\mathbf{x}=\mathbf{Cv}+{\mathbf{U}}_{\bar{\mathbf{C}}}\mathbf{u},$$

$$\mathbf{d}$$

$$\mathbf{p}$$

$$\mathbf{G}{\mathbf{U}}_{\bar{\mathbf{C}}}=0\text{and}\mathbf{Y}{\mathbf{U}}_{\bar{\mathbf{C}}}=0,$$

$$\mathbf{Z}=\left[\begin{array}{c}\mathbf{G}\\ \mathbf{Y}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{U}}_{\mathbf{Z}}& {\mathbf{U}}_{\bar{\mathbf{Z}}}\end{array}\right]\left[\begin{array}{cc}{\mathbf{S}}_{\mathbf{Z}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_{\mathbf{Z}}^{\top }\\ {\mathbf{V}}_{\bar{\mathbf{Z}}}^{\top }\end{array}\right].$$

$$\hat{p}{}^n(\mathbf{d}|\mathbf{v})$$

$$\hat{p}{}^n(\mathbf{p}|\mathbf{v})$$

$$p(\mathbf{v},\mathbf{u})=N\Bigl{(}\begin{bmatrix}\mathbf{v}\\ \mathbf{u}\end{bmatrix};\begin{bmatrix}0\\ 0\end{bmatrix},\begin{bmatrix}\mathbf{\Sigma}_{\mathbf{v}}&\mathbf{\Sigma}_{\mathbf{v}\mathbf{u}}\\ \mathbf{\Sigma}_{\mathbf{v}\mathbf{u}}^{\top}{}&\mathbf{\Sigma}_{\mathbf{u}}\end{bmatrix}\Bigr{)},$$

$${\mathbf{\Sigma }}_{\mathbf{v}}={\mathbf{C}}^{†}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{C}}^{†}{}^{\top },{\mathbf{\Sigma }}_{\mathbf{u}}={\mathbf{U}}_{\bar{\mathbf{C}}}^{\top }{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{U}}_{\bar{\mathbf{C}}},\text{and}{\mathbf{\Sigma }}_{\mathbf{vu}}={\mathbf{C}}^{†}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{U}}_{\bar{\mathbf{C}}}.$$

$$\hat{p}{}^n(\mathbf{v})=N(\mathbf{v};0,{\mathbf{\Sigma }}_{\mathbf{v}}).$$

$${\mathbf{\Sigma }}_{\mathbf{v}}={\mathbf{C}}^{†}{\mathbf{\Sigma }}_{\mathbf{x}}{\mathbf{C}}^{†}{}^{\top },$$

Full text

1 Introduction

This paper considers the problem of assessing uncertainty in a prediction made for a particular system based on the combination of specific measurement data and expert domain knowledge. The focus of this work is environmental systems, such as river basins and groundwater systems, however, the analysis is likely to be applicable to a much wider set of domains. Correctly addressing such environmental prediction problems is fundamental to ensuring these important system are managed appropriately and sustainabl. Probabilistic Bayesian methods provide a theoretically consistent set of rules to combine such site specific data and prior knowledge through the use of a system model. However, these methods often fail when naively applied to a model that does not capture the full complexity of the system.

A possible solution is the incorporation of more detail and structural variations in the system model such that more sources of uncertainty are included. However, it is important to admit that this modelling ad. infinitum is not a solution as unlimited resources are never available for detailed model construction and execution, and simplifications must be made, at least at some level.

However, deciding on what should be included in a model of a system and what may be ignored is not straightforward. For instance in groundwater hydrology there is no consensus on what is an appropriate level of parameterization detail, e.g within the hydraulic conductivity field [(24), (48), (47), (51)]. In addition to this are the related decisions of what processes should be explicitly represented and what can be ignored. Or alternatively, when is it reasonable to lump many different processes together under a semi-physical or non-physical black box model that directly represents input-output relationships [(17), (46)].

What is necessary is an understanding of the effects of simplifying assumptions made during model development, and appropriate ways of dealing with them when calibrating the model and generating predictions.

A general, qualitative, notion of model adequacy was explored by gupta_towards_2012 in terms of the issues faced within the surface water, groundwater, unsaturated zone, and terrestrial hydrometeorology modelling communities. The work considered intermediate stages within the modelling process from the initial perceptual understanding through to the construction of a computational model and introduced a pluralistic definition of model adequacy. On one extreme is the engineering viewpoint that defines a structurally adequate model as one that can reproduce the input-output relationship of the system, with well characterized uncertainties (error models). On the other extreme is the physical science viewpoint that requires an adequate model to be consistent with the underlying physical system gupta_towards_2012 .

This paper will consider systems that may be exposed to future disturbances and limited data is available. Within these problems a well characterized regression model (the engineering viewpoint) cannot be constructed and some physical insight of the system is required gupta_towards_2012 . Furthermore, focus is given to environmental management problems, where the role of a model is to help inform a subsequent decision problem related to risk management e.g. engineering design freeze_hydrogeological_1990 , groundwater management doherty_groundwater_2013 , or climate change rougier_uncertainty_2014 . This requires the specification of a probability distribution over the predictions of interest such that decision theoretic methods can be applied to quantify the risk and determine the optimal management strategy berger_statistical_1985 ; smith_bayesian_2010 .

1.1 Related Work

The issue of mismatch between reality and a numerical representation of a system has received considerable attention from may different perspectives. The presence of simplifications may be identified in the data as additional misfit that is not consistent with measurements errors alone. These additional errors are referred to as: model error mclaughlin_reassessment_1996 ; white_quantifying_2014 ; ljung_stochastic_2014 , model structural error beven_concept_2005 , structural noise doherty_short_2010 , model inadequacy kennedy_bayesian_2001 ; gupta_towards_2012 , model discrepancy goldstein_probabilistic_2004 ; strong_when_2014 , modelization uncertainties tarantola_inverse_2005 , and others.

Within the system identification and control community ljung_system_1999 ; ljung_perspectives_2010 the concept is often referred to as system under-modelling ninness_estimation_1995 ; reinelt_comparing_2002 and two main probabilistic approaches have been developed: stochastic embedding goodwin_stochastic_1989 ; goodwin_quantifying_1992 ; ljung_stochastic_2014 and model error modelling ljung_model_1999 . These allow subjective information about the discrepancy between the transfer functions that describe the dynamics of the real system and that of a numerical model to be explicitly defined and combined with a time series dataset collected from the system. This in turn allows prediction uncertainty to include uncertainty due to modelling errors and measurement errors contained within the data. These approaches generally allow a very flexible black box representation of the system dynamics to be used.

In the statistical modelling community, computational simulators have been considered as an explicit approximation to the real physical system they are attempting to represent. The mismatch is represented as an additional error term and modelled probabilistically, generally with explicit space time correlation kennedy_bayesian_2001 ; craig_bayesian_2001 ; goldstein_probabilistic_2004 ; goldstein_bayes_2006 ; rougier_probabilistic_2007 ; rougier_uncertainty_2014 . Due to the computational complexity of the methods for high dimensional models, approximate Bayesian methods have also been developed goldstein_bayes_2006 . Such simulators may also be decomposed and internal error terms included to represent the structural errors within smaller submodels strong_when_2014 .

Other work within hydrology has examined model simplifications and how these effect the calibration and prediction processes. In particular it has been shown that model parameters that are designed to represent particular physical properties may be forced to undertake surrogate roles to compensate for the simplifications during the model calibration process. This surrogacy has the potential to introduce additional biases into the predictions that are unaccounted for by typical probabilistic methods. This has been analyzed by explicitly modelling the simplification process in mclaughlin_distributed_1988 ; doherty_short_2010 ; doherty_use_2011 ; watson_parameter_2013 ; white_quantifying_2014 .

Furthermore, several strategies have been considered to overcome the approximations inherent in a simplified model. This includes the generalized likelihood uncertainty estimation (GLUE) method beven_equifinality_2001 that modifies the data likelihood function such that less information is extracted from the data. For a dynamical system, the assumption of a deterministic dynamic transition function can be relaxed and a stochastic process or transition function used instead vrugt_improved_2005 ; clark_unraveling_2006 (such models are often referred to as data assimilation methods). Similarly the assumption that the model parameters are time invariant may be related such that they can change over time to better match the observed system behavior reichert_analyzing_2009 .

It is noted that many of these methods rely on a noticeable discrepancy between the measured data and model predicted values. However, the adverse effects of model simplifications may not necessarily cause any additional misfit between the data and model output, and thus may go unnoticed during model calibration and prediction white_quantifying_2014 . Thus it is critical to have a theoretical understanding of model adequacy that goes beyond data fit.

1.2 Approach and Contributions

This paper considers a subjective Bayesian framework and formally defines what is an adequate or appropriate model by explicitly considering: (i) the nature of the predictions, (ii) the available data and prior knowledge, (iii) the model simplification strategy, and (iv) the calibration scheme. It is shown that these are intrinsically linked and performing model calibration should take into account the nature of the data, predictions and any inherent limitations of the simplified computational model. It is noted that this departs from classical guidance on calibration and inversion that focuses on data misfit only mclaughlin_reassessment_1996 ; tarantola_inverse_2005 ; aster_parameter_2012 .

A dual model approach is used, where a “reference” mclaughlin_distributed_1988 or “reality” white_quantifying_2014 model is used to describe how the system is believed to function. Such a high fidelity reference model has been used in the past to characterize the performance of a simplified model. Here, the approach is extended to explicitly determine when and how a simplified model can be used to generate an accurate, or at least conservative, estimate of uncertainty in a prediction of interest. To allow concrete insights to be produced, only linear(ized) problems are considered and the beliefs are restricted to be Gaussian, where the second order moment is sufficient to capture uncertainty.

Within this framework, model simplifications are represented as a subspace projection that restricts the flexibility of the simplified model. This prevents the model from being able to fully represent the complexity of how the real system is believed to function. It is shown that if the model simplifications are ignored and standard probabilistic calibration methods are used, the model may generate overconfident and non-conservative predictions. Generally, this is avoided only when the simplification strategy is optimal. The key contribution of this paper is the characterization of two new calibration and prediction schemes for suboptimally simplified models that avoid this under estimation of uncertainty. The first scheme allows a simplified model to be used in a way that replicates the performance of an optimal scheme through the appropriate modification of the prior, or regularization, term. In the second scheme, focus is given on linking the predictions with the available data such that insight is gained into when the simplifications can be overcome by gathering the appropriate set of data.

Specifically, this paper fully defines the general problem of assessing uncertainty with a simplified model in Section 2. Section 3 summarizes the optimal benchmark solution based on the high fidelity reference model. The general representation of model simplifications based on linear projection is defined in Section 4. A typical probabilistic calibration scheme is considered in Section 5 that ignores the effects of model simplifications within the calibration and prediction process. This section defines the concept of an optimal simplification that determines when such a scheme is adequate for assessing uncertainty. Section 6 considers suboptimally simplified models and introduces two new calibration and prediction schemes, and defines the conditions when they are adequate for assessing predictive uncertainty. Discussions and directions for future research are covered in Section 7. A worked example involving a simplified groundwater prediction problem is given in Appendix A.

2 Problem Definition

Consider a system where a modeller is tasked with generating a probability distribution over a prediction of a specific feature of the system to inform a future decision process. For instance in a groundwater system, the prediction of interest may be the reduction in water levels that will be experienced by an aquifer due to an increase in future extractions. The prediction of interest will be denoted by the vector $\mathbf{p}\in\mathbb{R}^{D_{\mathbf{p}}}$.

To aid the modeller, a limited number of measurements have been made on the system and will be represented by $\mathbf{d}\in\mathbb{R}^{D_{\mathbf{d}}}$. Furthermore, it is considered that expert prior knowledge exists on how the system functions, for instance the physical processes involved, and allows the available data and prediction of interest to be related to the properties of the underlying system.

This information shall be considered to define a high fidelity or reference model for the system with parameters denoted by the vector $\mathbf{x}\in\mathbb{R}^{D_{\mathbf{x}}}$. Furthermore, the relationship between the data $\mathbf{d}$, predictions $\mathbf{p}$, and the underlying parameters $\mathbf{x}$ that the modeller believes exists shall be represented by the pair of equations

[TABLE]

It is considered that these functions include the operation of known physical laws, such as conservation of mass, energy and momentum, while the vector $\mathbf{x}$ includes a detailed representation of the forcing terms, initial conditions, material properties, etc. These requirements result in a very high dimensional vector $\mathbf{x}$ and complex functions $\mathcal{G}$ and $\mathcal{Y}$.

Furthermore, it is considered that epistemic uncertainty exists as to the value of $\mathbf{x}$ for the system under investigation. In addition, the errors on the right of Eq. 1 and Eq. 2 allow for the inclusion of additional uncertainty, for instance to capture measurement errors introduced by the data acquisition method, or uncertainty that may exist in how a physical process actually functions tarantola_inverse_2005 .

This allows the modeller’s complete belief about how the system behaves to be captured by three probability distribution functions

[TABLE]

These in turn can be used to generate a posterior distribution over the predictions of interest, denoted as $p(\mathbf{p}|\mathbf{d})$, using standard rules of probability theory.

It is noted that a major issue with pursuing this type of subjective Bayesian approach is the requirement that these probability distributions must represent every detail that is believed to be present in the physical system. This is not possible and simplification must be made.

2.1 Simplifications

It is considered that the process of modelling, for instance as outlined in gupta_towards_2012 , transforms the detailed knowledge of the modeller and produces the desired simplified computational representation of the problem that can be solved within the typical computational constraints. Specifically, the output of the modelling process is (or should be) an approximate representation of the modeller’s beliefs that can be processed numerically to generate a prediction posterior distribution that is somehow similar to the optimal, but computationally intractable, distribution $p(\mathbf{p}|\mathbf{d})$ generated by the high fidelity model. The modelling process is depicted in Fig. 1 where the outputs are denoted by the set of approximate beliefs

[TABLE]

These probability distributions are defined over a simplified description of the system, denoted by the parameter vector $\mathbf{v}\in\mathbb{R}^{D_{\mathbf{v}}}$. Furthermore, it is considered that the approximate posterior $\hat{p}(\mathbf{p}|\mathbf{d})$ should not just be similar to posterior $p(\mathbf{p}|\mathbf{d})$, it should be in some sense conservative such that the uncertainty is not underestimated doherty_groundwater_2013 .

Thus, the key question that is addressed in this paper is now specified as:

How can the subjective beliefs of a modeller, defined for the high fidelity reference model, be transformed to generate approximate probabilities that define a simplified computational model such that the computed prediction posterior $\hat{p}(\mathbf{p}|\mathbf{d})$ is a conservative approximation of the optimal reference posterior $p(\mathbf{p}|\mathbf{d})$? *

This question explicitly links the typically separate problems of: model construction, calibration and prediction.

Furthermore, answers will be provided under the restriction that the uncertainties are described by multivariate Gaussian distributions and the dependencies are linear. This will allow more specific and concrete insights to be gained into the effects of simplifications and how they may be overcome. Extensions to more general non-Gaussian and nonlinear systems is left for future work.

3 Optimal Inference and Prediction

The optimal Bayesian prediction scheme, which does not consider any simplifications, is now briefly reviewed berger_statistical_1985 ; tarantola_inverse_2005 ; aster_parameter_2012 . This is based on the high fidelity reference model.

Let the predictions of interest, the available data, and the parameters of the high fidelity model be defined by the random vectors: $\mathbf{p}\in\mathbb{R}^{D_{\mathbf{p}}}$, $\mathbf{d}\in\mathbb{R}^{D_{\mathbf{d}}}$, and $\mathbf{x}\in\mathbb{R}^{D_{\mathbf{x}}}$, where $D_{\mathbf{p}}$, $D_{\mathbf{d}}$ and $D_{\mathbf{x}}$ are their respective dimensions. It is considered that the dimension of $\mathbf{x}$ is excessively large, while the data is limited, and the number of predictions is small (perhaps only one), such that $D_{\mathbf{x}}\gg D_{\mathbf{d}},D_{\mathbf{p}}$.

Prior knowledge of the parameter vector is represented by the Gaussian probability density function $p(\mathbf{x})$, with mean of zero and known second order moment $\mathbf{\Sigma}_{\mathbf{x}}$ such that

[TABLE]

It is noted that the requirement of a zero mean distribution simplifies the analysis and can be met in general with a transformation into increments.

Now, the believed relationship between the system properties $\mathbf{x}$, and the data vector $\mathbf{d}$ is considered to be linear with coefficient matrix $\mathbf{G}$. This can be considered as a linearized approximation of Eq. 1, however the approximation introduced by the linearization is considered out of scope in this paper. In addition, the uncertainty that is believed to exist in this relationship is represented by the error term $\boldsymbol{\delta}$ such that $\mathbf{d}=\mathbf{G}\mathbf{x}+\boldsymbol{\delta}$. Furthermore, the uncertainty in the value of $\boldsymbol{\delta}$ is represented by a zero mean Gaussian density with known covariance matrix of $\mathbf{\Sigma}_{\boldsymbol{\delta}}$. In the simplest case, the matrix $\mathbf{\Sigma}_{\boldsymbol{\delta}}$ represents errors within the measurement process. This allows the modeller’s conditional belief of the measured data given the underlying parameters to be defined by the Gaussian density

[TABLE]

The prediction of interest $\mathbf{p}$ is also considered to be linearly dependent on $\mathbf{x}$ with a coefficient matrix $\mathbf{Y}$, and an additive error denoted by $\boldsymbol{\rho}$, such that $\mathbf{p}=\mathbf{Y}\mathbf{x}+\boldsymbol{\rho}$. In addition, it is considered that the uncertainty in the value of the error is a zero mean Gaussian density with known covariance $\mathbf{\Sigma}_{\boldsymbol{\rho}}$. It is noted that this requires that $\boldsymbol{\delta}$ and $\boldsymbol{\rho}$ are independent. If any correlation exists, it must be included in $\mathbf{x}$. Finally, the conditional belief of the prediction given the system properties can now be defined as

[TABLE]

Note that if the system is believed to be well described by a deterministic model, the covariance of the prediction error could be considered negligible.

The calibration or inversion stage seeks to determine the posterior belief over the parameters $\mathbf{x}$ given the available data $\mathbf{d}$

[TABLE]

Here, the posterior density $p(\mathbf{x}|\mathbf{d})=N(\mathbf{x};\boldsymbol{\mu}{}_{\mathbf{x}|\mathbf{d}},\mathbf{\Sigma}_{\mathbf{x}|\mathbf{d}})$ is Gaussian with a mean and covariance given by

[TABLE]

The coefficient matrix of the data $\mathbf{d}$ in Eq. 8 is sometimes referred to as the optimal estimator, or gain matrix, and will be denoted by $\mathbf{E}=\mathbf{\Sigma}_{\mathbf{x}}\mathbf{G}^{\top}{}(\mathbf{G}\mathbf{\Sigma}_{\mathbf{x}}\mathbf{G}^{\top}{}+\mathbf{\Sigma}_{\boldsymbol{\delta}})^{-1}$.

To propagate the posterior belief into the prediction of interest, requires the integration over the underlying parameters

[TABLE]

The prediction posterior density $p(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\boldsymbol{\mu}{}_{\mathbf{p}|\mathbf{d}},\mathbf{\Sigma}_{\mathbf{p}|\mathbf{d}})$ is Gaussian with mean and covariance given by

[TABLE]

The overall calibration and prediction scheme is now summarized below. This sets the benchmark for the other schemes considered in this paper.

Definition 1** (Optimal Scheme).**

The optimal scheme is denoted by the set of linear Gaussian probability density functions $p(\mathbf{x})$, $p(\mathbf{d}|\mathbf{x})$ and $p(\mathbf{p}|\mathbf{x})$, parameterized by the matrices $\mathbf{\Sigma}_{\mathbf{x}}$, $\mathbf{G}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\mathbf{Y}$, and $\mathbf{\Sigma}_{\boldsymbol{\rho}}$, as defined in Eqs. 5 to 7.

For a given dataset $\mathbf{d}$, the posterior belief over the predictions of interest is the Gaussian density $p(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\boldsymbol{\mu}{}_{\mathbf{p}|\mathbf{d}},\mathbf{\Sigma}_{\mathbf{p}|\mathbf{d}})$, with mean and covariance defined by the functions ${\mu}(\cdot)$ and $\mathit{\Sigma}(\cdot)$ respectively

[TABLE]

*Note that the covariance matrix of the posterior is not a function of the data $\mathbf{d}$. *

The structure of the above prediction scheme that encapsulates the prior densities $p(\mathbf{x})$, $p(\mathbf{d}|\mathbf{x})$, $p(\mathbf{p}|\mathbf{x})$ that are based on the high fidelity model is depicted graphically in Fig. 2.

4 Model Simplifications

The simplified computational model is considered to contain an approximate system simulator, with exposed parameters $\mathbf{v}$ and linear structure analogous to that defined for the high fidelity model

[TABLE]

Here $\tilde{\mathbf{G}}$ and $\tilde{\mathbf{Y}}$ are linear data and prediction matrices of a simplified system simulator. The exact form of the approximations within these relationships are of critical interest and will be explored by explicitly representing the simplifications involved.

To achieve this consider initially that the parameters $\mathbf{v}$ are somehow physically inspired and have some meaning to the modeller in describing aspects of the system. This allows the parameters of the simplified model to be used to describe a restricted version of the high fidelity model with parameters $\mathbf{x}^{\prime}\in\mathbb{R}^{D_{\mathbf{x}}}$ determined by the parameter vector of the simplified model $\mathbf{v}^{\prime}\in\mathbb{R}^{D_{\mathbf{v}}}$ via a matrix $\mathbf{C}\in\mathbb{R}^{D_{\mathbf{x}}\times D_{\mathbf{v}}}$

[TABLE]

This explicit representation of model simplifications allows aspects of the system to be ignored doherty_short_2010 ; white_quantifying_2014 and allows for the incorporation of arbitrary linear parameterization schemes, such as spatial and temporal homogeneity mclaughlin_reassessment_1996 .

With this definition, the simplified data and prediction matrices may be rewritten in terms of the matrices of the high fidelity model

[TABLE]

The matrix $\mathbf{C}$ will be referred to as a simplification matrix and provides a way of linking the simplified simulator to the original high fidelity reference model of how the system functions.

Finally, it is important to note that physical intuition is not strictly required when defining the simple model as for any $\mathbf{G}$, $\mathbf{Y}$ and $\tilde{\mathbf{G}}$, $\tilde{\mathbf{Y}}$, such that \bigl{[}\begin{smallmatrix}\mathbf{G}\\ \mathbf{Y}\end{smallmatrix}\bigr{]} has full row rank, there is always a $\mathbf{C}$ that satisfies Eq. 15 and Eq. 16. However, this may not be unique and thus there are potentially many different ways to interpret the meaning of the simple model. Here, it is considered that the modeller specifies the interpretation by e.g. specifying $\mathbf{C}$ directly.

4.1 Representing Unmodelled Complexity

A simplification matrix $\mathbf{C}\in\mathbb{R}^{D_{\mathbf{x}}\times D_{\mathbf{v}}}$, divides the space $\mathbb{R}^{D_{\mathbf{x}}}$ into two perpendicular subspaces: $\operatorname{columnspace}(\mathbf{C})$ and $\operatorname{cokernel}(\mathbf{C})$. The column space contains the set of parameter vectors that can be explicitly represented by a low dimensional vector. The second subspace $\operatorname{cokernel}(\mathbf{C})$ contains vectors that includes some degree of complexity that cannot be represented by a low dimensional vector.

Before defining these further, the following assumption is made.

Assumption 1**.**

*A simplification matrix $\mathbf{C}$ has full column rank, i.e. $\operatorname{rank}(\mathbf{C})=D_{\mathbf{v}}$. *

It is noted that if a simplification matrix $\mathbf{C}$ does not have full column rank, then there will be parameter combinations of the simple model that have the same influence on the high fidelity model, which implies that the simple model is not as simple as it could be for the same expressive power.

Now, consider the singular value decomposition of $\mathbf{C}$ under the above assumption

[TABLE]

where $\mathbf{U}_{\mathbf{C}}$ is the collection of orthogonal unit vectors that span the column space of $\mathbf{C}$ and similarly $\mathbf{U}_{\bar{\mathbf{C}}}$ spans the subspace perpendicular to this, the cokernel of $\mathbf{C}$.

This enables the parameter vector of the high fidelity model to be expanded as

[TABLE]

where $\mathbf{v}$ is the parameter vector of the simple model and $\mathbf{u}$ is a random vector that captures all the unmodelled complexity. Finally, it is noted that this expansion is unique and only dependent on the definition of the simplification matrix $\mathbf{C}$.

Under the above expansion, the data and prediction is related to the two components $\mathbf{v}$ and $\mathbf{u}$ through

[TABLE]

From these equations, it is noted that the unmodelled components of the system will be expressed in the data whenever $\mathbf{G}\mathbf{U}_{\bar{\mathbf{C}}}\neq 0$ and in the predictions whenever $\mathbf{Y}\mathbf{U}_{\bar{\mathbf{C}}}\neq 0$. These additional error terms are denoted as $\boldsymbol{\eta}$ and $\boldsymbol{\epsilon}$ respectively and play a very different role to the other error terms $\boldsymbol{\delta}$ and $\boldsymbol{\rho}$ as they are correlated via $\mathbf{u}$. It is noted that the additional error $\boldsymbol{\eta}$ introduced by the model simplifications was explicitly considered in mclaughlin_reassessment_1996 and tarantola_inverse_2005 and was used to define a composite measurement error model for the sum $\boldsymbol{\delta}+{\boldsymbol{\eta}}$. However, the effects of the simplifications on the prediction, including the explicit correlation between $\boldsymbol{\eta}$ and $\boldsymbol{\epsilon}$, was not addressed.

In an ideal world it could be argued that the objective of any modelling exercise should be to construct a model that is a simplified version of how the system is believed to function that captures the full complexity of the available data and predictions of interest. That is, which has a $\mathbf{C}$ such that $\boldsymbol{\eta}=\boldsymbol{\epsilon}=0$. Such a simplification matrix always exists, and will be referred to as an optimal simplification.

Definition 2** (Optimal Simplification).**

Let $\mathbf{G}$ and $\mathbf{Y}$ be data and predictions matrices for a high fidelity system model. Then a simplification matrix $\mathbf{C}$ is optimal if

[TABLE]

*where $\mathbf{U}_{\bar{\mathbf{C}}}$ spans the cokernel of $\mathbf{C}$, as defined in Eq. 17. *

It is noted that the above definition of optimal simplification differs from that of doherty_short_2010 ; doherty_use_2011 ; watson_parameter_2013 in that it explicitly considers the nature of the predictions.

An optimal simplification matrix, $\mathbf{C}^{*}$, can be found for any given data and prediction matrices $\mathbf{G}$ and $\mathbf{Y}$, by first taking a singular value decomposition of the matrix \mathbf{Z}=\bigl{[}\begin{smallmatrix}\mathbf{G}\\ \mathbf{Y}\end{smallmatrix}\bigr{]}

[TABLE]

Now, an optimal simplification matrix is given by $\mathbf{C}^{*}=\mathbf{V}_{\mathbf{Z}}$.

It is considered that an optimally simplified model is difficult, or even achievable, to obtain in practice and the remainder of the paper will focus on understanding and overcoming the issues caused by suboptimal simplifications.

5 Naive Use of Simplified Models

Here, a prediction scheme is defined that first infers the parameters of the simple model and then propagates them into the prediction, but ignores the potential for unmodelled complexity to be expressed in the data or predictions. This represents a standard probabilistic Bayesian approach to model calibration and prediction. It is shown that this scheme is generally not conservative as the uncertainty in the predictions is underestimated.

5.1 Prior Information

The scheme will have the same basic structure as the optimal benchmark scheme. In particular it will include an explicit representation of a prior over the parameters $\hat{p}{}^{n}{}(\mathbf{v})$, a data likelihood $\hat{p}{}^{n}{}(\mathbf{d}|\mathbf{v})$ and conditional density for the prediction $\hat{p}{}^{n}{}(\mathbf{p}|\mathbf{v})$. Here the superscript $n$ denotes approximate distributions associated with this naive method.

The conditional densities $\hat{p}{}^{n}{}(\mathbf{d}|\mathbf{v})$ and $\hat{p}{}^{n}{}(\mathbf{p}|\mathbf{v})$ are defined directly with the simplified data and prediction matrices $\tilde{\mathbf{G}}$ and $\tilde{\mathbf{Y}}$ respectively. Furthermore, the covariances are set to those used in the optimal scheme i.e.

[TABLE]

This formulation of the data and prediction densities is equivalent to ignoring the approximations in Eq. 12 and Eq. 13 introduced by the model simplifications.

Now, to be somewhat rigorous, the specification of a prior distribution on the parameters $\mathbf{v}$ is performed by explicitly considering the nature of the simplification and the uncertainty in the parameters of the high fidelity model. This is accomplished by considering a transformation that propagates the prior distribution over the parameters of the high fidelity model $\mathbf{x}$ into the joint space formed by the parameter vector of the simple model $\mathbf{v}$ and the vector $\mathbf{u}$ that denotes the unmodelled complexity. This is defined in the following proposition.

Proposition 1** (Uncertainty Propagation).**

Let $p(\mathbf{x})=N(\mathbf{x};0;\mathbf{\Sigma}_{\mathbf{x}})$ be a zero mean Gaussian density, with covariance $\mathbf{\Sigma}_{\mathbf{x}}$, that captures the prior uncertainty in the parameters of a high fidelity model. Also, let $\mathbf{C}$ be a simplification matrix that describes a simplified model.

Then, the joint probability density function for $\mathbf{v}$ and $\mathbf{u}$, representing the modelled and unmodelled components of the simple model respectively, is the zero mean correlated Gaussian density

[TABLE]

where the covariance terms are given by

[TABLE]

*Here, $(\cdot)^{\dagger}{}$ is the pseudoinverse operator and $\mathbf{U}_{\bar{\mathbf{C}}}$ spans the cokernel of $\mathbf{C}$ as defined in Eq. 17. *

See Appendix B for a proof.

This proposition defines the covariance matrix $\mathbf{\Sigma}_{\mathbf{v}}$, that represents the uncertainty in the parameters of the simplified model, to be the orthogonal projection of $\mathbf{\Sigma}_{\mathbf{x}}$ onto the column space of $\mathbf{C}$. This is the most rigorous approach to specifying the prior over the parameters of the simple model as it preserves the uncertainty that is believed to exist and is consistent with the physical intuition used to construct the simplified model.

This marginal density is used to define the prior of the naive prediction scheme

[TABLE]

5.2 Naive Calibration and Prediction

The naive calibration and prediction scheme is now defined that generates a posterior via the standard rules of probability theory. This is similar to the optimal scheme, but without consideration to the simplifications within the model’s data and prediction equations.

Definition 3** (Naive Scheme).**

A naive prediction scheme is denoted by the simplified linear Gaussian density functions $\hat{p}{}^{n}{}(\mathbf{v})$, $\hat{p}{}^{n}{}(\mathbf{d}|\mathbf{v})$ and $\hat{p}{}^{n}{}(\mathbf{p}|\mathbf{v})$ parameterized by the matrices $\mathbf{\Sigma}_{\mathbf{v}}$, $\tilde{\mathbf{G}}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\tilde{\mathbf{Y}}$, and $\mathbf{\Sigma}_{\boldsymbol{\rho}}$.

Furthermore, the prior covariance matrix $\mathbf{\Sigma}_{\mathbf{v}}$ is defined in terms of a high fidelity model using 1, that is

[TABLE]

where $\mathbf{\Sigma}_{\mathbf{x}}$ denotes a covariance matrix that describes the uncertainty in the parameters of the high fidelity model, and $\mathbf{C}$ denotes the simplification matrix that links the parameters of the simple and high fidelity models.

For a given dataset $\mathbf{d}$, the posterior belief over the predictions of interest generated by this naive scheme is the Gaussian density $\hat{p}{}^{n}{}(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{n},\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{n})$, with mean and covariance defined as

[TABLE]

*where the functions ${\mu}(\cdot)$ and $\mathit{\Sigma}(\cdot)$ are given in 1. *

5.3 Performance

To assess the performance of the naive scheme the notion of conservativeness, discussed in Section 2, is first defined by considering the squared error in the mean prediction. This is performed for a pair of probability density functions in 4, and extended to calibration and prediction schemes in 5. A further generalization of these definitions, which explicitly incorporates the utility function of the subsequent decision problem, is proposed in Section 7.

Definition 4** (Conservative Density).**

Let $p(\mathbf{p})$ and $\hat{p}(\mathbf{p})$ be density functions defined over the random vector $\mathbf{p}$. Furthermore, let $\hat{\boldsymbol{\mu}}{}_{\mathbf{p}}$ be the mean of $\hat{p}(\mathbf{p})$. Then $\hat{p}(\mathbf{p})$ is defined as a conservative approximation of the reference density $p(\mathbf{p})$ if the approximate density’s expected mean squared error $E_{\hat{p}(\mathbf{p})}\{(\hat{\boldsymbol{\mu}}{}_{\mathbf{p}}-\mathbf{p})(\hat{\boldsymbol{\mu}}{}_{\mathbf{p}}-\mathbf{p})^{\top}{}\}$, does not increase when the expectation is instead taken with respect to reference density $p(\mathbf{p})$

[TABLE]

*Where $\mathbf{A}\succeq\mathbf{B}$ requires that $\mathbf{A}-\mathbf{B}$ is positive semi-definite. *

It is noted that condition Eq. 28 can be rewritten in terms of the means $\boldsymbol{\mu}{}_{\mathbf{p}}$, $\hat{\boldsymbol{\mu}}{}_{\mathbf{p}}$ and covariances $\mathbf{\Sigma}_{\mathbf{p}}$, $\hat{\mathbf{\Sigma}}_{\mathbf{p}}$ of the two densities, and yields the condition

[TABLE]

An example of a conservative density, that satisfies 4, is given in Fig. 3.

The notion of conservativeness is now generalized to a prediction scheme. This is performed by considering an expectation over typical data.

Definition 5** (Conservative Scheme).**

Let the densities $p(\mathbf{x})$, $p(\mathbf{d}|\mathbf{x})$, and $p(\mathbf{p}|\mathbf{x})$ denote the prior information of a reference scheme. Also, let $p(\mathbf{p}|\mathbf{d})$ denote the posterior density this scheme generates for a given dataset $\mathbf{d}$. Furthermore, let $\hat{p}(\mathbf{p}|\mathbf{d})$ denote a posterior density generated by an approximate scheme.

For the given dataset $\mathbf{d}$, the degree of conservativeness of the approximate posterior is denoted by the function

[TABLE]

Now, the approximate scheme is defined as conservative if the expectation of the function $\Omega(\mathbf{d})$ is positive semi-definite

[TABLE]

*Here, the density $p(\mathbf{d})$ denotes how probable different datasets are under the prior knowledge of the reference scheme and is given by $p(\mathbf{d})=\int p(\mathbf{x})p(\mathbf{d}|\mathbf{x})d\mathbf{x}$. *

It is noted that for the linear Gaussian densities considered in this work, the covariance of the posterior distributions are independent of the data and Eq. 30 can be simplified as

[TABLE]

Thus, a scheme is considered as conservative if it generates a posterior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}$ that is inflated with respect to $\mathbf{\Sigma}_{\mathbf{p}|\mathbf{d}}$ by an amount dependent on the average squared difference in the posterior means.

It is now of interest to consider the performance of the naive scheme and determine when it is conservative with respect to the optimal benchmark scheme. This is performed in the following proposition.

Proposition 2** (Performance of Naive Scheme).**

Suppose the matrices $\mathbf{\Sigma}_{\mathbf{x}}$, $\mathbf{G}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\mathbf{Y}$, $\mathbf{\Sigma}_{\boldsymbol{\rho}}$ define an optimal calibration and prediction scheme, with posterior density denoted by $p(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\boldsymbol{\mu}{}_{\mathbf{p}|\mathbf{d}},\mathbf{\Sigma}_{\mathbf{p}|\mathbf{d}})$ for arbitrary dataset $\mathbf{d}$.

Let $\mathbf{C}$ denote a simplification matrix such that the matrices $\mathbf{\Sigma}_{\mathbf{v}}=\mathbf{C}^{\dagger}{}\mathbf{\Sigma}_{\mathbf{x}}\mathbf{C}^{\dagger}{}{}^{\top}{}$, $\tilde{\mathbf{G}}=\mathbf{G}\mathbf{C}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\tilde{\mathbf{Y}}=\mathbf{Y}\mathbf{C}$, $\mathbf{\Sigma}_{\boldsymbol{\rho}}$ define a naive prediction scheme for a simplified model. Now, consider the posterior density generated by the naive scheme $\hat{p}{}^{n}{}(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{n},\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{n})$.

If $\mathbf{C}$ is an optimal simplification, then the generated mean and covariance of the naive prediction scheme are equivalent to those of the optimal prediction scheme, that is

[TABLE] 2. 2.

If $\mathbf{C}$ is a suboptimal simplification, then the naive scheme is conservative, provided the following condition holds

[TABLE]

where $\mathbf{E}{}^{n}{}=\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}(\tilde{\mathbf{G}}\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}+\mathbf{\Sigma}_{\boldsymbol{\delta}})^{-1}$ denotes the estimator matrix of the naive scheme and $\mathbf{U}_{\bar{\mathbf{C}}}$ spans the cokernel of $\mathbf{C}$ as defined in Eq. 17. 3. 3.

If $\mathbf{C}$ is a suboptimal simplification, condition Eq. 32 does not hold, $\mathbf{u}$ and $\mathbf{v}$ are considered independent, and the covariance in the unmodelled complexity is non-zero $\mathbf{\Sigma}_{\mathbf{u}}\neq 0$, then the naive scheme is strictly non-conservative.

See Appendix C for a proof.

Under a suboptimal simplification, the non-conservative nature of the prediction scheme is caused by the influence of unmodelled complexity on the data and/or predictions. If it influences the predictions, then there is a direct bias effect. If it influences the data, then this causes the estimated parameters of the simple model to become biased by forcing them to take on surrogate roles to compensate for the simplifications. Furthermore, this bias in the estimated parameters is propagated to the prediction. Neither of these effects are taken into account by the naive prediction scheme and the uncertainty is underestimated.

The condition in Eq. 32, which guarantees the scheme is conservative, requires that these two sources of errors in the prediction caused by the unmodelled complexity to exactly cancel each other out. This is unlikely to hold in practical scenarios without careful attention to the data, predictions, simplifications, and the type of prior knowledge available. This condition will form the basis of the data driven prediction scheme introduced in Section 6.

Finally, it is noted that for the general case when $\mathbf{v}$ and $\mathbf{u}$ are not independent and Eq. 32 does not hold, the naive scheme is still not guaranteed to be conservative. However there are more special cases which may produce conservative posterior densities. The conditions that delineate the strictly conservative from the non-conservative scenarios are likely to be of limited interest in realistic problems and have not been enumerate here.

5.4 Summary

This section has formally defined a probabilistic calibration and prediction scheme that embeds the standard separation of modelling, calibration and prediction. In particular, physical knowledge is used to define a prior probability distribution for the parameters included in the simplified model. Also, the conditional data and prediction probability distributions are defined by ignoring the simplifying assumptions used to construct the data and prediction equations. Furthermore, it is shown that a true characterization of predictive uncertainty typically requires an optimally simplified model. For suboptimally simplified models it is shown that this scheme generally underestimates the uncertainty in the predictions and is not conservative.

The dependency structure of the probability distributions for this naive scheme is depicted in Fig. 4. Also displayed is the dependency on the unmodelled complexity for an optimally simplified model. It is noted that under this condition the data and predictions are conditionally independent of the unmodelled complexity $\mathbf{u}$ given the model parameters $\mathbf{v}$.

Finally, it is noted that the non-conservative, or overconfident, nature of predictions generated by naively applying probabilistic Bayesian methods on simplified models is not a new finding. In particular, the result can be considered as a generalization of doherty_short_2010 ; doherty_use_2011 ; white_quantifying_2014 and is consistent with arguments of beven_equifinality_2001 .

6 Overcoming Model Simplifications

It is of interest now to understand what to do about the issue of overconfidence when a suboptimal model is used. Firstly, it is noted that this issue may not be a problem. For instance, the prediction uncertainty may turn out to be larger than expected, and this in itself may provide sufficient information to allow a given management decision to be made. However, the situation becomes more difficult when it is necessary to determine an accurate or at least conservative prediction probability distribution, for instance as required by doherty_groundwater_2013 . In this scenario the modeller may opt to:

• •

Expand the sources of uncertainty considered in the model through less simplistic modelling assumptions with the hope of producing an optimally simplified model. This is not just about representing greater spatial/temporal detail in the model hunt_are_2007 , but more importantly the goal should be the inclusion of more sources of uncertainty. As an example, if there exists uncertainty over the presence of some important structural feature, a multi-hypothesis, or multi-model, analysis may be needed poeter_all_2007 ; poeter_mma_2007 ; rojas_application_2010 .

• •

Directly model the additional errors introduced by the model simplifications. For instance by defining an additional probabilistic model for $p(\boldsymbol{\eta},\boldsymbol{\epsilon})$ introduced in Eq. 19 and Eq. 20 that captures the uncertainty the modeller has in how well the simplistic model represents the behavior of the system. Such approaches have been developed in kennedy_bayesian_2001 ; craig_bayesian_2001 ; goldstein_probabilistic_2004 ; goldstein_bayes_2006 ; rougier_probabilistic_2007 ; rougier_uncertainty_2014 .

• •

Within the context of a dynamical system, remove the deterministic assumption on how the system evolves over time. This allows the simplifications to be represented by an error model within the modelled transition function. Such a scheme is employed by vrugt_improved_2005 . Alternatively, the assumption that the model parameters are time invariant can be relaxed such that the parameters can change over time to better match the underlying system behavior and inject greater uncertainty in the predictions reichert_analyzing_2009 .

• •

Modify the way that the data is used through changing the likelihood model $p(\mathbf{d}|\mathbf{v})$ such that less information is propagated into the prediction via the model parameters. The generalized likelihood uncertainty estimation method beven_equifinality_2001 can be considered as an example of this approach.

In the remainder of this section, two additional approaches are defined, along with explicit conditions that determine when they are appropriate. Firstly a scheme is defined that allows the posterior of the ideal prediction scheme to be reproduced using a suboptimally simplified model through the appropriate modifications of the prior covariances. Secondly, structural considerations of the data and predictions are considered and an approach developed that allows the simplifications to be overcome through the use of the right calibration data.

6.1 Optimal Use of Simplified Models

The naive scheme developed previously was not conservative for a suboptimally simplified model. Here a prediction scheme is defined that is capable of reproducing the optimal scheme, even under suboptimal simplifications.

To start, consider the probability densities $\hat{p}{}^{o}{}(\mathbf{v})$, $\hat{p}{}^{o}{}(\mathbf{d}|\mathbf{v})$, $\hat{p}{}^{o}{}(\mathbf{p}|\mathbf{v})$ to be defined in a similar manner to the naive scheme, but parameterized in terms of new covariance matrices $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$, $\hat{\mathbf{\Sigma}}_{\boldsymbol{\delta}}^{o}$ and $\hat{\mathbf{\Sigma}}_{\boldsymbol{\rho}}^{o}$ such that

[TABLE]

Here the covariance matrices will be considered as adjustable such that they can be chosen in such a way that they compensate for the simplifications and allow the optimal posterior to be reproduced. Thus, it is of interest to select $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$, $\hat{\mathbf{\Sigma}}_{\boldsymbol{\delta}}^{o}$ and $\hat{\mathbf{\Sigma}}_{\boldsymbol{\rho}}^{o}$ such that the posterior density that is produced by their combination, $\hat{p}{}^{o}{}(\mathbf{p}|\mathbf{d})$, replicates the posterior of the optimal scheme, i.e.

[TABLE]

The posterior density $\hat{p}{}^{o}{}(\mathbf{p}|\mathbf{d})$ is Gaussian, with mean and covariance defined in a similar fashion to the optimal scheme. Thus, the above condition will be satisfied when the means and covariances are equivalent, which occurs when

[TABLE]

If these conditions can be met, it defines a set of optimal covariance matrices of the simplified probability density functions such that when they are used to generate a posterior, the performance of the optimal scheme is replicated.

6.1.1 Highly Parameterized Models

The above conditions will now be specialized for a highly parameterized model, where the calibration problem is under-constrained. Such models are recommended by hunt_are_2007 and white_quantifying_2014 , and the results will explicit determine when such approaches are appropriate. This will be performed by focusing attention only on the prior covariance matrix for the parameters of the simple model, such that $\hat{\mathbf{\Sigma}}_{\boldsymbol{\delta}}^{o}$ and $\hat{\mathbf{\Sigma}}_{\boldsymbol{\rho}}^{o}$ remain the same as those used in the naive and optimal schemes, i.e. $\hat{\mathbf{\Sigma}}_{\boldsymbol{\delta}}^{o}=\mathbf{\Sigma}_{\boldsymbol{\delta}}$ and $\hat{\mathbf{\Sigma}}_{\boldsymbol{\rho}}^{o}=\mathbf{\Sigma}_{\boldsymbol{\rho}}$. Such a scheme is now defined.

Definition 6** (Optimally Compensated Scheme).**

An optimally compensated scheme is denoted by the simplified linear Gaussian probability density functions $\hat{p}{}^{o}{}(\mathbf{v})$, $\hat{p}{}^{o}{}(\mathbf{d}|\mathbf{v})$ and $\hat{p}{}^{o}{}(\mathbf{p}|\mathbf{v})$ parameterized by the matrices $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$, $\tilde{\mathbf{G}}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\tilde{\mathbf{Y}}$, and $\mathbf{\Sigma}_{\boldsymbol{\rho}}$.

To fully define the prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$, let $\mathbf{C}$ denote a simplification matrix that links the simplified model to a high fidelity model with data and prediction matrices denoted by $\mathbf{G}$ and $\mathbf{Y}$ such that $\tilde{\mathbf{Y}}=\mathbf{Y}\mathbf{C}$ and $\tilde{\mathbf{G}}=\mathbf{G}\mathbf{C}$. Furthermore, let $\mathbf{\Sigma}_{\mathbf{x}}$ denote the covariance matrix that describes the uncertainty in the parameters of this high fidelity model. The prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ is now defined as

[TABLE]

where $\mathbf{R}=\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}$, and \mathbf{Z}=\bigl{[}\begin{smallmatrix}\mathbf{G}\\ \mathbf{Y}\end{smallmatrix}\bigr{]}, \tilde{\mathbf{Z}}=\bigl{[}\begin{smallmatrix}\tilde{\mathbf{G}}\\ \tilde{\mathbf{Y}}\end{smallmatrix}\bigr{]}.

For a given dataset $\mathbf{d}$, the posterior belief over the predictions of interest generated by this scheme is the Gaussian density $\hat{p}{}^{o}{}(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{o},\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{o})$, with mean and covariance defined as

[TABLE]

*where the functions ${\mu}(\cdot)$ and $\mathit{\Sigma}(\cdot)$ are given in 1. *

This scheme is now shown to be equivalent to the optimal scheme when applied to a suboptimal but highly parameterized model where the number of free parameters is at least as large as the number of linearly independent measurements and predictions of interest.

Proposition 3** (Performance of Compensated Scheme).**

Suppose the matrices $\mathbf{\Sigma}_{\mathbf{x}}$, $\mathbf{G}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\mathbf{Y}$, $\mathbf{\Sigma}_{\boldsymbol{\rho}}$ define an optimal calibration and prediction scheme, with posterior density denoted by $p(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\boldsymbol{\mu}{}_{\mathbf{p}|\mathbf{d}},\mathbf{\Sigma}_{\mathbf{p}|\mathbf{d}})$ for arbitrary dataset $\mathbf{d}$.

Let $\mathbf{C}$ denote a suboptimal simplification matrix and let the matrices $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}=\mathbf{R}\mathbf{\Sigma}_{\mathbf{x}}\mathbf{R}^{\top}{}$, $\tilde{\mathbf{G}}=\mathbf{G}\mathbf{C}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\tilde{\mathbf{Y}}=\mathbf{Y}\mathbf{C}$, and $\mathbf{\Sigma}_{\boldsymbol{\rho}}$ define an optimally compensated scheme, where $\mathbf{R}=\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}$ and \mathbf{Z}=\bigl{[}\begin{smallmatrix}\mathbf{G}\\ \mathbf{Y}\end{smallmatrix}\bigr{]}, \tilde{\mathbf{Z}}=\mathbf{Z}\mathbf{C}=\bigl{[}\begin{smallmatrix}\tilde{\mathbf{G}}\\ \tilde{\mathbf{Y}}\end{smallmatrix}\bigr{]}. Now, consider the posterior density generated by the optimally compensated scheme $\hat{p}{}^{o}{}(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{o},\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{o})$. If the simplification is chosen such that

[TABLE]

Then, the posterior mean and covariance generated by the optimally compensated scheme are equivalent to those of the optimal scheme

[TABLE]

See Appendix D for a proof.

This demonstrates that given a model simplified in a suboptimal fashion, the optimal performance can be recovered through the adjustment of the assumed prior uncertainties. The main rank condition in Eq. 38 is fairly easy to achieve in practice and is met when, e.g. the number of parameters is not smaller than the number of linearly independent measurements and predictions, and the simplified model matrix $\tilde{\mathbf{Z}}=\mathbf{Z}\mathbf{C}$ has full (row) rank.

A model may be highly parameterized due to the degrees of spatial variability in e.g. the modelled material properties that is included in the model (as in white_quantifying_2014 ). However, the above result can apply to models of dynamical systems that use stochastic transition models (sometimes referred to as data assimilation methods). In these models the additional dynamic error terms that are incorporated at each time increment can be similarly viewed as a set of additional model parameters.

The general structure of the optimally compensated calibration and prediction scheme is depicted with a Bayesian network in Fig. 5. It is noted that the difference between this scheme and the naive scheme lies only in how the prior covariance matrix for the parameters of the simple model is specified. They are both defined as a transformation of the covariance matrix $\mathbf{\Sigma}_{\mathbf{x}}$ that captures the uncertainty in the parameters for the high fidelity model. Recall that the naive scheme uses

[TABLE]

while the optimally compensated scheme uses

[TABLE]

The matrix $\mathbf{R}$ is dependent on the simplification matrix $\mathbf{C}$ but it is also dependent on the data and predictions matrices $\mathbf{G}$ and $\mathbf{Y}$ of the high fidelity models. It is this explicit dependency on the data and predictions that allows the prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ to compensate for the errors in the data and prediction equations of the simplified model.

To gain additional intuition as to how the compensating prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ is related to $\mathbf{\Sigma}_{\mathbf{v}}$ used by the naive scheme, consider for the moment that $\mathbf{v}$ and $\mathbf{u}$ are independent. This enables the prior covariance matrix for the parameters of the high fidelity model to be decomposed as $\mathbf{\Sigma}_{\mathbf{x}}=\mathbf{C}\mathbf{\Sigma}_{\mathbf{v}}\mathbf{C}^{\top}{}+\mathbf{U}_{\bar{\mathbf{C}}}\mathbf{\Sigma}_{\mathbf{u}}\mathbf{U}_{\bar{\mathbf{C}}}^{\top}{}$. With this, the prior matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ can be rewritten as

[TABLE]

where $\mathbf{\Sigma}_{+}=\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}\mathbf{U}_{\bar{\mathbf{C}}}\mathbf{\Sigma}_{\mathbf{u}}\mathbf{U}_{\bar{\mathbf{C}}}^{\top}{}\mathbf{Z}^{\top}{}\tilde{\mathbf{Z}}^{\dagger}{}{}^{\top}{}$ and is positive semi-definite matrix. This means that for $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ to compensate for the simplifications, the parameters must be allowed greater flexibility than would be given by the naive use of $\mathbf{\Sigma}_{\mathbf{v}}$ as $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}\succeq\mathbf{\Sigma}_{\mathbf{v}}$.

6.1.2 Weighted Least Squares Inversion

In addition to the Bayesian arguments provided above, now consider a regularized weighted least squares formulation commonly used for calibration and inversion problems mclaughlin_reassessment_1996 ; tarantola_inverse_2005 . With the assumption of linear models, the mean vector $\hat{\boldsymbol{\mu}}_{\mathbf{v}|\mathbf{d}}^{o}$ is also the maximum a posteriori estimate of the posterior parameter distribution of the simplified model $\hat{p}{}^{o}{}(\mathbf{v}|\mathbf{d})\propto\hat{p}{}^{o}{}(\mathbf{d}|\mathbf{v})\hat{p}{}^{o}{}(\mathbf{v})$, and may be equivalently defined as the solution to the regularized weighted least squares optimization problem

[TABLE]

Note that this explicitly employs the simplified forward model $\tilde{\mathbf{G}}\mathbf{v}$. Furthermore, the prediction generated by this estimate $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{o}=\tilde{\mathbf{Y}}\hat{\boldsymbol{\mu}}_{\mathbf{v}|\mathbf{d}}^{o}$ is the mean and mode of the prediction posterior $\hat{p}{}^{o}{}(\mathbf{p}|\mathbf{d})$. Due to the optimality of the compensated scheme, this point prediction is also the optimal minimum error variance prediction that the high fidelity model would generate.

This demonstrates that a simplified, but highly parameterized model, can be calibrated through the careful selection of a Tikhonov regularizer to generate the optimal prediction. Additionally, this optimal regularizer can also be used to characterize the predictive uncertainty of this point prediction through the use of the covariance formula in Eq. 37.

6.1.3 Summary

The introduced optimally compensated prediction scheme allows suboptimal model simplifications to be overcome through modifying the prior probability distribution for the parameters of the simplified model. This explicitly links: (i) the modelling problem, through the selection of a simplification matrix $\mathbf{C}$; (ii) the nature of the available data via $\mathbf{G}$; (iii) the predictions of interest defined by $\mathbf{Y}$; and (iv) how the simple model is calibrated, as the required regularization term $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ is explicitly dependent on all three terms.

It is noted however, that this explicit dependency on the high fidelity model will make the generation of $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ problematic in practice as the matrix $\mathbf{G}$ and $\mathbf{Y}$ will not generally be available. Nevertheless, it formally defines how a prior regularization term for a suboptimal model should be selected.

6.2 Data Driven Predictions

It was noted that the main barrier to applying the optimally compensated scheme defined above is the difficulty in determining the covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ that compensates for the simplifications. Here a new scheme is proposed that overcomes this issue by moving away from optimality and focusing on conservativeness.

In particular, a specific class of predictions will be considered that allow the simplifications to be hidden behind the data. Thus, the objective is to examine how model simplifications can be overcome through the collection of the right data. The intuition here is that data of a similar type to the predictions is often available, for example groundwater head measurements are often available when it is of interest to predict heads, similarly stream flow measurements are often available when predictions of stream flows are of interest, etc.

The class of predictions that will be considered here incorporates two important principles:

• •

The effects of all unmodelled components of the system on the predictions are captured by the data.

• •

Any component of the system that effects the predictions, and is not captured by the data, is explicitly represented in the simplified model.

From these principles it is considered that the prediction of interest can be explicitly decomposed into an intermediate data term that is dependent on $\mathbf{g}=\mathbf{G}\mathbf{x}$ and a term wholly dependent on the parameters of the simple model, that is

[TABLE]

where $\bar{\mathbf{A}}$ and $\bar{\mathbf{B}}$ are arbitrary matrices. Predictions with this structure are conditionally independent of $\mathbf{u}$ given $\mathbf{g}$ and $\mathbf{v}$, i.e. $p(\mathbf{p}|\mathbf{g},\mathbf{v},\mathbf{u})=p(\mathbf{p}|\mathbf{g},\mathbf{v})$. This structure is represented by the Bayesian network in Fig. 6 and requires the prediction matrix to have the following form

[TABLE]

It will be shown in the sequel that even with this restricted form it is not possible to overcome the simplifications and further constraints must be imposed.

6.2.1 Naive Scheme

Now, consider the naive calibration and prediction scheme introduced in 3. Recall that this method will be conservative for a suboptimal simplification when condition Eq. 32, introduced in 2, is satisfied. Specifically, this requires

[TABLE]

where $\mathbf{E}{}^{n}{}=\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}(\tilde{\mathbf{G}}\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}+\mathbf{\Sigma}_{\boldsymbol{\delta}})^{-1}$. For a prediction matrix $\mathbf{Y}$ with structure consistent with Eq. 42, this condition can be rewritten as

[TABLE]

Now, under non-trivial conditions, this holds when

[TABLE]

Condition Eq. 43a requires the random vector $\bar{\mathbf{b}}=\bar{\mathbf{B}}\mathbf{v}$ to be uncorrelated with $\tilde{\mathbf{g}}=\tilde{\mathbf{G}}\mathbf{v}$ under the prior covariance $\mathbf{\Sigma}_{\mathbf{v}}$. Furthermore, condition Eq. 43b requires $\tilde{\mathbf{G}}\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}\gg\mathbf{\Sigma}_{\boldsymbol{\delta}}$ and occurs when the data is perfect such that $\mathbf{\Sigma}_{\boldsymbol{\delta}}\approx 0$, or when the prior knowledge in the subspace that is informed by the data (i.e. $\operatorname{rowspace}(\tilde{\mathbf{G}})$) is very weak such that $\tilde{\mathbf{G}}\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}\approx\infty$.

These conditions ensure that: (i) the subset of parameters of the simple model that directly influence the predictions, represented by $\bar{\mathbf{b}}=\bar{\mathbf{B}}\mathbf{v}$, cannot be estimated from the data and (ii) the simple model can exactly reproduce the data.

It is important to note however that these conditions are unlikely to hold in practical scenarios and the naive scheme is not guaranteed to be conservative, even for the predictions with the structure considered in Eq. 41.

6.2.2 Data Driven Scheme

To overcome the above issues, it is proposed to ensure conservativeness through the judicial removal of information such that the conditions defined in Eq. 43 can be satisfied. This will occur in two main areas:

• •

An easily computable inflation of the prior covariance matrix $\mathbf{\Sigma}_{\mathbf{v}}$.

• •

The application of a preprocessing or filtering step that may discard or combine some components of the data.

Additionally, the class of predictions must be restricted further than those allowed by Eq. 41.

To start, let $\mathbf{F}\in\mathbb{R}^{D_{\mathbf{d}^{\prime}}\times D_{\mathbf{d}}}$ denote a data filtering matrix such that

[TABLE]

where $D_{\mathbf{d}^{\prime}}\leq D_{\mathbf{d}}$. Similarly define the transformed data matrices as $\mathbf{G}^{\prime}=\mathbf{F}\mathbf{G}$, and $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}\tilde{\mathbf{G}}$, and the transformed data error covariance as $\mathbf{\Sigma}_{\boldsymbol{\delta}}^{\prime}=\mathbf{F}\mathbf{\Sigma}_{\boldsymbol{\delta}}\mathbf{F}^{\top}{}$. The purpose of the filter $\mathbf{F}$ will be elaborated on later. Note however, filtering is optional and the identity matrix can be used $\mathbf{F}=\mathbf{I}$.

Now, consider the singular value decomposition of $\tilde{\mathbf{G}}^{\prime}$

[TABLE]

The last expression $\mathbf{U}_{1}\mathbf{S}_{1}\mathbf{V}_{1}^{\top}{}$ only includes the nonzero singular values, and will be referred to as the compact SVD. With this expansion, the vectors $\mathbf{v}_{1}=\mathbf{V}_{1}^{\top}{}\mathbf{v}$ and $\mathbf{v}_{2}=\mathbf{V}_{2}^{\top}{}\mathbf{v}$ define the rowspace and nullspace components of the parameters of the simple model for the given filtered data matrix $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}\tilde{\mathbf{G}}$. To meet condition Eq. 43b, the prior covariance will be inflated such that all information pertaining to the rowspace of $\tilde{\mathbf{G}}^{\prime}$, i.e. the subspace spanned by $\mathbf{V}_{1}$, is removed. In addition any correlation in the random vectors $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ will be ignored. Such a calibration and prediction scheme that embodies these properties is now defined.

Definition 7** (Data Driven Scheme).**

For a given filtering matrix $\mathbf{F}\in\mathbb{R}^{D_{\mathbf{d}^{\prime}}\times D_{\mathbf{d}}}$ with $D_{\mathbf{d}^{\prime}}\leq D_{\mathbf{d}}$, a data driven scheme is denoted by the set of simplified linear Gaussian probability density functions $\hat{p}{}^{d}(\mathbf{v})$, $\hat{p}{}^{d}(\mathbf{d}^{\prime}|\mathbf{v})$ and $\hat{p}{}^{d}(\mathbf{p}|\mathbf{v})$ parameterized by the matrices $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}$, $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}\tilde{\mathbf{G}}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}^{\prime}=\mathbf{F}\mathbf{\Sigma}_{\boldsymbol{\delta}}\mathbf{F}^{\top}{}$, $\tilde{\mathbf{Y}}$, and $\mathbf{\Sigma}_{\boldsymbol{\rho}}$.

To fully define the prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}$, let the covariance matrix $\mathbf{\Sigma}_{\mathbf{v}}$ denote the uncertainty in the parameters of the simplified model defined in 1 as $\mathbf{\Sigma}_{\mathbf{v}}=\mathbf{C}^{\dagger}{}\mathbf{\Sigma}_{\mathbf{x}}\mathbf{C}^{\dagger}{}{}^{\top}{}$. Furthermore, let the matrices $\mathbf{V}_{1}$ and $\mathbf{V}_{2}$ each denote a set of orthonormal column vectors that form a basis for the rowspace and nullspace of $\tilde{\mathbf{G}}^{\prime}$, e.g. as defined in Eq. 45. The prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}$ is now defined by the limit

[TABLE]

For a given dataset $\mathbf{d}$ the posterior belief over the predictions of interest generated by the data driven scheme is the Gaussian density $\hat{p}{}^{d}(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{d},\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{d})$, with mean and covariance given by

[TABLE]

*where the functions ${\mu}(\cdot)$ and $\mathit{\Sigma}(\cdot)$ are given in 1. *

Finally, a more specific class of predictions is considered than that defined in Eq. 42, such that condition Eq. 43a will hold by construction. Specifically, the prediction must depend on the uncorrupted version of the filtered data, denoted by $\mathbf{g}^{\prime}=\mathbf{F}\mathbf{G}\mathbf{x}$, or on the components of the parameters of the simple model that line in the nullspace of $\tilde{\mathbf{G}}^{\prime}$, denoted by $\mathbf{v}_{2}=\mathbf{V}_{2}^{\top}{}\mathbf{v}$. The prediction cannot depend on any rowspace components directly. Predictions with this structure can be written in the form

[TABLE]

where $\mathbf{A}$ and $\mathbf{B}$ are arbitrary matrices. Furthermore, the dependency structure of this restricted class of predictions is depicted in Fig. 7. The key difference between this and the structure depicted in Fig. 6 is the explicit separation of $\mathbf{v}$ into $\mathbf{v}_{1}$ and $\mathbf{v}_{2}$ and the restriction that $\mathbf{v}_{1}$ cannot have a direct dependency on the prediction. That is, the prediction $\mathbf{p}$ is conditionally independent of $\mathbf{u}$ and $\mathbf{v}_{1}$ given $\mathbf{g}^{\prime}$ and $\mathbf{v}_{2}$, or equivalently $p(\mathbf{p}|\mathbf{g}^{\prime},\mathbf{v}_{1},\mathbf{v}_{2},\mathbf{u})=p(\mathbf{p}|\mathbf{g}^{\prime},\mathbf{v}_{2})$. A prediction with this restricted structure has a prediction matrix of the form

[TABLE]

It is noted that the matrices $\mathbf{A}$ and $\mathbf{B}$ are arbitrary, the important characteristic is the structure that links the data and predictions to the simplification. This structure not only ensures that the simplifications are hidden behind the data, but also that any surrogate roles that the parameter vector $\mathbf{v}_{1}$ is forced to undertake during calibration cannot adversely affect the predictions.

It is now demonstrated that the data driven scheme is conservative for predictions that have this special structure.

Proposition 4** (Performance of Data Driven Scheme).**

Suppose the matrices $\mathbf{\Sigma}_{\mathbf{x}}$, $\mathbf{G}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}$, $\mathbf{Y}$, $\mathbf{\Sigma}_{\boldsymbol{\rho}}$ define an optimal scheme with posterior denoted by $p(\mathbf{p}|\mathbf{d})=N(\mathbf{p};\boldsymbol{\mu}{}_{\mathbf{p}|\mathbf{d}},\mathbf{\Sigma}_{\mathbf{p}|\mathbf{d}})$ for arbitrary dataset $\mathbf{d}$.

Furthermore, let $\mathbf{C}$ denote a suboptimal simplification matrix and let the matrices $\mathbf{F}$, $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}$, $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}\mathbf{G}\mathbf{C}$, $\mathbf{\Sigma}_{\boldsymbol{\delta}}^{\prime}=\mathbf{F}\mathbf{\Sigma}_{\boldsymbol{\delta}}\mathbf{F}^{\top}{}$, $\tilde{\mathbf{Y}}=\mathbf{Y}\mathbf{C}$, and $\mathbf{\Sigma}_{\boldsymbol{\rho}}$ denote a data driven scheme, where $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}$ is as specified in 7. Additionally, let $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{d}$ and $\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{d}$ denote the mean and covariance of the posterior produced by this scheme for the dataset $\mathbf{d}$.

If the filtering matrix $\mathbf{F}$ is selected such that $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}\tilde{\mathbf{G}}$ has full row rank, then the estimator matrix $\mathbf{E}{}^{d}$ of the data driven scheme is equivalent to the pseudoinverse of $\tilde{\mathbf{G}}^{\prime}$ and can be expressed in terms of the compact SVD

[TABLE]

In addition, the posterior mean and covariance can be rewritten as

[TABLE]

where $\mathbf{W}=\mathbf{I}-\mathbf{V}_{1}\mathbf{V}_{1}^{\top}{}$. 2. 2.

If $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}\tilde{\mathbf{G}}$ has full row rank, and the high fidelity prediction matrix has the form

[TABLE]

then the data driven scheme is conservative with respect to the optimal scheme.

See Appendix E for a proof.

Before providing additional intuition as to the importance of this result, and the role of the filtering term $\mathbf{F}$, it is first demonstrated that the data driven scheme is a generalization of the truncated singular value decomposition calibration method.

6.2.3 Truncated SVD Inversion

The truncated singular value decomposition inversion or calibration method is a commonly used technique aster_parameter_2012 ; moore_role_2005 ; white_quantifying_2014 that generates the solution to the parameter estimation problem using a truncated decomposition of the data matrix. It will be demonstrated that the truncated SVD scheme can be replicated with the appropriate selection of a filtering matrix. To start, consider the SVD of the simplified data matrix

[TABLE]

Now, let $\tilde{\mathbf{U}}_{t}$ denote the columns of $\tilde{\mathbf{U}}$ that correspond to the largest $k$ nonzero singular values. Furthermore, consider the filtering matrix defined by

[TABLE]

With this definition, the filtered data matrix $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}_{\text{TSVD}}\tilde{\mathbf{G}}$ becomes

[TABLE]

where $\tilde{\mathbf{V}}_{t}$, $\tilde{\mathbf{S}}_{t}$ similarly denote truncated versions of $\tilde{\mathbf{V}}$, $\tilde{\mathbf{S}}$. Now, as $\mathbf{F}_{\text{TSVD}}\tilde{\mathbf{G}}$ has full row rank (only nonzero singular values are included), the results of 4(1) allow the estimator matrix for the scheme to be given by the pseudoinverse of $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}_{\text{TSVD}}\tilde{\mathbf{G}}$ and thus can be written as

[TABLE]

The estimated parameter vector of the simple model is now related to the data vector via the inverse of the truncated data matrix $\tilde{\mathbf{G}}$, i.e.

[TABLE]

This demonstrates that the truncated SVD inversion method can be replicated with the selection of an appropriate filtering matrix.

It is important to note however that this does not mean that the predictions generated by the truncated SVD method will be conservative. To guarantee this, the prediction must have the form specified in condition Eq. 49 of 4, which in this case requires the prediction matrix to have the form

[TABLE]

where $\tilde{\mathbf{V}}_{\bar{t}}$ contains the columns of $\tilde{\mathbf{V}}$ that were removed by the truncation. If this is not obeyed, the predictions may be over confident. Furthermore, this condition is implicitly dependent on the truncation point $k$, such that this may hold for some values and not others. This dependency explains in part the results of the simulation studies performed by white_quantifying_2014 which demonstrated the difficulty in choosing an appropriate truncation point such that predictive uncertainty is accurately estimated by a simplified model.

6.2.4 Selection of Data Filtering

A key element of the defined data driven prediction scheme is the filtering matrix $\mathbf{F}$. For the results of 4 to hold, and the scheme to be conservative, two main requirements must be met:

$\mathbf{F}$ must be selected such that $\tilde{\mathbf{G}}^{\prime}=\mathbf{F}\tilde{\mathbf{G}}$ has linearly independent rows. This is a fairly trivial requirement, and simply requires that dependent, e.g. duplicated, measurements be combined, e.g. by averaging. Note that the information content is preserved as the error covariance is also transformed. In addition, it requires measurements that are insensitive to the model parameters be dropped. 2. 2.

$\mathbf{F}$ must be chosen such that the induced separation of the parameter vector $\mathbf{v}$ into those components that are estimatable $\mathbf{v}_{1}=\mathbf{V}_{1}\mathbf{v}$ and lie in the rowspace of $\mathbf{F}\tilde{\mathbf{G}}$ from those that are not $\mathbf{v}_{2}=\mathbf{V}_{2}\mathbf{v}$ and lie in the nullspace of $\mathbf{F}\tilde{\mathbf{G}}$, forces all components of the parameter vector that have a direct influence on the prediction to be contained in the vector $\mathbf{v}_{2}$ and are thus not updated during the calibration process. This ensures that these parameters do not take on any surrogate roles, and that the predictions are not corrupted.

In addition to these two requirements, the filtering matrix can be selected to improve the predictive performance by removing data components that are only weakly informative to the parameters. This can be considered identical to selecting an appropriate truncation point on a SVD calibration scheme such that the highly informative prior knowledge can be used instead of data that only imposes weak constraints, e.g. as recommended by moore_role_2005 .

However before searching for an optimal filter, it is important to note that for a given type of data, predictions, and simplification, defined by the triple $\mathbf{G},\mathbf{Y},\mathbf{C}$, there may not be an $\mathbf{F}$ that allows the prediction to be cast in the form required by 4 and the use of the data driven scheme is not guaranteed to be conservative. This reinforces the fact that choosing a model simplification and performing calibration with a particular dataset, must be considered an integrated task that is explicitly dependent on the types of predictions that are required.

6.2.5 Summary

This section has introduced a calibration and prediction scheme that is a generalization of the typically used truncated SVD scheme. Furthermore it has been demonstrated that the scheme is guaranteed to produce a conservative prediction posterior when the unmodelled complexity is hidden behind the data. The general structure of the prediction problem, that obeys this condition, is depicted in graphically in Fig. 8. This represents the probability densities that define the data driven scheme and includes the dependency structure of the predictions (as required by condition Eq. 49 and depicted in Fig. 7) as well as the structure of the inflated prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}$ (as specified in Eq. 46).

The important contribution here is the derived structural condition between the data, predictions and the simplifications that determines when it is adequate to use the data driven scheme, with a simplified model, to assess predictive uncertainty. It is noted that unlike the optimally compensated scheme, this scheme can be used without needing the explicit numerical values within the data and prediction matrices of the high fidelity reference model. It is noted however that, structural information from them is still required to ensure 4 is met. Although this still may not be straightforward in practice, structural considerations are generally significantly easier to handle than numerical ones.

7 Discussion and Conclusions

This paper has considered what constitutes a simplified, but useful, model. In particular it has examined how simplified models can be used to combine data and expert knowledge within a calibration or inversion process to generate a prediction with a conservative estimate of uncertainty. The concept of an optimal simplified model was defined that determines when standard probabilistic calibration methods are adequate to quantify predictive uncertainty. The main contribution is the introduction of two new calibration and prediction schemes, along with conditions that explicitly define when they are appropriate to generate a conservative estimate of uncertainty, for suboptimal models. These conditions explicitly relate the nature of the calibration data, the predictions of interest as well as the simplifications within the model.

The first scheme allows the optimal posterior distribution to be generated by allowing the simplifications to be overcome by adjusting how the beliefs of the modeller are used to define a prior term. The scheme is only applicable for highly parameterized models. Furthermore, it requires the prior covariance for the parameters of the simple model to be generated using the data and prediction matrices of a high fidelity reference model. This is a significant limitation as if these matrices could be generated for a practical problem they could be directly used to produce a prediction posterior without the need to consider a simplified model. Nevertheless, the value of this scheme lies in defining the ideal calibration process for a simplified model and demonstrates it is possible to overcome suboptimal simplifications through the judicial selection of a prior regularization matrix.

The second data driven scheme is designed for predictions that are strongly related to the data, such that the unmodelled complexity effects both in the same way. This scheme does not require the data and prediction matrices of the high fidelity model to be available for model calibration in the way the optimally compensated scheme does. However, it is only guaranteed to be conservative if the predictions, data and simplification have the required structural form. The key insight provided by this scheme is understanding how the model simplifications can be overcome with the use of the right calibration data.

It was also demonstrated that this data driven scheme is a generalization of the popular truncated singular value decomposition inversion scheme aster_parameter_2012 . The generalization allows greater flexibility in filtering the data to ensure that the predictions are of the appropriate form for a given simplified model. Both calibration schemes have been applied to a prototypical groundwater prediction problem in Appendix A.

Finally, it is noted that each of the two newly defined calibration and prediction schemes have conditions that are linked to the data and prediction matrices of a high fidelity reference model. For any practical problem it is unlikely that these conditions can be directly assessed and further subjective judgment will be needed. To make this process easier, two areas of further work can be pursued. Firstly, more synthetic experimental analyses are required to demonstrate how the two schemes can be applied in more complex problems, e.g. as in white_quantifying_2014 . The second area is to understand the effect of partial non-satisfaction of these conditions, and determine how the schemes can be made robust to these.

It is also noted that the use of Bayesian networks, for instance as depicted in Figure 10, may be of great benefit to determine when the required structural conditions are likely to be obeyed, for instance, when is a given simplified model optimal, or when is it appropriate to use the data driven scheme. Such structural considerations may also help frame the arguments put forward in modelling projects (e.g. performed within environmental impact assessment studies) as to why a given modelling and calibration approach is adequate for a given prediction problem. For instance these argument can be framed using a two stage process. The first may conceptually link the complex system to an optimally simplified model, this stage would consider the features and characteristics of the system that ideally should be modelled. The second stage may then put forward arguments as to how any further simplifications used to produce a suboptimal numerical model will be handled by the calibration process.

Lastly, several other important areas of further research are identified:

• •

Nonlinear simulators. The analysis within this paper has required a linear relationship between the system properties and the data and predictions. This has allowed generic insights to be obtained, but is a significant limitation and further work is needed to relax this. One approach is to consider higher order expansions of the models such that some of the nonlinearities can be included, for example second order expansions have been considered in box_bias_1971 ; mclaughlin_distributed_1988 ; cooley_bias_2006 . Alternatively, more direct probabilistic formulations may also be possible that e.g. generalize the data driven scheme and only exploit the structural constraints with the problem.

• •

Nonlinear parameterizations. In the developed analytical framework only linear relationships between the parameters of the high fidelity model and the simplified model where considered. This should be extended to considered nonlinear parameterizations.

• •

Over constrained calibration problems. The results obtained for the optimally compensated prediction scheme required that the simplified model is under constrained by the data and provides no insight to aid the calibration of over constrained problems. Nevertheless, other approaches may be possible that reproduce the optimal, or at least a conservative, result. Similar modifications may also be possible for the data driven scheme.

• •

Conservativeness. This paper has employed a definition of conservativeness based on the mean and covariance of the model predictions. Generalizations of this to non-Gaussian distributions have been considered in e.g. bailey_conservative_2012 ; ajgl_conservativeness_2013 . However, it is perhaps more appropriate to consider a generalization of conservativeness that explicitly includes the subsequent decision problem (i.e. engineering design or environmental management). This can be performed using a decision theoretic approach berger_statistical_1985 that includes the utility function of this subsequent decision problem and defines an approximate density $\hat{p}(\mathbf{p})$ as conservative if the expected utility is not over estimated, e.g.

[TABLE]

Here $U(\mathbf{a},\mathbf{p})$ is the utility function that encodes the gain under action $\mathbf{a}$, when the consequence $\mathbf{p}$ occurs. Such a generalization would also allow the incorporation of the decision problem into how modelling and calibration should be performed. It is noted that the above generalization is equivalent to 4 when $U(\mathbf{a},\mathbf{p})$ has the form of a weighted squared difference, i.e. $U(\mathbf{a},\mathbf{p})=-(\mathbf{a}-\mathbf{p})^{\top}{}\mathbf{A}(\mathbf{a}-\mathbf{p})$ for arbitrary positive semidefinite matrix $\mathbf{A}$.

Appendix A Example: Groundwater Head Prediction

The two newly defined calibration and prediction schemes are now applied to a prototypical groundwater prediction problem similar to that considered in white_quantifying_2014 . The scenario is depicted in Fig. 9 and consists of a 1D confined aquifer with a single observation well. Of interest is the predicted hydraulic head within the aquifer to the right of the observation well. It is noted that this is a very rudimentary problem, however the objective here is to demonstrate the differences between the schemes and how the conditions that guarantee conservativeness relate to a specific problem.

A.1 Prior Beliefs

For this example, the expert modeller has the following beliefs about the system (taken from white_quantifying_2014 where possible)

A constant head boundary condition is believed to exist at the left side, corresponding to a discharge location. The head at this location, $h_{0}$, is believed to be normally distributed with mean $1.0$m and standard deviation $0.75$m above the upper confining layer. 2. 2.

No areal recharge or leakage is believed to take place within the domain of interest. 3. 3.

The thickness of the aquifer $b$ is known to be a constant over the domain with value $10$m. 4. 4.

The rate $q$ of water flowing through the aquifer known to be $0.5$m3/day. 5. 5.

The system is believed to be in steady state. Furthermore, a numerical simulation of the system with cell length $\ell=10$m is believed to be adequate to capture the spatial variation in the head field that is of interest. 6. 6.

The hydraulic conductivity of the aquifer is believed to be heterogeneous, with a mean value of $2.5$m/day. A set of 10 cells, each of length $\ell=10$m, is considered adequate to describe the heterogeneity, where the hydraulic conductivity of each cell, $K_{i}$, is believed to be log normally distributed such that $\log_{10}K_{i}$ has a mean $\log_{10}2.5\approx 0.398$, and a spatial correlation described by the exponential variogram with sill of $0.1\approx(0.316)^{2}$ and range $300$m. 7. 7.

The error in the data acquisition method is believed to be small and normally distributed with a mean of zero and standard deviation $0.1$m.

It is noted that these capture what is believed to realistically describe the system (or at least an optimally simplified model). Further simplifying assumptions will be considered in Section A.4.

Based on the above prior knowledge, the system properties can be defined as the vector that captures the constant boundary head, and the set of hydraulic conductivities

[TABLE]

Furthermore, the uncertainty in the system properties is fully captured by the mean and covariance matrix that defines the Gaussian distribution $p(\mathbf{x})=N(\mathbf{x};\boldsymbol{\mu}{}_{\mathbf{x}},\mathbf{\Sigma}_{\mathbf{x}})$.

To define the data and prediction likelihood functions, consider the vector $\mathbf{x}$ as known and define the mean with the following nonlinear functions, corresponding to the finite difference solution to Darcy’s equations bear_modeling_2010

[TABLE]

Note the variables $\ell$, $b$ and $q$ are all known constants. Furthermore, as the numerical simulation at this discretization is believed to be adequate, the data error $\boldsymbol{\delta}$ is completely captured by errors in the measurement process such that $\mathbf{\Sigma}_{\boldsymbol{\delta}}=(0.1)^{2}$. In addition the prediction error variance is taken to be zero, $\mathbf{\Sigma}_{\boldsymbol{\rho}}=0$.

It is noted that the prediction equation has been rewritten on line Eq. 59 to explicitly include the data equation as an additive term. This form will be important when judging whether the data driven scheme is appropriate for a particular system simplification.

A.2 Data Generation

The measured head at the observation well is $\mathbf{d}=2.5$m. This is generated by the data equations $\mathbf{d}=\mathcal{G}(\mathbf{x}_{t})$ with no added measurement error. In addition the system properties $\mathbf{x}_{t}$ are the same as the prior mean $\boldsymbol{\mu}{}_{\mathbf{x}}$, with the exception that the boundary head is set to $h_{0}=1.5$m. This difference corresponds to $2/3$ of the prior standard deviation.

A.3 Linearized Solution

From the nonlinear functions $\mathcal{G}(\mathbf{x})$ and $\mathcal{Y}(\mathbf{x})$, a pair of linear functions are constructed by linearizing about the prior mean with a first order Taylor expansion mclaughlin_reassessment_1996 ; carrera_inverse_2005 .

[TABLE]

and similarly for $\mathcal{Y}$. Furthermore, a transform into zero mean increments is performed such that $\Delta\mathbf{x}=\mathbf{x}-\boldsymbol{\mu}{}_{\mathbf{x}}$, $\Delta\mathbf{d}=\mathbf{d}-\mathcal{G}(\boldsymbol{\mu}{}_{\mathbf{x}})$, and $\Delta\mathbf{p}=\mathbf{p}-\mathcal{Y}(\boldsymbol{\mu}{}_{\mathbf{x}})$. Under the linearized approximation the data and prediction equations become

[TABLE]

where the data and prediction matrices are defined by the Jacobian matrices evaluated at the prior mean $\mathbf{G}=\nabla_{\mathbf{x}}\mathcal{G}(\boldsymbol{\mu}{}_{\mathbf{x}})$ and $\mathbf{Y}=\nabla_{\mathbf{x}}\mathcal{Y}(\boldsymbol{\mu}{}_{\mathbf{x}})$.

The solutions to the full non-linear prediction problem and the approximate linearized version defined above are given in Fig. 11(a). The non-Gaussian posterior density of the nonlinear problem has been produced using a standard Metropolis MCMC algorithm brooks_handbook_2011 to generate samples that are then smoothed using a kernel density estimator scott_multivariate_1992 . Overlaid with this is the Gaussian posterior produced from the linearized problem. It is noted that the posterior under the non-linear equations is slightly more peaked and non-symmetric when compared to the posterior generated by the linearized scheme.

It is noted that the non-linearity is only present as the log conductivity is considered to be normally distributed. If, for instance, the hydraulic resistance (inverse conductivity) was represented directly, the equations would become linear. However, this is not pursued here and the error caused by assuming a linear approximation of the nonlinear equations is considered out of scope.

A.4 Assumed Simplifications

Consider the following simplifying assumptions

1. A1:

The constant flow boundary condition is assumed known and set to the prior value of $1.0$m. 2. A2:

Cells 1-5 are grouped together and assumed a homogeneous unit. This set of cells will be referred to as Zone A. 3. A3:

Cells 6-10 are grouped together and assumed a homogeneous unit. This set of cells will be referred to as Zone B.

This enables the parameters of the simple model to be defined by the vector of log conductivities of the two zones

[TABLE]

This simplification scheme corresponds to a $\mathbf{C}$ matrix of

[TABLE]

The use of increments allows the state of the high fidelity model to be defined in terms of an increment $\Delta\mathbf{v}$ such that

[TABLE]

Furthermore, the simplified data equation becomes

[TABLE]

and similarly for the prediction equation. Now, the linearized form of the simplified data and prediction equations become

[TABLE]

Here the simplified data and prediction matrices can be written in terms of the high fidelity model $\tilde{\mathbf{G}}=\mathbf{G}\mathbf{C}$ and $\tilde{\mathbf{Y}}=\mathbf{Y}\mathbf{C}$, where $\mathbf{G}$ and $\mathbf{Y}$ are the Jacobian matrices of the high fidelity model $\mathbf{G}=\nabla_{\mathbf{x}}\mathcal{G}(\boldsymbol{\mu}{}_{\mathbf{x}})$, $\mathbf{Y}=\nabla_{\mathbf{x}}\mathcal{Y}(\boldsymbol{\mu}{}_{\mathbf{x}})$.

It is now of interest to determine whether this simplification is optimal. If it is, the naive prediction scheme when used with the simplified model will produce the same results as the optimal scheme applied to the high fidelity model. For the simplification to be optimal, the unmodelled complexity should have no effect on the data or predictions. This occurs when the following conditions hold

[TABLE]

where the columns of $\mathbf{U}_{\bar{\mathbf{C}}}$ define an orthonormal basis for the cokernel of $\mathbf{C}$, as defined in Eq. 17. The unmodelled complexity consists of $h_{0}$ and the two 4D parameters that describe the small scale heterogeneity in the conductivity of the two zones. It is noted that this condition does not hold and thus the simplification is not optimal. Furthermore, the failure of this condition to hold is completely due to assumption A1 that considers the boundary condition $h_{0}$ is known. If this assumption is removed and $h_{0}$ included in the simplified model, while at the same time retaining the homogeneity assumptions of A2 and A3, then the simplification would become optimal (under the linearization considered). This optimal simplification will not be considered and it will be expected that the naive method will not be conservative.

Finally, the structure of the high fidelity prediction problem and the simplified problem is depicted graphically with two Bayesian networks in Fig. 10. It is noted that the suboptimality of the simplification can be observed in Fig. 10(b) as the data and predictions are not conditionally independent of the unmodelled components ($h_{0}$) given the model parameters ($K_{A}$ and $K_{B}$).

A.5 Prediction Results

Here, the performance of the two new calibration and prediction schemes are considered, with the resulting posterior density functions for the predicted head displayed in Fig. 11(b). The results are also compared with the naive scheme which represents a typical probabilistic calibration scheme of a simplified model. It is important to note that no “ground truth” value for the prediction is given for comparison purposes as the objective is to produce the full probability distribution and not make a single point prediction.

A.5.1 Naive Scheme

The naive scheme considers the uncertainty in the parameters included by the simple model, but ignores any errors introduced in the simplification. As the simplification is suboptimal, it should not be expected that the prediction posterior will be conservative. Furthermore, the results in Fig. 11(b) demonstrate that this is true for this scenario as the posterior is overly narrow.

To understand the cause of the overconfidence in this scenario, consider the prior covariance in the parameters $\Delta\mathbf{v}$ as defined in 1

[TABLE]

Due to the large range of 300m used for the variogram of the log conductivity field, the correlation between the two parameters maintained by the simple model is considerable. In addition, the known flow rate and the assumption of known boundary head $h_{0}$, allows the head measurement to be used to very accurately, but incorrectly, estimate the conductivity $K_{A}$. Furthermore, due to the large correlation that $K_{A}$ has with $K_{B}$, this scheme produces an overconfident prediction probability density function for the parameter $K_{B}$, which in turn causes the predicted head to be also overconfident and non-conservative.

A.5.2 Optimally Compensated Scheme

Now consider the optimally compensated calibration and prediction scheme. Firstly it is noted that the simplified problem is highly parameterized as $D_{\mathbf{v}}\geq D_{\mathbf{p}}+D_{\mathbf{d}}$ as $D_{\mathbf{v}}=2$ and $D_{\mathbf{p}}=D_{\mathbf{d}}=1$. The major difference between the naive scheme and the optimally compensated scheme is the modification of the prior covariance used for the parameters of the simple model. In particular, the modified prior of this scheme is not just defined in terms of the simplification matrix, but also in terms of the data and prediction matrices, via the matrix $\mathbf{R}=\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}$, where \mathbf{Z}=\bigl{[}\begin{smallmatrix}\mathbf{G}\\ \mathbf{Y}\end{smallmatrix}\bigr{]} and \tilde{\mathbf{Z}}=\bigl{[}\begin{smallmatrix}\tilde{\mathbf{G}}\\ \tilde{\mathbf{Y}}\end{smallmatrix}\bigr{]}. The optimally compensating prior covariance becomes

[TABLE]

It is noted that the only difference between the values of the naive prior covariance $\mathbf{\Sigma}_{\mathbf{v}}$ and $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$, is the inflation of the marginal variance of the first component corresponding to $\log K_{A}$. This causes less information to be propagated from the data into $\log K_{B}$, and enables the prediction posterior to replicate the results of the high fidelity model.

The key difficulty with this scheme is that it requires the data and prediction matrices of the high fidelity model in order to produce the inflated prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ that compensates for the suboptimal simplifications.

A.5.3 Data Driven Scheme

To apply the data driven calibration and scheme, it is first of interest to determine if it is appropriate for the problem, in particular, does the prediction have the right structural relationship with the data and the simplifications as specified in 4. This will be checked first using the algebraic expression, and then using the structural characteristics of the network depicted in Fig. 10(b).

The required algebraic condition requires that for some $\mathbf{A}$ and $\mathbf{B}$ the prediction matrix can be written as

[TABLE]

For this problem, the data is scalar, and thus the filtering matrix will be set to unity $\mathbf{F}=1$.

Also, from the prediction equation in Eq. 59, it is clear that the prediction matrix $\mathbf{Y}$ is equal to the data matrix $\mathbf{G}$ plus an additional term, which will be denoted by $\mathbf{T}$. This allows the Jacobian matrix $\mathbf{Y}$ to be written as

[TABLE]

Thus, the matrix $\mathbf{A}$ in Eq. 61 will be set to unity $\mathbf{A}=1$.

Now it is of interest to consider the second term in Eq. 61. For the prediction to have the required form, the matrix $\mathbf{T}=\mathbf{Y}-\mathbf{G}$ must be equivalent to $\mathbf{B}\mathbf{V}_{2}^{\top}{}\mathbf{C}^{\dagger}{}$, and thus must satisfy two conditions:

[TABLE]

The first condition (63) requires the difference between the data and the prediction to be independent of the unmodelled complexity. For this scenario this can be checked algebraically using the high fidelity prediction and data matrices and is satisfied. Now consider the second condition (64), this requires the difference between the data and the prediction to be independent on the rowspace components of the $\tilde{\mathbf{G}}$. This can be easily checked using just the simplified data and prediction matrices and is also satisfied. Thus, it is guaranteed that the data driven scheme is conservative for this problem.

As an alternative to this algebraic check, the conditional independence requirements can be check using the structure of the Bayesian network. In particular, the requirement of Eq. 61 is equivalent to the requirement that the prediction $\mathbf{p}=h_{10}$ is conditionally independent of all unmodelled components, $\mathbf{u}=[h_{0},\mathbf{K}_{A}^{-},\mathbf{K}_{B}^{-}]$, given the uncorrupted data $\mathbf{g}=h_{5}$ and the nullspace components of the parameter vector $\mathbf{v}_{2}=K_{B}$. From the structure of the network, these two nodes are the only parents of the prediction and thus, the required independence requirement is satisfied. This graphical approach provides much greater intuition into when the requirements are met.

When the data driven scheme is applied, it is noted that the available prior knowledge is modified in two ways. Firstly, the prior information on $K_{A}$ is ignored and this parameter is estimated from the data only. Secondly, the correlation between $K_{A}$ and $K_{B}$ is also ignored, this causes the posterior over the parameter $K_{B}$ to be the same as the prior. These modifications are encapsulated in the data driven prior covariance matrix

[TABLE]

The posterior generated by this scheme is conservative and is depicted in Fig. 11(b). It is noted that the posterior is slightly shifted than that produced by the optimal scheme and has a larger variance.

A.6 Summary

This section has demonstrated the operation of the naive, optimally compensated, and data driven calibration and prediction schemes on a prototypical groundwater problem. A high fidelity model was considered to represent the true belief of an modeller and a suboptimal simplification developed to represent a computational model.

The key difference between the three schemes is in the specification of the prior covariance matrix for the parameters of the simplified model. In particular, the naive scheme generates this by directly projecting the prior distribution that exists in the parameters of the high fidelity model into the parameters of the simplified model. The optimally compensated scheme performs the projection using full knowledge of the data and prediction matrices of the high fidelity model. Lastly, the data driven scheme uses the same matrix as the naive scheme but throws away some of the information it contains.

The conditions under which the schemes are conservative (or optimal) have also been highlighted. The naive scheme should only be applied to optimally simplified models, this was not satisfied and the posterior shown to be overconfident and non-conservative. The optimally compensated scheme should only be applied to highly parameterized models. And the data driven scheme is conservative only when the predictions and data have a similar dependency to the unmodelled complexity. These last two conditions were shown to hold for the problem considered, and the generated posterior distributions were also conservative.

Appendix B Proof of Proposition 1

Proof.

For a given simplification matrix, the parameters of the high fidelity model $\mathbf{x}$ can be expanded using Eq. 18 as

[TABLE]

Furthermore, from the decomposition in Eq. 17, the matrix $\begin{bmatrix}\mathbf{C}&\mathbf{U}_{\bar{\mathbf{C}}}\end{bmatrix}$ can be rewritten as

[TABLE]

This is nonsingular and has an inverse that simplifies to

[TABLE]

where $\mathbf{C}^{\dagger}{}=\mathbf{V}_{\mathbf{C}}\mathbf{S}_{\mathbf{C}}^{-1}{}\mathbf{U}_{\mathbf{C}}^{\top}{}$. Thus, the transformed covariance matrix in the space of $\mathbf{v}$, $\mathbf{u}$ has the form

[TABLE]

Appendix C Proof of Proposition 2

Proof.

Under the conditions of part (1) the matrix $\mathbf{C}$ denotes an optimal simplification and thus $\mathbf{G}\mathbf{U}_{\bar{\mathbf{C}}}=\mathbf{Y}\mathbf{U}_{\bar{\mathbf{C}}}=0$, where $\mathbf{U}_{\bar{\mathbf{C}}}$ is an orthonormal basis for the cokernel of $\mathbf{C}$, e.g. as given in Eq. 17. Using the expansion of the parameter vector as $\mathbf{x}=\mathbf{C}\mathbf{v}+\mathbf{U}_{\bar{\mathbf{C}}}\mathbf{u}$, the covariance $\mathbf{\Sigma}_{\mathbf{x}}$ becomes

[TABLE]

Substituting Eq. 65 into the posterior mean and covariance of the optimal scheme, defined in Eq. 10 and Eq. 11, and noting that $\mathbf{G}\mathbf{U}_{\bar{\mathbf{C}}}=\mathbf{Y}\mathbf{U}_{\bar{\mathbf{C}}}=0$, the following forms are obtained

[TABLE]

Noting that $\tilde{\mathbf{G}}=\mathbf{G}\mathbf{C}$ and $\tilde{\mathbf{Y}}=\mathbf{Y}\mathbf{C}$ these expressions are equivalent to mean and covariance those of the naive scheme, i.e. $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{n}=\boldsymbol{\mu}{}_{\mathbf{p}|\mathbf{d}}$ and $\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{n}=\mathbf{\Sigma}_{\mathbf{p}|\mathbf{d}}$ for all $\mathbf{d}$. This proves part (1) of the proposition.

To proceed, to parts (2) and (3), for the naive scheme to be conservative, the mean and covariance of the posterior must obey Eq. 30, which simplifies to

[TABLE]

Also, the expectation over $\mathbf{d},\mathbf{p}$ is equivalent to an expectation over the independent variables $\mathbf{x},\boldsymbol{\delta},\boldsymbol{\rho}$, where $\mathbf{d}=\mathbf{G}\mathbf{x}+\boldsymbol{\delta}$ and $\mathbf{p}=\mathbf{Y}\mathbf{x}+\boldsymbol{\rho}$. Thus to prove part (2) of the proposition, it will be now shown that this holds under the special condition of $\mathbf{Y}\mathbf{U}_{\bar{\mathbf{C}}}=\tilde{\mathbf{Y}}\mathbf{E}{}^{n}{}\mathbf{G}\mathbf{U}_{\bar{\mathbf{C}}}$, defined in Eq. 32.

To start, note that $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{n}=\tilde{\mathbf{Y}}\mathbf{E}{}^{n}{}\mathbf{d}$, $\mathbf{d}=[\mathbf{G}\mathbf{x}+\boldsymbol{\delta}]$ and $\mathbf{x}=\mathbf{C}\mathbf{v}+\mathbf{U}_{\bar{\mathbf{C}}}\mathbf{u}$, thus the difference $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{n}-\mathbf{p}$ can be written

[TABLE]

Where, Eq. 68 has used condition Eq. 32. Now, as $\mathbf{v}=\mathbf{C}^{\dagger}{}\mathbf{x},\boldsymbol{\delta},\boldsymbol{\rho}$ are all independent and zero mean, the expected squared difference of $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{n}-\mathbf{p}$ can be written in terms of the covariance matrices

[TABLE]

Expanding the right hand side with the naive estimator matrix $\mathbf{E}{}^{n}{}=\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}(\tilde{\mathbf{G}}\mathbf{\Sigma}_{\mathbf{v}}\tilde{\mathbf{G}}^{\top}{}+\mathbf{\Sigma}_{\boldsymbol{\delta}})^{-1}$ and simplifying, produces the result

[TABLE]

Thus, the special condition $\mathbf{Y}\mathbf{U}_{\bar{\mathbf{C}}}=\tilde{\mathbf{Y}}\mathbf{E}{}^{n}{}\mathbf{G}\mathbf{U}_{\bar{\mathbf{C}}}$ is sufficient for Eq. 66 to be satisfied (with equality), ensuring the naive scheme is conservative. This proves part (2).

To prove part (3), consider equation Eq. 67 above and assume $\mathbf{v}$ and $\mathbf{u}$ are independent and $\tilde{\mathbf{Y}}\mathbf{E}{}^{n}{}\mathbf{G}\mathbf{U}_{\bar{\mathbf{C}}}\neq\mathbf{Y}\mathbf{U}_{\bar{\mathbf{C}}}$, then the expected squared difference of $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{n}-\mathbf{p}$ becomes

[TABLE]

Defining $\mathbf{M}=\tilde{\mathbf{Y}}\mathbf{E}{}^{n}{}\mathbf{G}\mathbf{U}_{\bar{\mathbf{C}}}-\mathbf{Y}\mathbf{U}_{\bar{\mathbf{C}}}\neq 0$, this simplifies to

[TABLE]

Now, as $\mathbf{M}$ is a non-zero matrix and $\mathbf{\Sigma}_{\mathbf{u}}$ is a non-zero positive semi-definite covariance matrix, $\mathbf{M}\mathbf{\Sigma}_{\mathbf{u}}\mathbf{M}^{\top}{}$ must be non-zero and positive semi-definite. Thus, the required condition Eq. 66 does not hold as

[TABLE]

and the scheme is strictly not conservative. This proves part (3).

Appendix D Proof of Proposition 3

Proof.

To prove optimality of the compensated scheme it is necessary to show the conditions defined by Eq. 34 are satisfied by the prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}$ defined by this scheme. It is noted that these are satisfied when

[TABLE]

These conditions can be combined and rewritten as

[TABLE]

where \mathbf{Z}=\bigl{[}\begin{smallmatrix}\mathbf{G}\\ \mathbf{Y}\end{smallmatrix}\bigr{]} and \tilde{\mathbf{Z}}=\bigl{[}\begin{smallmatrix}\tilde{\mathbf{G}}\\ \tilde{\mathbf{Y}}\end{smallmatrix}\bigr{]}. Now with the definition of $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{o}=\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}\mathbf{\Sigma}_{\mathbf{x}}\tilde{\mathbf{Z}}^{\dagger}{}{}^{\top}{}\mathbf{Z}^{\top}{}$, condition Eq. 69 becomes

[TABLE]

and holds when $\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}=\mathbf{Z}$.

Now, to demonstrate $\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}=\mathbf{Z}$ is satisfied, the condition $\operatorname{rank}(\mathbf{Z}\mathbf{C})=\operatorname{rank}(\mathbf{Z})$ implies that $\operatorname{columnspace}(\tilde{\mathbf{Z}})=\operatorname{columnspace}(\mathbf{Z}\mathbf{C})=\operatorname{columnspace}(\mathbf{Z})$ (seber_matrix_2008, , 3.16). Furthermore, it is noted that $\operatorname{columnspace}(\tilde{\mathbf{Z}})=\operatorname{columnspace}(\tilde{\mathbf{Z}}^{\dagger}{}{}^{\top}{})$ (seber_matrix_2008, , 7.52(l)), and thus

[TABLE]

This implies that $\mathbf{Z}^{\top}{}(\mathbf{I}-\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\dagger}{})=0$ (seber_matrix_2008, , 2.34). Furthermore, as $\tilde{\mathbf{Z}}$ is real, $(\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\dagger}{})^{\top}{}=\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\dagger}{}$. Thus, taking the transpose of $\mathbf{Z}^{\top}{}(\mathbf{I}-\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\dagger}{})=0$, produces

[TABLE]

and hence $\tilde{\mathbf{Z}}\tilde{\mathbf{Z}}^{\dagger}{}\mathbf{Z}=\mathbf{Z}$.

Appendix E Proof of Proposition 4

Proof.

To start, note that the inverse of the prior covariance matrix $\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}$ simplifies as follows

[TABLE]

Now, the estimator matrix is defined as $\mathbf{E}{}^{d}=\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}\tilde{\mathbf{G}}^{\prime}{}^{\top}{}(\tilde{\mathbf{G}}^{\prime}\hat{\mathbf{\Sigma}}_{\mathbf{v}}^{d}\tilde{\mathbf{G}}^{\prime}{}^{\top}{}+\mathbf{\Sigma}_{\boldsymbol{\delta}}^{\prime})^{-1}$. Using the matrix inversion identity (seber_matrix_2008, , 15.1(a)), the inverse definition in Eq. 71, the compact SVD of $\tilde{\mathbf{G}}^{\prime}=\mathbf{U}_{1}\mathbf{S}_{1}\mathbf{V}_{1}^{\top}{}$, and noting that as $\tilde{\mathbf{G}}^{\prime}$ has full row rank $\mathbf{U}_{1}^{-1}{}=\mathbf{U}_{1}^{\top}{}$, the estimator matrix simplifies to

[TABLE]

This can be used to directly define the mean of the prediction posterior

[TABLE]

With similar substitutions as above, and using the Woodbury identity (seber_matrix_2008, , 15.3(b)(i)), the covariance of the prediction posterior simplifies to

[TABLE]

where $\mathbf{W}=\mathbf{V}_{2}\mathbf{V}_{2}^{\top}{}=\mathbf{I}-\mathbf{V}_{1}\mathbf{V}_{1}^{\top}{}$ represents a projection on to the nullspace of $\tilde{\mathbf{G}}^{\prime}$. This concludes the proof of part (1) of the proposition.

For the scheme to be conservative, the mean and covariance of the posterior must obey

[TABLE]

Also, the expectation over $\mathbf{d},\mathbf{p}$ is equivalent to an expectation over the independent variables $\mathbf{x},\boldsymbol{\delta},\boldsymbol{\rho}$, where $\mathbf{d}=\mathbf{G}\mathbf{x}+\boldsymbol{\delta}$ and $\mathbf{p}=\mathbf{Y}\mathbf{x}+\boldsymbol{\rho}$.

Now, the difference between the prediction $\mathbf{p}=\mathbf{Y}\mathbf{x}+\boldsymbol{\rho}$ and the mean $\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{d}$ simplifies to

[TABLE]

Furthermore, as $\mathbf{v}=\mathbf{C}^{\dagger}{}\mathbf{x}$, $\boldsymbol{\delta}$ and $\boldsymbol{\rho}$ are all independent and zero mean, then the difference $\mathbf{p}-\hat{\boldsymbol{\mu}}_{\mathbf{p}|\mathbf{d}}^{d}$ has an expected value of zero and a covariance given which is equivalent to $\hat{\mathbf{\Sigma}}_{\mathbf{p}|\mathbf{d}}^{d}$

[TABLE]

Thus, the required condition Eq. 72 is satisfied (with equality) and the prediction scheme is conservative. This proves part (2) of the proposition.