Adaptive Anchored Inversion for Gaussian
Random Fields Using Nonlinear Data

We begin with an initial approximation to the posterior $p(\theta\,|\,z,\bm{\mu},\bm{\Sigma})$, denoted by $f^{(0)}(\theta)$, which is taken to be a sufficiently diffuse normal distribution. The first iteration updates $f^{(0)}$ to $f^{(1)}$. In general, the $k$th iteration updates the approximation $f^{(k-1)}$ to $f^{(k)}$ as follows.

1. 1.

   Drawing a random sample of the parameters.

   Take $n$ random samples, denoted by $\theta_{1},\dotsc,\theta_{n}$, from $f^{(k-1)}$. For $i=1,\dotsc,n$, compute the prior density $s_{i}=\pi(\theta_{i})$ according to (7), and the “proposal density” $t_{i}=f^{(k-1)}(\theta_{i})$; let $w_{i}=\frac{s_{i}/t_{i}}{\sum_{j=1}^{n}s_{j}/t_{j}}$.

2. 2.

   Evaluating the forward model.

   For $i=1,\dotsc,n$, sample $\bm{y}_{i}$ from $p(\bm{y}\,|\,\theta_{i},\bm{\mu},\bm{\Sigma})$ according to the density (6), and evaluate $z_{i}=\mathcal{M}(\bm{y}_{i})$. The sampling of $\bm{y}_{i}$ is typically accomplished by a geostatistical simulation algorithm outlined in section 3, while evaluation of $z_{i}$ usually entails running a numerical forward model with the simulated field $\bm{y}_{i}$ as input.

3. 3.

   Reducing the dimension of the conditioning data by principal component analysis (PCA).

   Let $\bm{Z}={\left[\vbox{\hbox{\kern-1.15198pt\raise 0.0pt\hbox{\kern 46.08006pt}\kern-46.08006pt\vbox{\offinterlineskip\vbox{}}\kern-1.15198pt}}\right]}^{T}$, where $z_{1},\dotsc,z_{n}$ are column vectors of length $n_{z}$. By PCA, we find a $n_{z}\times m$ matrix $\bm{A}$, where $m\leq n_{z}$, such that the matrix $\bm{Z}\bm{A}$ “explains” specified proportions (e.g. 99%) of the variations in $\bm{Z}$. Subsequently, we transform each simulated $z_{i}$, $i=1,\dotsc,n$, to ${\bm{A}}^{T}z_{i}$, and also transform the observation vector $z$ similarly. This transformation reduces the dimension of $Z$ from $n_{z}$ to $m$. To keep the notation simple, we shall continue to use the symbol $z$ for this transformed variable.

4. 4.

   Estimating the joint density of $(\text{\tt parameter, data})$.

   Based on the weighted sample $\{(\theta_{1},z_{1}),\dotsc,(\theta_{n},z_{n})\}$ with weights $w_{i}$, approximate the density function of the random vector $(\theta,Z)$ by a multivariate normal mixture:

   $$p(\theta,z)\approx\sum_{i=1}^{n}w_{i}\,p_{i}(\theta,z)=\sum_{i=1}^{n}w_{i}\,\phi(\theta,z\,|\,\bm{m}_{i},\bm{V}_{i}),$$

   (9)

   where $\bm{m}_{i}$ and $\bm{V}_{i}$ are the mean vector and covariance matrix of the $i$th mixture component.

   We use a normal-kernel density estimator with two tuning parameters: one is bandwidth, which is a standard tuning parameter; the other is a localization parameter, which specifies how large a fraction of the entire sample is used to calculate the covariance matrix of each mixture component. These two tuning parameters are determined simultaneously by a likelihood-based optimization procedure. Because this step is a completely modular component in the entire algorithm (and alternative normal-kernel density estimators may be used without any change to other parts of the algorithm), we refer the reader to Zhang [2012] for details.

5. 5.

   Conditioning on the data.

   For $i=1,\dotsc,n$, partition $\bm{V}_{i}$ as $\left[\vbox{\hbox{\kern-1.15198pt\raise 0.0pt\hbox{\kern 51.84009pt}\kern-51.84009pt\vbox{\offinterlineskip\vbox{}}\kern-1.15198pt}}\right],$ where $\bm{V}_{\theta\theta i}$ is the covariance matrix of $\theta$, $\bm{V}_{zzi}$ is the covariance matrix of $Z$, and $\bm{V}_{\theta zi}$ as well as $\bm{V}_{z\theta i}$ are cross-covariance matrices. Similarly, decompose $\bm{m}_{i}$ into $\bm{m}_{\theta i}$ and $\bm{m}_{zi}$, representing the $\theta$ part and $Z$ part, respectively. Use the notations $p_{i}(\theta)$, $p_{i}(z)$, $p_{i}(\theta\,|\,z)$, and $p_{i}(z\,|\,\theta)$ as defined for relation (8). Then

   $$p_{i}(z)=\phi(z\,|\,\bm{m}_{zi},\bm{V}_{zzi})$$

   (10)

   and

   $$p_{i}(\theta\,|\,z)=\phi(\theta\,|\,\bm{m}_{\theta i}^{\prime},\bm{V}_{\theta\theta i}^{\prime}),$$

   (11)

   where

   $$\bm{m}_{\theta i}^{\prime}=\bm{m}_{\theta i}+\bm{V}_{\theta zi}\bm{V}_{zzi}^{-1}(z-\bm{m}_{zi}),\quad\bm{V}_{\theta\theta i}^{\prime}=\bm{V}_{\theta\theta i}-\bm{V}_{\theta zi}\bm{V}_{zzi}^{-1}\bm{V}_{z\theta i}.$$

   (12)

   Here $z$ is the observed forward data. Substituting (10)-(11) into (8), we get the conditional density

   $$p(\theta\,|\,z,\bm{\mu},\bm{\Sigma})\approx\sum_{i=1}^{n}v_{i}\,\phi(\theta\,|\,\bm{m}_{\theta i}^{\prime},\bm{V}_{\theta\theta i}^{\prime}),\quad\text{where }v_{i}=\frac{w_{i}\,p_{i}(z)}{\sum_{j=1}^{n}w_{j}\,p_{j}(z)}.$$

   (13)

   The posterior approximation $f^{(k-1)}$ is now updated to $f^{(k)}(\theta)=p(\theta\,|\,z)$. This completes an iteration of the algorithm.

Note $f^{(k)}$ is a normal mixture, just like $f^{(k-1)}$, and is ready to be updated in the next iteration.

Suppose the combined model and measurement error in $z$ is additive, as represented in (1). Assume $\epsilon$ is independent of $\theta$, $\bm{Y}$, and $Z$. Further assume the error is normal, $\epsilon\sim N(\bm{0},\bm{\Sigma}_{\epsilon})$. Then the $\bm{V}_{zzi}$ in (10) and (12) should be replaced by $\bm{V}_{zzi}+\bm{\Sigma}_{\epsilon}$.

The algorithm is able to accommodate more general errors. For example, non-additive errors and errors that might occur in certain internal steps of the forward model $\mathcal{M}$ can be built into the forward evaluation $z_{i}=\mathcal{M}(\bm{y}_{i})$ of step 2 by introducing randomness according to a description of the errors. This way, multiple forms of errors can be accommodated. An overall, explicit formula for the errors as a whole is not required. The passage from $\theta$ to $z$ in steps 1–2 embodies our best knowledge of the parameter-data connection: part in the anchor parameterization $p(\bm{y}\,|\,\theta,\bm{\mu},\bm{\Sigma})$, part in the forward model $\mathcal{M}(\bm{y})$, and part in errors in our implementation of $\mathcal{M}$ and in the measurement of $z$.

Going one step further, it is possible to relax the requirement that the forward model $\mathcal{M}$ be “deterministic”. If the forward model is stochastic [see an example in Besag et al., 1991], the parameter-data connection, as depicted in (4), includes one more layer of randomness. As far as the algorithm is concerned, the situation is not different from inclusion of model errors.

The overall idea in this algorithm may be summarized as “simulate $p(\theta,Z)$, then condition on the observed $z$.”

The algorithm does not require one to be able to draw a random sample from the prior of $\theta$. Instead, a convenient initial approximation $f^{(0)}$ starts the procedure. One only needs to be able to evaluate the prior density at any particular value of $\theta$. This provides flexibilities in specifying both the prior $\pi(\theta)$ and the initial approximation $f^{(0)}$.

Some statistics of the distribution $f^{(k)}$ may be examined semi-analytically, taking advantage of its being a normal mixture. More often, one is interested in the unknown field $\bm{Y}$ or a function thereof rather than in the parameter $\theta$. These may be investigated based on a sample of $\bm{Y}$ based on the posterior of $\theta$, as explained in section 3.

Typically, by far the most expensive operation in this algorithm is evaluating the forward model $\mathcal{M}(\bm{y})$. This occurs once for each sampled value of $\theta$. In contrast, sampling $\theta$ or generating field realizations are easy, because the distributions involved are normal or normal mixtures.

The versatility of normal mixture in approximating complex densities is well documented [Marron and Wand, 1992; West, 1993; McLachlan, 2000]. This approximation requires every component of $(\theta,Z)$ to be a continuous variable defined on $(-\infty,\infty)$, hence some transformations may be necessary in the definition of $\theta$ and $Z$. (See section 7 for examples.)

The dimension reduction achieved by principal component analysis in step 3 can be significant, for example with time series data. This facilitates the use of high-dimensional (i.e. large $n_{z}$) data without worrying (too much) about the correlation between the data components. The dimension reduction is also a stability feature. For example, components of $Z$ that are almost constant in the simulations will not cause trouble.

The global structure of this algorithm bears similarities to that of West [1993], with important differences. While West [1993] focuses on situations where the likelihood function is known, the algorithm here completely forgoes the requirement of a (or the) likelihood function. This distinction has far-reaching implications for the applicability of the proposed method.

This is an important feature of the proposed algorithm: it does not require a likelihood function. This feature distinguishes itself from the Markov chain Monte Carlo (MCMC) methods. There is a variant of MCMC called “approximate Bayesian computation (ABC)” [see Marjoram et al., 2003] that applies MCMC in problems without known likelihood functions. ABC uses a measure of the simulation-observation mismatch in lieu of likelihood. For multivariate data, there are significant difficulties in how to define this measure. The algorithm described here is very different from ABC.

In fact, this algorithm is not intrinsically tied to the anchor parameterization or anchored inversion. It is a general algorithm for approximating the posterior distribution where the likelihood is unknown but the model (which corresponds to the likelihood) can be simulated.

The choice of anchorset is a model selection problem, because any anchorset represents a particular parameterization, or “model”, for the field $Y(x)$. To facilitate comparison of models, we measure the performance of a model by a kind of integrated predictive likelihood defined as

where $f^{(k)}(\theta)$ is the approximate posterior density of $\theta$ obtained in the $k$th iteration. This measure is analogous to Bayes factors [Kass and Raftery, 1995; George, 2006] or, in particular, “posterior” Bayes factors [Aitkin, 1991]. However, the current setting differs from common (Bayesian) model selection in several ways. For example, the likelihood $p(z\,|\,\theta,\bm{\mu},\bm{\Sigma})$ is unknown, and is being estimated; the ultimate subject is not really the model parameter $\theta$, but rather the field $Y$; the algorithm is iterative, which makes use of the data $z$ in each iteration and corrects for the repeated use of data by importance weighting.

As an approximation, we assume the predictive distribution, $\int p(z\,|\,\theta,\bm{\mu},\bm{\Sigma})\,f^{(k)}(\theta)\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta$, is normal. The predictions $z_{1},\dotsc,z_{n}$ generated in step 2 of the algorithm is a random sample of this distribution. Based on the sample mean $\overline{z}$ and sample covariance matrix $\bm{S}_{z}$ of $z_{1},\dotsc,z_{n}$, an approximation of $\log L_{I}$ could be computed by

where $z$ is the observation.

The measure $L$ is subject to a fair amount of uncertainty arising from several sources. First, it is computed based on random samples $\theta_{1},\dotsc,\theta_{n}$ and $z_{1},\dotsc,z_{n}$. Second, each $z$ is a random draw conditional on its corresponding $\theta$ (via a random draw of $\bm{y}$). Third, the normal assumption for $Z$ may not be fully justifiable. Fourth, the sample size ($n$) is rather moderate in view of the dimensions of the random variables ($n_{\theta}$ and $n_{z}$). In experiments we used $n$ in the thousands and $n_{\theta}+n_{z}$ in the tens or hundreds. In fact, when $n_{\theta}+n_{z}$ is relatively large (say about 100), it is not rare that the empirical covariance $\bm{S}_{z}$ is not invertible. In view of these subtleties that harm the robustness of $L$, we further simplify and take

where $[j]$ indicates the $j$th dimension of $Z$. This approximation ignores correlations between the components of $Z$.

Denote the anchorset at the beginning of the $k$th iteration by $\bm{H}$. Towards the end of this iteration (somewhere in step 5), we consider whether to switch to an alternative anchorset, say $\bm{H}_{*}$, in the next iteration. The decision to switch or not will be made after comparing the predicted performances of the model in the next iteration using anchors $\bm{H}$ and $\bm{H}_{*}$, respectively. Let the two performances be denoted by $L_{k+1}^{*}(\bm{H})$ and $L_{k+1}^{*}(\bm{H}_{*})$, respectively.

The model (5), the prior (7), and the relation (6) suggest that $\theta_{*}$ conditional on $\theta$ has a normal distribution, which we denote by $\pi(\theta_{*}\,|\,\theta,\bm{\mu},\bm{\Sigma})$ to stress the fact that it is derived from the “prior” relations. Analogously, we write $\pi(\theta\,|\,\theta_{*},\bm{\mu},\bm{\Sigma})$ and $\pi(\theta,\theta_{*}\,|\,\bm{\mu},\bm{\Sigma})$.

The distribution $f^{(k-1)}(\theta)$ induces a distribution for $\theta_{*}$, which we denote by $f_{*}^{(k-1)}(\theta_{*})$. In fact, the sampling in steps 1–2 of the $k$th iteration also generates a sample of the alternative anchors, $\{\theta_{*i}=\bm{H}_{*}\bm{y}_{i}\}_{i=1}^{n}$, and the sample is a random sample of $f_{*}^{(k-1)}(\theta_{*})$ by the definition of $f_{*}^{(k-1)}(\theta_{*})$. Alternatively, we can consider the sample $\{\theta_{*i}\}$ to be obtained via $f_{*}^{(k-1)}(\theta_{*})=\int_{\Theta}\pi(\theta_{*}\,|\,\theta,\bm{\mu},\bm{\Sigma})\,f^{(k-1)}(\theta)\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta$.

Suppose we continue to use $\bm{H}$ in the $(k+1)$th iteration. In step 5 of the $k$th iteration, $f^{(k-1)}(\theta)$ is updated to $f^{(k)}(\theta)$. The distribution $f^{(k)}(\theta)$ will be used in steps 1–2 of the $(k+1)$th iteration to generate a random sample of $Z$, which is then used to calculate the performance indicator $L_{k+1}^{*}(\bm{H})$. Due to the high computational cost of the forward model, however, we do not want to actually generate this sample of $Z$ before committing to using the anchors $\bm{H}$ in the $(k+1)$th iteration. The idea then is to use the sample $z_{1},\dotsc,z_{n}$ created in step 2 of the $k$th iteration, which is based on $f^{(k-1)}(\theta)$, as an importance sample from the distribution of $Z$ that would be generated in step 2 of the $(k+1)$th iteration, which would be based on $f^{(k)}(\theta)$. The importance weights are proportional to $f^{(k)}(\theta_{i})/f^{(k-1)}(\theta_{i})$, $i=1,\dotsc,n$. This weighted sample of $Z$ then provides weighted sample means $\overline{z[j]}$ and weighted sample variances $s^{2}_{z[j]}$ to be used in (14) for predicting $L_{k+1}^{*}(\bm{H})$.

Now suppose we switch to $\bm{H}_{*}$ in the $(k+1)$th iteration. First, at the end of step 5 of the $k$th iteration, we need to obtain $f_{*}^{(k)}(\theta_{*})$ instead of $f^{(k)}(\theta)$. Second, the distribution $f_{*}^{(k)}(\theta_{*})$ would give rise to a sample of $Z$ in step 2 of the $(k+1)$th iteration, but we again would like to use the sample $z_{1},\dotsc,z_{n}$ created in step 2 of the $k$th iteration as a weighted replacement. We address these two steps now.

1. (1)

   We would like to view the sample $\{\theta_{*i}\}$ as an importance sample from $\pi(\theta_{*})$, then we would be able to conduct steps 4–5 for $\theta_{*}$ (instead of $\theta$) and obtain $f_{*}^{(k)}(\theta_{*})$. The question is: What are the importance weights?

   If the anchorsets $\bm{H}$ and $\bm{H}_{*}$ have no common anchors, then

   $$\pi(\theta_{*})=\int_{\Theta}\pi(\theta_{*}\,|\,\theta,\bm{\mu},\bm{\Sigma})\,\pi(\theta)\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta=\int_{\Theta}\pi(\theta_{*}\,|\,\theta,\bm{\mu},\bm{\Sigma})\,\frac{\pi(\theta)}{f^{(k-1)}(\theta)}f^{(k-1)}(\theta)\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta.$$

   (15)

   This suggests that the importance weights are proportional to $\frac{\pi(\theta)}{f^{(k-1)}(\theta)}$.

   If $\bm{H}$ and $\bm{H}_{*}$ share some common anchors, then write $\theta=(\theta_{\square},\theta_{\diamond})$ and $\theta_{*}=(\theta_{\square},\theta_{\bigstar})$. Let $f^{(k-1)}(\theta_{\square})$ and $f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})$ denote the marginal and conditional distributions induced by $f^{(k-1)}(\theta)$. We have

   $$\begin{split}\pi(\theta_{*})&=\pi(\theta_{\square})\,\pi(\theta_{\bigstar}\,|\,\theta_{\square},\bm{\mu},\bm{\Sigma})\\ &=f^{(k-1)}(\theta_{\square})\frac{\pi(\theta_{\square})}{f^{(k-1)}(\theta_{\square})}\int_{\Theta_{\diamond}}\pi(\theta_{\bigstar}\,|\,\theta_{\diamond},\theta_{\square},\bm{\mu},\bm{\Sigma})\,\frac{\pi(\theta_{\diamond}\,|\,\theta_{\square},\bm{\mu},\bm{\Sigma})}{f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})}f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta_{\diamond}.\end{split}$$

   (16)

   Noticing that the sample $\{\theta_{*i}\}$ is obtained via $f^{(k-1)}(\theta_{\square})\,f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})\,\pi(\theta_{\bigstar}\,|\,\theta_{\diamond},\theta_{\square},\bm{\mu},\bm{\Sigma})$, the above suggests the importance weights are proportional to $\frac{\pi(\theta_{\square})}{f^{(k-1)}(\theta_{\square})}\frac{\pi(\theta_{\diamond}\,|\,\theta_{\square},\bm{\mu},\bm{\Sigma})}{f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})}=\frac{\pi(\theta)}{f^{(k-1)}(\theta)}$.

   Relations (15)–(16) show that the weights are proportional to $\pi(\theta_{i})/f^{(k-1}(\theta_{i})$ when $\{\theta_{*i}\}$ are taken as an importance sample from $\pi(\theta_{*})$. These are the same weights (i.e. the $w_{i}$s) had we continued to use $\theta$ in steps 4–5. In other words, the weights do not depend on the specific choice of $\bm{H}_{*}$. Using these weights and the sample $\{\theta_{*i},z_{i}\}$, steps 4–5 of the $k$th iteration produce $f_{*}^{(k)}(\theta_{*})$.

2. (2)

   After deriving $f_{*}^{(k)}(\theta_{*})$, we would view $\{z_{i}\}$, obtained in step 2 of the $k$th iteration, as an importance sample of the $Z$ that would be generated starting with $f_{*}^{(k)}(\theta_{*})$ in steps 1–2 of the $(k+1)$th iteration. The importance weights are equal to those of $\{\theta_{*i}\}$ obtained in the $k$th iteration as an importance sample from $f_{*}^{(k)}(\theta_{*})$.

   If the anchorsets $\bm{H}$ and $\bm{H}_{*}$ have no common anchors, then

   $$\begin{split}f_{*}^{(k)}(\theta_{*})&=\frac{f_{*}^{(k)}(\theta_{*})}{\pi(\theta_{*})}\int_{\Theta}\pi(\theta_{*}\,|\,\theta,\bm{\mu},\bm{\Sigma})\,\pi(\theta)\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta\\ &=\frac{f_{*}^{(k)}(\theta_{*})}{\pi(\theta_{*})}\int_{\Theta}\pi(\theta_{*}\,|\,\theta,\bm{\mu},\bm{\Sigma})\frac{\pi(\theta)}{f^{(k-1)}(\theta)}f^{(k-1)}(\theta)\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta.\end{split}$$

   (17)

   This suggests the importance weights are proportional to

   $$\frac{f_{*}^{(k)}(\theta_{*})/\pi(\theta_{*})}{f^{(k-1)}(\theta)/\pi(\theta)}.$$

   (18)

   If $\bm{H}$ and $\bm{H}_{*}$ share some common anchors, then, using the notation introduced above, we have

   $$\begin{split}f_{*}^{(k)}(\theta_{*})&=f_{*}^{(k)}(\theta_{\square})\,f_{*}^{(k)}(\theta_{\bigstar}\,|\,\theta_{\square})\\ &=f^{(k-1)}(\theta_{\square})\frac{f_{*}^{(k)}(\theta_{\square})}{f^{(k-1)}(\theta_{\square})}\frac{f_{*}^{(k)}(\theta_{\bigstar}\,|\,\theta_{\square})}{\pi(\theta_{\bigstar}\,|\,\theta_{\square})}\int_{\Theta_{\diamond}}\pi(\theta_{\bigstar}\,|\,\theta_{\diamond},\theta_{\square},\bm{\mu},\bm{\Sigma})\,\pi(\theta_{\diamond}\,|\,\theta_{\square})\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta_{\diamond}\\ &=f^{(k-1)}(\theta_{\square})\frac{f_{*}^{(k)}(\theta_{\square})}{f^{(k-1)}(\theta_{\square})}\frac{f_{*}^{(k)}(\theta_{\bigstar}\,|\,\theta_{\square})}{\pi(\theta_{\bigstar}\,|\,\theta_{\square})}\\ &\phantom{\mathrel{=}}\;\;\int_{\Theta_{\diamond}}\pi(\theta_{\bigstar}\,|\,\theta_{\diamond},\theta_{\square},\bm{\mu},\bm{\Sigma})\,\frac{\pi(\theta_{\diamond}\,|\,\theta_{\square})}{f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})}f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})\mathop{\mathrm{\mathstrut d}}\nolimits\!\theta_{\diamond},\end{split}$$

   (19)

   suggesting the importance weights are proportional to

   $$\frac{f_{*}^{(k)}(\theta_{\square})}{f^{(k-1)}(\theta_{\square})}\frac{f_{*}^{(k)}(\theta_{\bigstar}\,|\,\theta_{\square})}{\pi(\theta_{\bigstar}\,|\,\theta_{\square})}\frac{\pi(\theta_{\diamond}\,|\,\theta_{\square})}{f^{(k-1)}(\theta_{\diamond}\,|\,\theta_{\square})}=\frac{f_{*}^{(k)}(\theta_{*})/\pi(\theta_{\bigstar}\,|\,\theta_{\square})}{f^{(k-1)}(\theta)/\pi(\theta_{\diamond}\,|\,\theta_{\square})}=\frac{f_{*}^{(k)}(\theta_{*})/\pi(\theta_{*})}{f^{(k-1)}(\theta)/\pi(\theta)}.$$

   (20)

   The weights turn $\{z_{i}\}$ into a weighted sample in computing $\overline{z[j]}$ and $s^{2}_{z[j]}$ for predicting $L^{*}(\bm{H}_{*})$ by (14).

Next, if it turns out that $L^{*}(\bm{H}_{*})>L^{*}(\bm{H})$, we take $f_{*}^{(k)}(\theta_{*})$ as the output of the $k$th iteration, with the anchorset $\bm{H}_{*}$. Otherwise, the output is $f^{(k)}(\theta)$ with the anchorset $\bm{H}$.

A remaining question is, how do we identify, in a way that is flexible and efficient, candidate alternative anchorsets (i.e. $\bm{H}_{*}$’s) to examine?

Suppose the current anchorset, represented by $\bm{H}$, contains $n_{\theta}$ anchors. Recall the definition of anchorset: the model domain is partitioned into $n_{\theta}$ sub-regions; the mean value of $Y$ in each sub-region is designated an anchor. Let us call each sub-region the “support” of the corresponding anchor.

We consider $n_{\theta}$ alternative anchorsets, denoted by $\bm{H}_{1},\dotsc,\bm{H}_{n_{\theta}}$, each containing $n_{\theta}+1$ anchors. The anchorset $\bm{H}_{j}$, $j=1,\dotsc,n_{\theta}$, is identical to the original anchorset $\bm{H}$ except that the support of the $j$th anchor in $\bm{H}$ is split into two sub-regions, defining two new anchors.

The $2n_{\theta}$ newly-created “small” anchors resulting from splitting each anchor in $\bm{H}$ form an “umbrella” anchorset, denoted by $\bm{H}_{*}$, which defines anchor vector $\theta_{*}=\bm{H}_{*}\bm{Y}$. Similarly defined are anchor vectors $\theta=\bm{H}\bm{Y}$, $\theta_{1}=\bm{H}_{1}\bm{Y}$,…, $\theta_{n_{\theta}}=\bm{H}_{n_{\theta}}\bm{Y}$. It can be seen that each of $\theta$, $\theta_{1}$,…, $\theta_{n_{\theta}}$ is a linear function of $\theta_{*}$. In steps 3–5 of the algorithm, we estimate the density of $\theta_{*}$. (To estimate $f^{(k)}(\theta_{*})$, extract the sample $\{\bm{H}_{*}\bm{y}_{i}\}$, then use the weights $\{w_{i}\}$ in steps 3–5.) Because the estimate is a normal mixture, the densities of $\theta$, $\theta_{1}$,…, $\theta_{n_{\theta}}$ follow immediately. Based on these densities, we predict the $L^{*}$ for each of these $n_{\theta}+1$ anchorsets and choose the anchorset that achieves the largest $L^{*}$. Steps 3–5 of the algorithm are then repeated to estimate $f^{(k)}$ for the chosen anchorset.

By this strategy, steps 3–5 of the algorithm are performed twice in order to compare and choose from $n_{\theta}+1$ anchorsets, including the original and the alternatives. As iterations proceed, more alternative anchorsets are considered in each iteration as the size of (i.e. number of anchors contained in) the currently used anchorset grows. The computational cost, however, barely increases. This efficient strategy takes advantage of the fact that the weights $\{w_{i}\}$ remain the same for all anchorsets.

An obvious extension to this strategy allows an alternative anchorset to split more than one anchor (say up to two) in the original anchorset. This extension examines more alternative anchorsets without substantially increasing the computational cost, because the umbrella anchorset $\bm{H}_{*}$ remains unchanged.

As this adaptive procedure suggests, the inversion exercise begins with a low number of anchors and gradually increases the number of anchors in iterations. In particular, we may begin with an anchorset that bisects each dimension of the space; that is, use 2, 4, and 8 anchors in 1D, 2D, and 3D spaces, respectively. This completely relieves the user from the burden of specifying an initial anchorset. More anchors are introduced in the iterations as more computation has been invested in the task. New anchors tend to be introduced in regions of the model domain where a more detailed representation of the $Y$ field is predicted to bring about most improvement to the model’s capability in reproducing the observation $z$. The partition of the model domain by the anchorset is not uniform—it is tailored to the unknown field $Y$, the forward process $\mathcal{M}$, and existing anchors in an automatic, adaptive, and evolving fashion.

In steps 3–5, the original $\theta$ is replaced by the umbrella $\theta_{*}$ (obtained from the sample $\{\bm{y}_{i}\}$ via $\bm{H}_{*}\bm{y}_{i}$), which is usually twice as long as $\theta$. At the end of step 5, $p(\theta_{*}\,|\,z,\bm{\mu},\bm{\Sigma})$ is obtained.

Step 5 now does not conclude an iteration. Display (13) is followed by the following new step.

1. 6.

   Selecting an anchorset and deriving its distribution.

   From the estimated density $p(\theta_{*}\,|\,z,\bm{\mu},\bm{\Sigma})$, given by (13), derive the density of each alternative anchorset as well as the original anchorset, then predict $L^{*}$ for each of these anchorsets.

   Pick the anchorset that achieves the largest $L^{*}$. Use this anchorset to repeat steps 3–5 and derive its density $p(\theta\,|\,z,\bm{\mu},\bm{\Sigma})$. This density becomes $f^{(k)}(\theta)$, where $\theta$ corresponds to the picked anchorset.

   This concludes an iteration of the algorithm.

The mean vector $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$ of the field may be parameterized in conventional geostatistical forms. Because of the existence of anchors, sophisticated trend models, such as polynomials in spatial coordinates, are usually unnecessary. One may model the mean $\bm{\mu}$ by a global constant:

The stationary covariance may be parameterized as

where $\eta^{2}$ is variance and $\varphi$ contains parameters pertaining to range, smoothness, nugget, geometric anisotropy, and so on. We leave specifics of this parameterization open. See section 7.1 for an example of a particular form of this parameterization.

A requirement imposed by the proposed inversion method is that all unknown parameters are defined on supports that can be mapped onto $(-\infty,\infty)$. This usually means that the support of each parameter should be in one of the following forms: $(-\infty,\infty)$, $(-\infty,c)$, $(c,\infty)$, and $(c_{1},c_{2})$. Once this requirement is satisfied, all parameters are transformed to be defined on $(-\infty,\infty)$, hence the posterior distribution can be approximated by a normal mixture. Note the support of each anchor is already $(-\infty,\infty)$ because the support of $Y$ is $(-\infty,\infty)$.