Uncertainty Quantification and Experimental Design for Large-Scale Linear Inverse Problems under Gaussian Process Priors ^††thanks: \fundingThis work was funded by the Swiss National Science Fundation (SNF) through project no. 178858.

We now consider a sequential data assimilation problem where data is made available in discrete stages. At each stage, a set of observations described by a (discretized) operator $\underline{G}_{i}$ is made and one observes a realization $\bm{y}_{i}$ of

where $\bm{W}_{i}$ is some set of points in $D$. Then, the posterior mean and covariance after each stage of observation may be obtained by performing a low rank update of their counterparts at the previous stage. Indeed, (Travelletti and Ginsbourger,, 2022, Corollary 6) provides an extension of Chevalier et al., 2014c ; Emery, (2009); Gao et al., (1996); Barnes and Watson, (1992) to linear operator observations, and gives:

At each stage $i$, these formulae require computation of the $p_{i}\times m$ matrix $\lambda_{i}\left(\bm{X}\right)$, which involves a $p_{i}\times p_{i}$ matrix inversion, where $p_{i}$ is the dimension of the operator $\underline{G}_{n}$ describing the current dataset to be included. This allows computational savings by reusing already computed quantities, avoiding inverting the full dataset at each stage, which would require a $p_{tot}^{2}$ matrix inversion, where $p_{tot}=\sum_{i=1}^{n}p_{i}$.

In order for these update equations to bring computational savings, one has to be able to store the past covariances $K_{\bm{W}_{n}\bm{W}_{n}}^{(n-1)}$ (Chevalier et al.,, 2015). This makes their application to large-scale sequential Bayesian inverse problems difficult, since the covariance matrix on the full discretization may become too large for storage above a certain number of discretization points. The next section presents our main contributions to overcome this limitation. They rely on an implicit representation of the posterior covariance that allows the computational savings offered by the kriging update formulae to be brought to bear on large scale inverse problems.

We consider the same sequential data assimilation setup as in the previous section, and for the sake of simplicity we assume that $\bm{W}_{1},...,\bm{W}_{n}=\bm{X}$ and use the lighter notation $m^{(i)}:=m^{(i)}_{\bm{X}}$ and $K^{(i)}:=K^{(i)}_{\bm{X}\bm{X}}$. The setting we are interested in here is the one where $\bm{X}$ is so large that the covariance matrix gets bigger than the available computer memory.

Our key insight is that instead of building the full posterior covariance $K^{(n)}$ at each stage $n$, one can just maintain a routine that computes the product of the current posterior covariance with any other low rank matrix. More precisely, at each stage $n$, we provide a routine $\textrm{CovMul}_{n}$ (Algorithm 1), that allows to compute the product of the current covariance matrix with any thin matrix $A\in\mathbb{R}^{m\times q},~{}q\ll m$:

where thin is to be understood as small enough so that the result of the multiplication can fit in memory.

This representation of the posterior covariance was inspired by the covariance operator of Gaussian measures. Indeed, if we denote by $C_{\mu^{(n)}}$ the covariance operator of the Gaussian measure associated to the posterior distribution of the GP at stage $n$, then

Hence, the procedure $\textrm{CovMul}_{n}$ may be thought of as computing the action of the covariance operator of the Gaussian measure associated to the posterior on the Dirac delta functionals at the discretization points.

This motivates us to think in terms of an updatable covariance object, where the inclusion of new observations (the updating) amounts to redefining a right-multiplication routine. It turns out that by grouping terms appropriately in Theorem 2.1 such a routine may be defined by only storing low rank matrices at each data acquisition stage.  

Hence, in order to compute products with the posterior covariance at stage $n$, one only has to store $n$ matrices $\bar{K}_{i}$, each of size $m\times p_{i}$ and $n$ matrices $R_{i}^{-1}$ of size $p_{i}\times p_{i}$, where $p_{i}$ is the number of observations made at stage $i$ (i.e. the number of lines in $\underline{G}_{i}$). In turn, each of these objects is defined by multiplications with the covariance matrix at previous stages, so that one may recursively update the multiplication procedure $\textrm{CovMul}_{n}$. Algorithms 1, LABEL:, 2, and 3 may be used for multiplication with the current covariance matrix, update of the representation and update of the posterior mean.

To use Algorithm 1, one should be able to compute products with the prior covariance matrix $K^{(0)}$. This may be achieved by chunking the set of grid points into $n_{c}$ subsets $\bm{X}=\left(\bm{X}_{1},...,\bm{X}_{n_{c}}\right)$, where each $\bm{X}_{i}$ contains a subset of the points. Without loss of generality, we assume all subsets to have the same size $m_{c}$. We may then write the product as

Each of the subproducts may then be performed separately and the results gathered together at the end. The individual products then involve matrices of size $m_{c}\times m$ and $m\times q$. One can then choose the number of chunks so that these matrices can fit in memory. Each block $K^{(0)}_{\bm{X}_{i}\bm{X}}$ may be built on-demand provided $K^{(0)}_{\bm{X}\bm{X}}$ is defined through a given covariance function.
 
This ability of the prior covariance to be built quickly on-demand is key to our method. The fact that the prior covariance matrix does not need to be stored allows us to handle larger-than-memory posterior covariances by expressing products with it as a multiplication with the prior and a sum of multiplications with lower rank matrices.
 

For the sake of comparison, assume that all n datasets have the same size $p_{c}$ and let $p=np_{c}$ denote the total data size. The cost of computing products with the current posterior covariance matrix at some intermediate stage is given by:

Using this recursively, we can then compute the cost of creating the implicit representation of the posterior covariance matrix at stage $n$:

This can then be compared with a non-sequential approach where all datasets would be concatenated into a single dataset of dimension $p$. More precisely, define the $p\times m$ matrix $\underline{G}$ and the $p$-dimensional vector $\bm{y}$ as the concatenations of all the measurements and data vectors into a single operator, respectively vector. Then computing the posterior mean using Eq. 2 with those new observation operators and data vector the cost is, to leading order in $p$ and $m$:

In this light, we can now sum up the two main advantages of the proposed sequential approach:

• •

  the cubic cost $\mathcal{O}\left(p^{3}\right)$ arising from the inversion of the data covariances is decreased to $\mathcal{O}\left(p^{2}p_{c}\right)$ in the sequential approach

• •

  if a new set of observations has to be included, then the direct approach will require the $\mathcal{O}\left(m^{2}p\right)$ computation of the product $K\underline{G}^{T}$, which can become prohibitively expensive when the number of prediction points is large, whereas the sequential approach will only require a marginal computation of $\mathcal{O}\left(m^{2}p_{c}\right)$.

Aside from the computational cost, our implicit representation also provides significant memory savings compared to an explicit approach where the full posterior covariance matrix would be stored. The storage requirement for the implicit-representation as a function of the number of discretization points $m$ is shown in Figure 2.

To illustrate the various methods presented in Section 2.2, we here apply them to a toy two-dimensional example: we consider a GP discretized on a large square grid $\{0,...,M\}\times\{0,...,M\}$ and try to learn it through various types of linear operator data.
 
More precisely, we will here allow two types of observations: pointwise field values and Fourier coefficient data. Indeed, the field values at the nodes $Z_{kl},~{}k,l\in\{1,...,M\}$ are entirely determined by its discrete Fourier transform (DFT):

Note that the above may be viewed as the discretized version of some operator and thus fits our framework for sequential conditioning. One may then answer questions such as how much uncertainty is left after the first 10 Fourier coefficient $F_{kl},~{}k,l=1,...,10$ have been observed? Or which Fourier coefficient provide the most information about the GP?

For example, Fig. 3 compares the posterior mean after observation of various Fourier coefficients, compared to observing pointwise values along a space-filling sequence. The GP model is a constant mean Matérn $5/2$ random field on a square domain $\left[-1,1\right]^{2}$. The domain is discretized onto a $400\times 400$ grid. Note that, by nature, each Fourier coefficient involves the field values at all points of the discretization grid and thus direct computation of the posterior mean requires the full $400^{2}\times 400^{2}$ covariance matrix (which would translate to roughly 100 GB of storage). This makes this situation suitable to demonstrate the techniques presented in this section.

In this example, the Fourier coefficients are ordered by growing $l_{\infty}$ norm. One observes that Fourier coefficients provide very different information than pointwise observations and decrease uncertainty in a more spatially uniform way, as shown in Fig. 4.

The extent to which Fourier coefficients provide more or less information than pointwise observations (depending on ordering) depends on the correlation length of the field’s covariance function. Indeed, for long correlation lengths, the low frequency Fourier coefficients contain most of the information about the field’s global behaviour, whereas for small correlation length, high frequency coefficients are needed to learn the small-scale structure.

When using Gaussian process priors to solve inverse problems, one has to select the hyperparameters of the prior. There exists different approaches for optimizing hyperparameters. We here only consider maximum likelihood estimation (MLE).

We restrict ourselves to GP priors that have a constant prior mean $m_{0}\in\mathbb{R}$ and a covariance kernel $k$ that depends on a prior variance parameter $\sigma_{0}^{2}$ and other correlation parameters $\bm{\theta}_{0}\in\mathbb{R}^{t}$:

where $r(.,.;\bm{\theta}_{0})$ is a correlation function, such that $r(x,x;\bm{\theta}_{0})=1,\forall x\in D$. The maximum likelihood estimator for the hyperparameters may then be obtained by minimizing the negative marginal log likelihood (nmll) of the data, which in the discretized setting of Section 2 may be writen as (Rasmussen and Williams,, 2006):

Since only the quadratic term depends on $m_{0}$, we can adapt concentration identities (Park and Baek,, 2001) to write the optimal $m_{0}$ as a function of the other hyperparameters:

where $\bm{1}_{m}$ denotes the $m$-dimensional column vector containing only $1$’s. Here we always assume $R$ to be invertible. The remaining task is then to minimize the concentrated nmll:

Note that the main computational challenge in the minimization of Eq. 8 comes from the presence of the $m\times m$ matrix $K_{\bm{X}\bm{X}}$. In the following, we will only consider the case of kernels that depend on a single length scale parameter: $\bm{\theta}_{0}=\lambda_{0}\in\mathbb{R}$, though the procedure described below can in principle be adapted for multidimensional $\bm{\theta}_{0}$.

In practice, for kernels of the form Eq. 7 the prior variance $\sigma_{0}^{2}$ may be factored out of the covariance matrix (for known noise variance), so that only the prior length scale $\lambda_{0}$ appears in this large matrix. One then optimizes these parameters separately, using chunking (Section 2.3) to compute matrix products. Since $\sigma_{0}$ only appears in an $p\times p$ matrix which does not need to be chunked (the data size $p$ being moderate in real applications), one can use automatic differentiation libraries such as Paszke et al., (2019) to optimize it by gradient descent. On the other hand, there is no way to factor out $\lambda_{0}$ out of the large matrix $K_{\bm{X}\bm{X}}$, so we resort to a brute force approach by specifying a finite search space for it. To summarize, we proceed here in the following way:

1. (i)

   (brute force search) Discretize the search space for the length scale by only allowing $\lambda_{0}\in\Lambda_{0}$, where $\Lambda_{0}$ is a discrete set (usually equally spaced values on a reasonable search interval);

2. (ii)

   (gradient descent) For each possible value of $\lambda_{0}$, minimize the (concentrated) $\mathcal{L}$ over the remaining free parameter $\sigma_{0}$ by gradient descent.

We ran the above approach on the Stromboli dataset with standard stationary kernels (Matérn $3/2$, Matérn $5/2$, exponential). In agreement with Linde et al., (2014), the observational noise standard deviation is $0.1~{}[mGal]$. The optimization results for different values of the length scale parameter are shown in Fig. 6. The best estimates of the parameter values for each kernel are shown in Table 1. The table also shows the practical range $\bar{\lambda}$ which is defined as the distance at which the covariance falls to $5\%$ of its original value.

We asses the robustness of each kernel by predicting a set of left out observations using the other remaining observations. Fig. 7 displays RMSE and negative log predictive density for different proportion of train/test splits.

Note the the above procedure is more of a quality assurance than a rigorous statistical evaluation of the model, since all datapoints were already used in the fitting of the hyperparameters. Due to known pathologies of the exponential kernel (MLE for length scale parameter going to infinity), we choose to use the Matérn 3/2 model for the experiments of Section 3.2 and Section 3.3. The maximum likelihood estimator of the prior hyperparameters for this model are $\hat{m}_{0}^{MLE}=2139.1~{}[kg/m^{3}]$, $\hat{\sigma}_{0}^{MLE}=284.65~{}[kg/m^{3}]$ and $\hat{\lambda}_{0}^{MLE}=651.6~{}[m]$.

Our implicit representation also allows for efficient sampling from the posterior by using the residual kriging algorithm (Chilès and Delfiner,, 2012; de Fouquet,, 1994), which we here adapt to linear operator observations. Note that in order to sample a Gaussian process at $m$ sampling points, one needs to generate $m$ correlated Gaussian random variables, which involves covariance matrices of size $m^{2}$, leading to the same computational bottlenecks as described in Section 2. On the other hand, the residual kriging algorithm generates realizations from the posterior by updating realizations of the prior, as we explain next.

As before, suppose we have a GP $Z$ defined on some compact Euclidean domain $D$ and assume $Z$ has continuous sample paths almost surely. Furthermore, say we have $p$ observations described by linear operators $\ell_{1},...,\ell_{p}\in C(D)^{*}$. Then the conditional expectation of $Z$ conditional on the $\sigma$-algebra $\Sigma:=\sigma\left(\ell_{1}\left(Z\right),...,\ell_{p}\left(Z\right)\right)$ is an orthogonal projection (in the $L^{2}$-sense (Williams,, 1991)) of $Z$ onto $\Sigma$. This orthogonality can be used to decompose the conditional law of $Z$ conditional on $\Sigma$ into a conditional mean plus a residual. Indeed, if we let $Z^{{}^{\prime}}$ be another GP with the same distribution as $Z$ and let $\Sigma^{{}^{\prime}}:=\sigma\left(\ell_{1}\left(Z^{{}^{\prime}}\right),...,\ell_{p}\left(Z^{{}^{\prime}}\right)\right)$, then we have the following equality in distribution:

Compared to direct sampling of the posterior, the above approach involves two main operations: sampling from the prior and conditioning under operator data. When the covariance kernel is stationary and belongs to one of the usual families (Gaussian, Matérn), methods exist to sample from the prior on large grids (Mantoglou and Wilson,, 1982); whereas the conditioning part may be performed using our implicit representation.

As a last example of application where our implicit update method provides substantial savings, we consider a sequential data collection task involving an inverse problem. Though sequential design criterion for inverse problems have already been considered in the litterature (Attia et al.,, 2018), most of them only focus on selecting observations to improve the reconstruction of the unknown parameter field, or some linear functional thereof.

We here consider a different setting. In light of recent progress in excursion set estimation (Azzimonti et al.,, 2016; Chevalier et al.,, 2013), we instead focus on the task of recovering an excursion set of the unknown parameter field $\rho$, that is, we want to learn the unknown set $\Gamma^{*}:=\{x\in D:\rho\left(x\right)\geq T\}$, where $T$ is some threshold. In the present context of Stromboli, high density areas are related to dykes (previous feeding conduits of the volcano), while low density values are related to deposits formed by paroxysmal explosive phreato-magmatic events (Linde et al.,, 2014). To the best of our knowledge such sequential experimental design problems for excursion set learning in inverse problems have not been considered elsewhere in the litterature.

We here consider a sequential setting, where observations are made one at a time and at each stage we have to select which obsevation to make next in order to optimally reduce the uncertainty on our estimate of $\Gamma^{*}$. Building upon Picheny et al., (2010); Bect et al., 2012a ; Azzimonti et al., (2019); Chevalier et al., 2014a , there exists several families of criteria to select the next observations. Here, we restrict ourselves to a variant of the weighted IMSE criterion (Picheny et al.,, 2010). The investigation of other state-of-the-art criteria is left for future work. We note in passing that most Bayesian sequential design criteria involve posterior covariances and hence tend to become intractable for large-scale problems. Moreover in a sequential setting, fast updates of the posterior covariance are crucial. Those characteristics make the problem considered here particularly suited for the implicit update framework intoduced in Section 2.2.

The weighted IMSE criterion selects next observations by maximizing the variance reduction they will provide at each location, weighted by the probability for that location to belong to the excursion set $\Gamma^{*}$. Assuming that $n$ data collection stages have already been performed and using the notation of Section 2.1, the variant that we are considering here selects the next observation location by maximizing the weighted integrated variance reduction (wIVR):

where $s$ is some potential observation location, $K^{(n+1)}$ denotes the conditional covariance after including a gravimetric observation made at $s$ (this quantity is independent of the observed data) and $G_{s}$ is the forward operator corresponding to this observation. Also, here $p_{n}$ denotes the coverage function at stage $n$ (we refer the reader to Appendix B for more details on Bayesian set estimation). After discretization, applying Theorem 2.1 turns this criterion into:

where we have assumed that all measurements are affected by $\mathcal{N}(0,\tau^{2})$ distributed noise.

Note that for large-scale problems, the wIVR criterion in the form given in Eq. 13 becomes intractable for traditional methods because of the presence of the full posterior covariance matrix $K^{(n)}_{\bm{X}\bm{X}}$ in the parenthesis. The implicit representation presented in Section 2.2 can be used to overcome this difficulty. Indeed, the criterion can be evaluated using the posterior covariance multiplication routine Lemma 2.2 (where the small dimension $q$ is now equal to the number of candidate observations considered at a time, here $1$ but batch acquisition scenarios could also be tackled). New observations can be seamlessly integrated along the way by updating the representation using Algorithm 2.

Experiments and Results: We now study how the wIVR criterion can help to reduce the uncertainty on excursion sets within the Stromboli volcano. We here focus on recovering the volume of the excursion set instead of its precise location. To the best of our knowledge, in the existing literature such sequential design criteria for excursion set recovery have only been applied to small-scale inverse problems and have not been scaled to larger, more realistic problems where the dimensions at play prevent direct access to the posterior covariance matrix.

In the following experiments, we use the Stromboli volcano inverse problem and work with a discretization into cubic cells of $50~{}[m]$ side length. We use a Matérn $3/2$ GP prior with hyperparameters trained on real data (Table 1) to generate semi-realistic ground truths for the experiments. We then simulate numerically the data collection process by computing the response that results from the considered ground truth and adding random observational noise. When computing sequential designs for excursion set estimation, the threshold that defines the excursion set can have a large impact on the accuracy of the estimate. Indeed, different thresholds will produce excursion sets of different sizes, which may be easier or harder to estimate depending on the set estimator used. For the present problem, Fig. 9 shows the distribution of the excursion volume under the considered prior for different excursion thresholds.

It turns out that the estimator used in our experiments (Vorob’ev expectation) behaves differently depending on the size of the excursion set to estimate. Indeed, the Vorob’ev expectation tends to produce a smoothed version of the true excursion set, which in our situation results in a higher fraction of false positives for larger sets. Thus, we consider two scenarios: a large scenario where the generated excursion sets have a mean size of $10\%$ of the total inversion volume and a small scenario where the excursion sets have a mean size of $5\%$ of the total inversion volume. One should note that those percentages are in broad accordance with the usual size of excursion sets that are of interest in geology. The chosen thresholds are $2500~{}[kg/m^{3}]$ for the large excursions and $2600~{}[kg/m^{3}]$ for the small ones.

The experiments are run on five different ground truths, which are samples from a Matérn $3/2$ GP prior (see previous paragraphs). The samples were chosen such that their excursion set for the large scenario have volumes that correspond to the $5\%$, $27.5\%$, $50\%$, $72.5\%$ and $95\%$ quantiles of the prior excursion volume distribution for the corresponding threshold. Fig. 10 shows the prior excursion volume distribution together with the volumes of the five different samples used for the experiments. Fig. 11 shows a profile of the excursion set (small scenario) for one of the five samples used in the experiments. The data collection location from the 2012 field campaign (Linde et al.,, 2014) are denoted by black dots. The island boundary is denoted by blue dots. Note that, for the sake of realism, in the experiments we only allow data collection at locations that are situated on the island (data acquired on a boat would have larger errors); meaning that parts of the excursion set that are outside of the island will be harder to recover.

Experiments are run by starting at a fixed starting point on the volcano surface, and then sequentially choosing the next observation locations on the volcano surface according to the wIVR criterion. Datapoints are collected one at a time. We here only consider myopic optimization, that is, at each stage, we select the next observation site $s_{n+1}$ according to:

where ties are broken arbitrarily. Here $\bm{S}_{c}$ is a set of candidates among which to pick the next observation location. In our experiments, we fix $\bm{S}_{c}$ to consist of all surface points within a ball of radius $150$ meters around the last observation location. Results are summarized in Figs. 12 and 13, which shows the evolution of the fraction of true positives and false positives as a function of the number of observations gathered.

We see that in the large scenario (Fig. 12) the wIVR criterion is able to correctly detect 70 to 80% of the excursion set (in volume) for each ground truth after 450 observations. For the small scenario (Fig. 13) the amount of true positives reached after 450 observations is similar, though two ground truths are harder to detect.

Note that in Figs. 12 and 13 the fraction of false negatives is expressed as a percentage of the volume of the complementary of the true excursion set $D\setminus\Gamma^{*}$. We see that the average percentage of false positives after $450$ observations tends to lie between $5$ and $15\%$, with smaller excursion sets yielding fewer false positives. While the Vorob’ev expectation is not designed to minimize the amount of false positives, there exists conservative set estimators (Azzimonti et al.,, 2019) that specialize on this task. We identify the extension of such estimators to inverse problems as a promising venue for new research.

In both figures we also plot the fraction of true positives and false positives that result from the data collection plan that was used in Linde et al., (2014). Here only the situation at the end of the data collection process is shown. We see that for some of the ground truths the wIVR criterion is able to outperform static designs by around 10%. Note that there are ground truths where it performs similarly to a static design. We believe this is due to the fact that for certain ground truths most of the information about the excursion set can be gathered by spreading the observations across the volcano, which is the case for the static design that also considers where it is practical and safe to measure.

Limiting Distribution: The dashed horizontal lines in Figs. 12 and 13 show the detection percentage that can be achieved using the limiting distribution. We define the limiting distribution as the posterior distribution one were to obtain if one had gathered data at all allowed locations (everywhere on the volcano surface). This distribution may be approximated by gathering data at all points of a given (fine grained) discretization of the surface. In general, this is hard to compute since it requires ingestion of a very large amount of data, but thanks to our implicit representation (Section 2.2) we can get access to this object, thereby, allowing new forms of uncertainty quantification.

In a sense, the limiting distribution represents the best we can hope for when covering the volcano with this type of measurements (gravimetric). It gives a measure of the residual uncertainty inherent to the type of observations used (gravimetric). Indeed, it is known that a given density field is not identifiable from gravimetric data alone (see Blakely, (1995) for example). Even if gravity data will never allow for a perfect reconstruction of the excursion set, we can use the limiting distribution to compare the performance of different sequential design criteria and strategies. It also provides a mean of quantifying the remaining uncertainty under the chosen class of models. A sensible performance metric is then the number of observations that a given criterion needs to approach the minimal level of residual uncertainty which is given by the limiting distributions.

As a last remark, we stress that the above results and the corresponding reconstruction qualities are tied to an estimator, in our case the Vorob’ev expectation. If one were to use another estimator for the excursion set, those results could change significantly.

Posterior Volume Distribution: Thanks to our extension of the residual kriging algorithm to inverse problems (see Section 3.2), we are able to sample from the posterior at the end of the data collection process. This opens new venues for uncertainty quantification in inverse problems. For example, we can use sampling to estimate the posterior distribution of the excursion volume and estimate the residual uncertainty on the size of the excursion set.

Fig. 15 shows the empirical posterior distribution of the excursion volume for each of the ground truths considered in the preceding experiments. When compared to the prior distribution, Fig. 10, one sees that the wIVR criterion is capable of significantly reducing the uncertainty on the excursion volume. This shows that though the location of the excursion set can only be recovered with limited accuracy, as shown in Figs. 12 and 13, the excursion volume can be estimated quite well. This is surprising given that the criterion used (wIVR) is a very crude one and was not designed for that task. On the other hand, there exist more refined criterion, like the so-called SUR strategies (sequential uncertainty reduction) (Chevalier et al., 2014b, ; Bect et al.,, 2019), among which some were specifically engineered to reduce the uncertainty on the excursion volume (Bect et al., 2012b, ). Even though those criterion are more computationally challenging than the wIVR one, expecially in the considered framework, we identify their application to large Bayesian inverse problem as a promising avenue for future research.