{'fraction_non_alphanumeric': 0.039631966678064096, 'fraction_numerical': 0.016179167573280282, 'mean_word_length': 4.399714669352131, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bayesian multiscale deep generative model for the solution of high-dimensional inverse problems

Estimation of spatially-varying parameters for computationally expensive forward models governed by partial differential equations is addressed. A novel multiscale Bayesian inference approach is introduced based on deep probabilistic generative models. Such generative models provide a flexible representation by inferring on each scale a low-dimensional latent encoding while allowing hierarchical parameter generation from coarse- to fine-scales. Combining the multiscale generative model with Markov Chain Monte Carlo (MCMC), inference across scales is achieved enabling us to efficiently obtain posterior parameter samples at various scales. The estimation of coarse-scale parameters using a low-dimensional latent embedding captures global and notable parameter features using an inexpensive but inaccurate solver. MCMC sampling of the fine-scale parameters is enabled by utilizing the posterior information in the immediate coarser-scale. In this way, the global features are identified in the coarse-scale with inference of low-dimensional variables and inexpensive forward computation, and the local features are refined and corrected in the fine-scale. The developed method is demonstrated with two types of permeability estimation for flow in heterogeneous media. One is a Gaussian random field (GRF) with uncertain length scales, and the other is channelized permeability with the two regions defined by different GRFs. The obtained results indicate that the method allows high-dimensional parameter estimation while exhibiting stability, efficiency and accuracy.


Introduction
Inverse problems are important but challenging in many fields like geophysics, medical imaging, groundwater flows, and other. They address the estimation of model parameters from partial and noisy observations [1]. Two approaches for addressing inverse problems are typically employed. The deterministic methods convert parameter identification to an optimization problem that involves minimizing the misfit between model predictions and observations. Since limited observations are insufficient to identify the underlying parameters, regularization methods [2,3,4] are used to address this ill-posed problem. On the other hand, Bayesian inference approaches play a fundamental role in inverse problems allowing us to quantify the uncertainty of the solution and providing natural regularization via prior knowledge [5]. They treat parameters as random variables to highlight the uncertainty in their estimation. The non-uniqueness of the solution is addressed by computing the posterior of the parameters rather than a single point estimate. Variational inference (VI) and Monte Carlo (MC) methods are two main approximation methods to deal with the computation of the intractable posterior distribution. VI [6,7,8,9] is easy to implement but limited by the family of variational distributions. Most Bayesian approaches emphasize Markov Chain Monte Carlo (MCMC) methods that aim to generate samples from the posterior distribution that subsequently are used to produce statistics of the quantities of interest.
MCMC methods have the appealing property that they are asymptotically exact. Thus, many previous works have studied the MCMC method or its variants for Bayesian inverse problems (BIPs). However, there are two main difficulties for these methods. First, the dimensionality of the spatially-varying parameters can be high (e.g. equal to the number of grid points) leading to the so called curse of dimensionality. Second, these sampling-based approaches require multiple evaluations of the forward model (likelihood evaluation). Each evaluation involves a full forward simulation, which is computationally prohibiting for many practical problems governed by partial differential equations (PDEs).
For the first problem identified above, given prior information, parameterization methods are often used to provide a low-dimensional embedding of the unknown spatially-varying parameter. The common method in BIPs is the truncated Karhunen-Loève expansion (KLE) for the estimation of Gaussian random fields (GRFs) [10,11], where inference is performed over a small number of expansion coefficients. With limitations and strong assumptions on the mean and covariance functions, the KLE cannot reflect the true prior information, and is not a good choice for fields with nontrivial correlation structure. To address this, sparse grid interpolation [12,13] and wavelet-based [14] methods have been proposed. However, such methods still have difficulties in the parameterization of complex parameters such as multi-modal or non-Gaussian random fields [15,16].
DGMs are much more flexible and scalable, where the prior information is naturally incorporated into the training data without strong assumptions. Once the DGM is trained, one can sample latent variables from a low-dimensional simple distribution (like a Gaussian), and then generate the spatially-varying parameter using the pre-trained neural network.
A potential remedy of the requirement of MCMC methods for multiple calls to the forward model solver is to build a surrogate forward model, such as polynomial chaos [24,25], Gaussian process [26,27], or deep neural networks [28,29,30]. However, the surrogate model often introduces epistemic uncertainty that will result in broadening of the posterior for parameter estimation [31]. Furthermore, it is still a difficult task to construct an accurate surrogate for forward models with high-dimensional input using limited data. To reduce the computational burden of the simulation, multiscale [32,33,34] and multi-fidelity methods [35,36] have been applied to accelerate the Bayesian computation without sacri-ficing accuracy. The two-stage MCMC [37] designed a preconditioned Metropolis-Hastings algorithm to improve the acceptance rate in the fine-scale model. Inspired by the multilevel Monte Carlo, Multilevel MCMC methods [38,39,40,41] are proposed for BIPs to accelerate the estimation of the posterior distribution. All these methods are indeed promising for BIPs by leveraging the advantages of the accurate fine-scale model and the efficiency of the coarse-scale model.
In this work, we propose a multiscale deep generative model (MDGM) exploiting the multiscale nature of the parameter of interest. This extends existing DGMs and allows us to generate parameters on various scales with different discretization/resolution. Since GANs are notorious on training stability and mode collapse, and flow-based models [42,43] require an identical-dimensional latent space to the parameter, we derive the MDGM based on VAE. In the MDGM, we design a specific latent space that includes two latent variables, a low-dimensional latent variable that controls global and salient features and a higherdimensional latent variable that defines local and detailed features. Utilizing the hierarchical representation of the parameter and latent spaces, the multiscale inference is performed in the low-dimensional latent space rather than the original parameter space. This allows us to explore the posterior of the parameter from coarse-to fine-scales with a significant computational saving. Once most of the salient features are identified in the coarse-scale using a computationally inexpensive coarse-solver, the fine-scale estimation requires only few fine-scale simulations to refine the coarse-scale parameter estimation.

Bayesian inverse problems
In this section, we introduce the inverse problems of interest and briefly discuss the limitations of standard Bayesian inference approaches to inverse problems. We consider a spatially-varying parameter x(s) usually represented as a random field x(s, ω), where s is spatial location in the domain S and ω is a random event in the sample space Ω. This random field is discretized by a random vector x ∈ R M using standard finite element or finite difference discretization approaches. In our inverse problem setting, x(s), will be considered as our primary quantity of interest.
Let us consider a physical system governed by PDEs in a given spatial domain. We assume x(s, ω) to be an input parameter (e.g. material property) of this model. Of interest to this work are distributed properties with multiscale features. The forward model concerning this physical system is usually considered as a function F : R M → R D , which maps the unknown parameters x to the observable output D obs ∈ R D with a measurement noise ξ ∈ R D : The inverse problem is to infer the unknown parameters x based on these noisy data D obs .
In the particular problem we will focus in Section 4, our goal is to estimate the permeability field in a porous media flow using pressure measurements.
Without prior information about the measurement system and/or model evaluation, we assume that ξ is a zero-mean Gaussian noise with covariance matrix Σ, i.e., ξ ∼ N (0, Σ).
Since often dim(R D ) dim(R M ), the inverse problem is highly ill-posed and identification of the parameter x is highly-sensitive to this noise. The Bayesian paradigm [5] provides a general and natural way to treat the unknown parameter x as random variable to highlight the uncertainty in the inference process. Given the observation data D obs , one calculates the posterior probability π(x|D obs ) via Bayes' formula as follows: where π(x) is the prior distribution, and L e (D obs |x) is the likelihood function which evaluates the discrepancy between the forward predictions and observations. For the assumed case of Gaussian noise, we can define the likelihood function as As the parameter x of interest is high-dimensional, the normalization constant in Eq. (2) involves computing a high-dimensional integral that is often an intractable process. Thus approximate inference for the posterior π(x|D obs ) is performed using the unnormalized density, i.e., π(x|D obs ) ∝ L e (D obs |x)π(x).

Multiscale inference with MDGM
Without a closed-form expression, the posterior distribution in Eq. (4) must be computed numerically. To this end, MCMC [10,44] or other approximation methods like Ensemble Kalman filter (EnKF) [16,30] are often employed. However, there are still two main difficulties for these methods. MCMC and EnKF implementations will often fail to directly approximate the posterior of the high-dimensional spatially-varying parameter x. Moreover, for complex parameters (e.g. channelized permeability), the prior information cannot be easily cast as an explicit probability distribution. However, one often has access to a historical [16,45,46] where x (i) can be seen as samples from the underlying prior distribution π(x). One could use X to approximate π(x) with its empirical measure.
However, in this work, we will use this dataset to approximate the prior distribution with a generative model as follows: where p(z) is a simple distribution (e.g. Gaussian) for the latent variable z ∈ R d , and p θ (x|z) is a generative model parameterized by θ (decoder). In a DGM like a VAE, one can choose a Gaussian distribution N (0, I) for p(z), and N (µ θ (z), σ 2 I) for p θ (x|z), where µ θ (z) is the output of the decoder neural network, and σ is a hyperparameter that does not depend on the latent variable z. It is common practice in the literature [15,16,47,48] to ignore the noise and approximate the density p θ (x|z) with the point estimate µ θ (z). Once the model is trained, we can thus define a mapping from the latent space to the original parameter space, i.e. x = µ θ (z).
The parameters θ of the generative model p θ (x|z) can be computed using the given training dataset X by minimizing the Kullback--Leibler (KL) divergence D KL (π(x)||p(x|θ)) [49], . This leads to the equivalent problem of maximizing the marginal loglikelihood: The above marginalization is potentially very difficult to compute involving an intractable integration. Using Expectation-Maximization is also intractable as that will require the posterior p θ (z|x) that is also computationally intractable. The marginal likelihood for the can be reformulated using a variational density q φ (z (i) |x (i) ) (encoder) parameterized by φ. The details of these calculations will be given in Section 3.3.
To approximate the posterior distribution in Eq. (4) using the MCMC method, we are interested to generate realizations sampled from the underlying prior distribution π(x). To sample realizations from π(x) using the generative model, one can sample z i from the simple and low-dimensional distribution p(z) and then obtain the realization x (i) using the decoder model µ θ (z). Alternatively, instead of approximating the posterior π(x|D obs ) in Eq. (4), one can instead evaluate the low-dimensional posterior p(z|D obs ) using the following unnormalized density: where p(z) is an explicit distribution, and the likelihood can be evaluated using the decoder model µ θ (z) and the forward model F(x). The evaluation of the posterior of the lowdimensional latent variable z using MCMC or EnKF is computationally tractable [15,16].
To further improve the efficiency of the inference process, we will introduce a multiscale version of the above highlighted generative model to perform inference in each scale l = 1, 2, . . . , L from the coarsest-scale (l = 1) to the desired finest-scale (l = L). This multiscale scheme based on the MDGM contains a hierarchical simple distribution p(z l ) at each scale l and a conditional distribution p θ l (x l |z l ) that can generate the spatially-varying parameter x in each scale. Correspondingly, one can assess the posterior p(z l |D obs ) using MCMC with p(z l ), µ θ l (z l ), and F l (x l ). The details of this multiscale model are given next.

Multiscale dataset
Our physical systems of interest are governed by a system of PDEs, and the spatiallyvarying property of interest is a material property appearing e.g. in the constitutive equations. The forward problem defines the well-posed solution of the PDEs (with some boundary conditions) given appropriate material properties. Such problems are often solved in a discretized fashion with finite element or finite difference or spectral approximations for different levels of discretization of the spatial domain S. In this work, we are interested in a hierarchical parameterization of the spatially-varying parameters x l with different spatial discretization or resolutions at each scale l. If the forward model is performed in the 2-D space, the parameter random fields x l at the l-th scale are treated as images, e.g.
where H l , W l denote the number of the pixels in the horizontal and vertical directions, respectively.
For notational convenience, we assume that the finest scale parameters x L represent our "true parameter model". The noisy observations D obs in our numerical studies are taken from this discretization level. In the inverse problem of interest, our task is to compute x L given a finite number of observations. For the solution of this inverse problem, prior knowledge can provide useful information for x L before any observations. As prior information for our model, we assume that we are given a dataset X = {x To obtain images for training the generative model at different discretization levels, we will need to obtain a multiscale training dataset. This can be accomplished by upscaling the fine-scale training dataset [50,51,52]. For example, with an upscaling (deterministic) operator U : M l and M l−1 are the dimensions of x l and x l−1 , respectively. The dataset {x The datasets {x . . , L in different scales are obtained by adopting recursively U in Eq. (8) starting with the finest-scale l = L. we assume x (i) l is sampled from l-th scale underlying prior distribution π l (x l ). The operator U used in this paper is deterministic, which leads to a one-to-one correspondence between the elements in {x . One can choose different operators U depending on the particular parameter of interest. In this paper, we employ the arithmetic average 1 : where n e denotes the number of elements in the fine-scale l corresponding to one element in the coarse-scale (l − 1). The value at the coarse-grid element e is the mean of the values in the spatially corresponding elements e i in the fine-scale. Spatial correspondence between two adjacent scales with 50% coarsening in each direction is illustrated in Fig. 1.
An example illustrating this deterministic upscaling for channelized permeability using Eq. (9) is given in Fig. 2. The coarse-scale image provides a blurry representation of the fine-scale image but overall its features are consistent with those of the fine-scale image. It can be noticed that the coarse-scale image manifests itself with a checkerboard pattern that 1 http://www.epgeology.com/static-modeling-f39/how-upscale-permeability-t6045.html misses a lots of local information. With the scales of interest pre-determined and the training dataset defined at each scale, we are ready to train the MDGM and perform inference on each scale using the proposed hierarchical multiscale framework.

Model specification
Given the training dataset {x (i) } N i=1 , we seek to learn a DGM that can approximate π(x) with Eq. (5). The DGM involves the original parameter space and the latent space and the mappings between these spaces. For any x sampled from the underlying distribution π(x), the corresponding z is sampled from the conditional distribution q φ (z|x), where q φ (z|x) is called the recognition model or probabilistic encoder model, and φ are its model parameters.
In the reverse direction, one can sample a z from a simple and low-dimensional distribution p(z), and obtain the corresponding x from the generative model or probabilistic decoder model p θ (x|z).
, select an appropriate distribution p(z) for the latent variable z and learn the models p θ (x|z) and q φ (z|x), respectively, such that π(x) can be approximated using Eq. (5), where θ and φ denote the parameters of the generative and recognition models, respectively.
For high-dimensional inversion tasks, direct inference in the fine-scale is prohibited due to the computational cost of the forward model. In addition, inference of the latent parameters z that lead to good estimates of x through a generative model requires a high-dimensional z. This will lead to long exploration costs for MCMC and requires a large number of forward model evaluations. To this end, given the multiscale dataset as discussed in Section 3.1, we propose a multiscale scheme for posterior estimation by introducing a hierarchy of generative models from coarse-to fine-scales. In this scheme, the coarse-scale generative models have low-dimensional latent spaces. This together with inexpensive forward model evaluations in the coarse-scales would allow MCMC to explore the posterior in coarse-scales with much reduced cost. The computational savings can be even higher if the latent representation on a given scale utilizes the latent information that was inferred in the immediately coarser-scale.
In our construct, the latent variables z l at level l of the hierarchy are given as z l = (z l−1 , z l ), where the latent variables z l−1 and z l are encoded from the coarse-scale parameter x l−1 and the fine-scale parameter x l , respectively. The latent variable z l−1 can generate x l−1 through the generative model at scale (l − 1). It also impacts the generation of x l in the fine-scale l by way of dominating its salient features since z l−1 captured the information from the immediately coarser-scale. In this setting with z l−1 encoded from x l−1 , it is anticipated that x l generated by p θ l (x l |z l−1 , z l ) would sustain most of the features of x l−1 . We extend the Definition 3.1 to a multiscale scenario as follows.
Once the MDGM is established, the parameters x l are encoded by the latent variable z l so that one can perform inference of z l . Given the prior distribution p(z l ), the decoder model µ θ (z l ), the forward model F l , and observations D obs , inference of the posterior of z l is performed as follows: Note that MCMC converges and captures prominent and valuable features quickly in the coarse-scale (l − 1) since z l−1 is low-dimensional and the forward model F l−1 is less expensive in comparison to F l . For an efficient Bayesian inference at each scale l, we are interested in using the posterior distribution at coarse-scale (l − 1) to provide an informative prior information or improve sampling efficiency in the next finer-scale. This avoids relying completely on inference in a high-dimensional latent space where direct computation of finescale details would increase the model complexity. The purpose of inference on fine-scale is to correct the details rather than run long exploration for capturing all appropriate features.
A related idea was implemented earlier using hierarchical structured sparse grids in [13].
In summary, the inverse problem is divided into a multiscale posterior estimation, with the inference of parameters proceeding from coarse-to fine-scale. The definition of the multiscale posterior estimation problem is given next.

Probabilistic generative model
The MDGM is used to generate parameters with different discretization/resolution. We construct such a model based on the variational autoenoder (VAE). The coarsest generative model that involves a single-scale (l = 1) is vanilla VAE. It employs the dataset {x (generated as discussed in Section 3.1) sampled from the underlying distribution π(x 1 ), i.e.
x (i) 1 i.i.d ∼ π(x 1 ). We consider below the probabilistic generative model on a single-scale before deriving the multiscale formulation in Section 3.4. For simplicity of the notation, we drop the subscript 1 in the equations below even though this model will be used in the scale l = 1.
Given the training dataset {x (i) } N i=1 , one can introduce a variational family q φ (z|x) to convert the intractable computation of maximizing the marginal log-likelihood in Eq. (6) into an optimization problem, where φ denotes the model parameters. It can be written as follows: where the last step is the application of Jensen's inequality. The above lower bound is called the variational lower bound. For a given training dataset, one can maximize the variational lower bound rather than the marginal log-likelihood. This is an optimization problem with respect to model parameter θ and φ. The variational lower bound in Eq. (11) can be written as follows, The minimization of the −L(θ, φ) balances the optimization of both the recognition and generative models. Thus the recognition model parameters φ are learned jointly with the generative model parameters θ [18]. The graphical model is shown in Fig. 3. Figure 3: The directed graphical model for the probabilistic model [18]. The latent variable z (i) of each configuration x (i) is obtained by the probabilistic recognition model q φ (z (i) |x (i) ). The variational approximation is indicated with dashed edges and the generative model p θ (x|z)p(z) with solid edges.
To evaluate L(θ, φ) in Eq. (12), there are three probability distributions to be identified. As mentioned in Definition 3.1, we shall select appropriate simple distributions for the latent variables z. For example, in this paper, we let p(z) ∼ N (0, I). The (12) is the expected log-likelihood. It encourages the reconstructed datax by the decoder to approximate the original data x. We assume where the mean µ(z) is the output of a decoder neural network and σ is a constant hyperparameter that does not depend on the decoder so that it can be ignored during optimization. LetΣ = σ 2 I, the first term in Eq. (12) then takes the form based on the minibatches: where n denotes the number of training samples of x (also referred to as the batch size in the training of deep neural networks). For each epoch, there are N n minibatches, each batch uniformly sampled from the dataset . m is the number of z samples from q φ (z|x) for an expectation approximation using the Monte Carlo method. One can also refer to Eq. (13) as the reconstruction error.
is taken as a Gaussian distribution, where the meanμ(x) and the varianceσ 2 (x) are outputs of the encoder network.
The KL-divergence can be analytically computed when both distributions are Gaussian [18].
The KL-divergence works as an objective function for the optimization problem with respect to the encoder parameters φ. Based on the minibatches, it can be written as: whereμ k (x) andσ k (x) denote the k-th element of the mean and standard deviation, respectively, which are outputs of the encoder network. The results of Eqs. (13) and (14) summarize the objective function for optimization of the encoder and decoder neural networks.
Inspired from β-VAE [53], we slightly modify the underlying loss function. This modification of the VAE results by adding an extra hyperparameter β to the KL divergence.
This hyperparamter can constrict the capacity of the latent bottleneck and encourage a disentangled representation. We would like the individual dimensions of the latent variable z to be interpretable or to correspond to some features of parameters thus disentangling the true variation of data [53,54]. A good interpretability and factorized representation for the latent variable will facilitate the feature exploration in MCMC. Note that [55] has shown that β-VAE is an important method in disentangled representation for learning the compositional and hierarchical visual concept. Choosing an appropriate hyperparameterβ for loss function is important in MDGM. We can define the loss functionL ≡ − n N L(θ, φ) in each training iteration, the loss function for the single-scale probabilistic generative model as follows: Remark 1. In the vanilla VAE, the objective function considers the reconstruction error and the KL-divergence to be of equivalent importance. The hyperparameter β is introduced to break this balance. More specifically, high values of β put more emphasis on the latent space approximation than on the reconstruction, expediting the learning of notable feature variations but bringing blurred minutiae. For example, in the channelized permeability experiment, we noted that high β is conducive to capturing the continuous channel structures while losing much fidelity in local variations.
To optimize the parameters θ and φ in neural networks, one could employ stochastic gradient descent (SGD) or other gradient-based optimization algorithms related to back propagation. Note that the second term in Eq. (15) is an expectation approximation using the Monte Carlo method that samples z (i,j) from q φ (z (i,j) |x (i) ). But the expectation computation involving sampling with a high variance will reflect on the gradient estimation, which leads to an unfavorable influence on optimization. To make it trainable and back propagate the gradient correctly, the reparameterization trick is introduced. The latent variables z are expressed by a differentiable transformation g φ ( , x) with respect to an auxiliary independent random variable . In the Gaussian distribution case, we let ∼ N (0, I), sampling z via such a g φ ( , x): where refers to an element-wise product. The illustration of forward and back propagation computation with the reparameterization trick is shown in Fig. 4. when n is large enough.
As discussed in Definition 3.1, we constructed the probabilistic models p θ (x|z) and q φ (z|x) to sample x given its corresponding latent variable z and to map the parameters

Multiscale deep generative model (MDGM)
In Section 3.3, we presented the single-scale parameter generative model in a probabilistic perspective and discussed how to use deep neural networks for its implementation. In this section, we enhance this model to a multiscale framework. Based on the parameter data generation procedure in Section 3.1, we start from the training data in the finest-scale and coarse grain them recursively to represent the training data in the coarsest-scale. We demonstrate the MDGM by explaining how to train it in the l-th scale as an example based on the assumption that we obtained the pre-trained encoder and decoder networks in the To construct the connection among various scales, we design a special latent space for the finer-scale. In particular, its latent variable z l inherits the latent variable z l−1 of the previous scale and augments it with an additional latent variable z l . These two variables are independent. This is in principle similar to the Bayesian approach followed in [13] where a hierarchical sparse grid approximation was used to represent an unknown parameter field at different scales.
As we discussed in Definition 3.2, we also need to train the recognition model q φ (z l |x l ) and generative model p θ (x l |z l ) in scale l with the same objective function as in Eq. (6). The difference with standard generative models is that the recognition model is divided into two parts in order to encode the parameter x with the independent latent variables z l−1 and z l . As shown in Fig. 5, the encoder network includes a pre-trained encoder network and an augmented encoder network. q φ (z l−1 |x l ) is computed by adopting the upscaling operator over x l firstly and then using the previous scale encoder network. It can be modeled by Since the upscaling operator U is deterministic, we can write the following: As in the single-scale probabilistic generative model in the last section, the variational lower bound L(θ, φ; x) consists of two parts as described in Eq.
and D KL q φ (z (i) |x (i) )||p θ (z (i) ) . The first part in the l-th scale is formulated as the reconstruction error introduced in Eq. (13), which involves the decoder model p θ l (x l |z l ) and its model parameter θ l . The second part is much distinct with p(z l ) being the target distribution in the KL-divergence. It is an isotropic multivariate Gaussian distribution as well, but it can be viewed as the joint distribution of (z l−1 , z l ), where z l−1 , z l are independent.
Thus p(z l−1 ) and p(z l ) are also isotropic multivariate Gaussian distributions. Combining with Eq. (18), we rewrite the KL-divergence as follows: where φ in the l-th scale is φ l containing φ l−1 and φ l , where φ l−1 denotes the parameters of the encoder network q(z l−1 |x l−1 ) in the (l − 1)-th scale and φ l denotes the parameters of the augmented encoder network q(z l |x l ) in the l-th scale. Recall that we need to employ the pre-trained encoder network q φ l−1 (z l−1 |x l−1 ) directly in the l-th scale training because it should ensure that the finer-scale generative model shares common latent variables z l−1 with the coarser-scale, so We also use gradient-based optimization algorithm to train the MDGM where the back propagation is only applied in the decoder network and augmented encoder network as blue arrows illustrated in Fig. 5. We can formulate the loss function for each training iteration like Eq. (15) as follows: whereμ lk (·) andσ lk (·) denote the k-th element in the mean and standard deviation, respectively, which are outputs of the l-th scale augmented encoder network, and µ l (·) is the output of the l-th scale decoder network. Here, m is still set to be 1 when the batch size n is large enough. d l is the dimension of the latent variable z l . The hyperparameterβ is determined by balancing the disentanglement of the parameter features and reconstruction.
In summary, assume that we obtained the probabilistic encoder and decoder models in the previous (l − 1) level. We can then train the augmented encoder network and decoder networks using the objective function in Eq. (21) so that the probabilistic encoder q φ l (z l |x l ) and the probabilistic decoder p θ l (x l |z l ) in the l-th scale are acquired. For details see Algorithm 2. The architectures of the augmented encoder and decoder networks are the same as those in the encoder and decoder networks for the single-scale model (refer to Appendix A).
We have introduced so far the probabilistic generative model in a single-scale in the last section and the concept of generative modeling across two scales in this section. Since these unsupervised learning models only need the unlabeled data, we only provide the finest-  encoder network, where z l−1 is sampled from q φ l−1 (z l−1 |x l−1 ) using the reparameterization trick. A z l−1 with high-variance will impact the training stability and convergence inflicting big noise in the latent variable z l . We address this problem in the training procedure by replacing z l−1 with its meanμ l−1 (x l−1 ).

Bayesian inversion with the DGM
In Bayesian inversion, the Metropolis-Hastings (MH) algorithm [56,57] is a popular MCMC method to approximate the posterior distribution [58]. As discussed in Definition 3.3, Algorithm 2 Training of the multiscale probabilistic generative model in the l-th scale ← Sample n noise variables from the Gaussian distribution p( ) 6: z n l ← Compute z l by the augmented encoder network in Eq. (16), and concatenate with z l−1 for each data 7: (21) with respect to θ l and φ l 8: θ l = θ −η∇ θ lL ← using gradient-based optimization algorithm (e.g. SGD or Adam) 9: φ l = φ−η∇ φ lL ← using gradient-based optimization algorithm (e.g. SGD or Adam) 10: end while 11: q φ l (z l |x l ) ← combine the pre-trained encoder q φ l−1 (z l−1 |x l−1 ) with q φ l (z l |x l ) Output: probabilistic encoder q φ l (z l |x l ), probabilistic decoder p θ l (x l |z l ).
we design a multiscale scheme aimed to realize the high-dimensional parameter estimation in inverse problems where the posterior distribution is approximated recursively. We utilize the previous coarser-scale estimation information as we proceed with inference in the next immediate finer-scale. Note that the coarsest-scale inference only involves single-scale information in which the parameterization employs the DGM given in Algorithm 1.
For the coarsest-scale posterior distribution, we are interested in combining the standard MH algorithm with the single-scale generative model to estimate the parameters x 1 ∈ R M 1 using the low-dimensional latent variable z 1 ∈ R d 1 , where M 1 is the dimension of the coarsestscale parameter x 1 that also determines the spatial discretization of the coarsest-scale forward model F 1 , and d 1 is the dimension of the coarsest-scale latent space. MCMC based on the pre-trained model p θ 1 (x 1 |z 1 ) allows us to explore the posterior approximation of z 1 .
Recall that the model p θ 1 (x 1 |z 1 ) is a Gaussian distribution with mean µ θ 1 (z 1 ) and standard derivation σ, where µ θ 1 (z 1 ) is the output of the decoder network in Algorithm 1 and σ is a hyperparameter that determines the noise level. In the test procedure, we can thus define x 1 = µ θ 1 (z 1 ) as a mapping from the latent space to the original parameter space. Then one can sample in the continuous latent space to produce various parameters x 1 using this mapping. In this way, the problem of inferring x 1 becomes equivalent to the one of inferring Given the observed data D obs , the forward function F 1 , the decoder model x 1 = µ θ 1 (z 1 ), and the prior distribution p(z 1 ), the unnormalized posterior distribution in the coarsest-scale is as follows: where the prior distribution p(z 1 ) is N (0, I).
For the above target distribution in the coarsest-scale, we only need to create a Markov chain with length N ite for the latent variable z 1 . The proposal distribution π q (z 1 |z (j) 1 ) can be the simple random walk [59] or a preconditioned Crank-Nicholson (pCN) algorithm [38,60,61]. To initialize the first state z from the prior distribution p(z 1 ). For each j ≥ 1, we first sample a candidate z 1 from the proposal distribution π q (·|z (j) 1 ) and then reconstruct the spatially-varying parameter x 1 by the decoder model x 1 = µ φ 1 (z 1 ). One can compute the likelihood function and the Metropolis acceptance ratio α using the forward model F 1 (x 1 ). The acceptance ratio α is defined as in MH algorithm: With probability α, the candidate sample z 1 is accepted, otherwise rejected. Specifically, in each iteration, the state includes not only the latent variable z 1 but also its corresponding spatially-varying parameter x 1 . The Markov chain typically takes some iterations to reach its stationary distribution, which depends on the target distribution. The first n b states of the Markov chain (burn-in stage) are thrown-away to ensure a valid approximation of the target distribution. Once the latent variable z 1 is estimated, its corresponding parameter samples {x (j) 1 } N ite j=n b can be viewed as samples from the posterior distribution π(x 1 |D obs ). A summary of the MH algorithm for Bayesian inverse modeling based on the deep generative model is given in Algorithm 3. This algorithm can be applied to solve any high-dimensional parameter estimation in the single-scale.
In the multiscale context, the coarsest posterior was approximated by Algorithm 3. But there are still two main problems to be resolved by invoking a finer-scale model. One is that the high-resolution parameter estimation that can present richer details is desired; the other is the inaccurate forward model introduces an epistemic uncertainty or model error, which will reflect on the parameter estimation. The model error δ l in the l-th scale is defined by This model error δ l will be zero if and only if we adopt the "true model" F L . The posterior distribution π(x L |D obs ) is the ultimate goal of our inverse problem. Draw z ∼ π q (·|z (j) ) 4: Compute the spatially-varying parameter x = µ(z ) by the decoder model 5: Compute the likelihood function by solving the forward model F(x )

6:
Compute the acceptance ratio if ρ < α then l . The acceptance ratio of the reject/accept scheme for the (j + 1)-th state in the MH-MDGM is still the same as that in the single-scale MH algorithm, i.e.
where the proposal π q (z l |z (j) l ) is a particular distribution since z l = (z l−1 , z l ) contains two independent random variables, which have different connotation in the MH-MDGM. The proposal is taken to have a factorized form as follows [38,62]: In the MDGM, z l is sampled from an isotropic Gaussian distribution, which is independent of z l−1 . The proposal distribution for z (j) l can be simplified as This can be a simple random walk or the pCN algorithm. Exploring z l provides additional information to that obtained at the (l−1)-scale allowing us to enrich the obtained local details and correct global features. One can set a big step size for z l in the proposal distribution to explore and detect various local patterns. In the MH-MDGM, the initial state z l can be sampled from the prior distribution p(z l ).
Based on the independence assumption, one can define a proposal distribution in the l-th scale for the latent variable z l−1 that connects adjacent scales as follows: Once z l−1 has been estimated in the previous-scale, some multiscale estimation methods like [13] will not correct and estimate it again in the finer-scale. However, note that the model error that resulted from the coarse-scale solver can impact the posterior estimation.
We  Message passing across scales depends on the encoder network π θ l−1 (z l−1 |x l−1 ). In order to generate x l for the fine-scale inference, we need to use {x which satisfies the detailed balance condition.
Proof. The transition kernel K(z l , z (j) l ) with the proposal distribution π q (z l |z (j) l ) in the MH algorithm is: where δ(z l , z l , it is obvious that the below detailed balance condition is satisfied: When l ). Based on Eqs. (27) and (28), the proposal distribution can be written as π q (z l |z Then we can write the following: and the detailed balance condition is satisfied.

Numerical Examples
The code and data for reproducing all examples reported in this paper can be found at https://github.com/zabaras/MH-MDGM. In this section, we discuss and compare the performance of the proposed multiscale method with that of the single-scale inference method.
Our focus is on the estimation of the permeability field in a single phase, steady-state Darcy flow. For any permeability field K on a 2D unit square domain S = [0, 1] 2 , the pressure field p and velocity field v are governed by the following equations: with boundary conditions v(s) ·n = 0, s ∈ Γ N ,  Draw z l−1 ∼ π q (·|z (j) l−1 ).

7:
Compute the spatially-varying parameter x l = µ θ l (z l ) by the decoder model 8: Compute the likelihood function by solving the forward model F l (x l )

9:
Compute the acceptance ratio α = min 1, π(z l |D obs )π q (z (j)   [63] with third-order Raviart-Thomas elements for the velocity v, and forth-order discontinuous elements for the pressure p. The computational cost for solving this forward model with different discretizations is reported in Table 1. To enforce the non-negative permeability constraint, we consider the log-permeability as the main parameter of interest in the inverse problem with x = log(K). The inverse problem for the above model is to infer the unknown log-permeability field given noisy pressure measurements at some sensor locations. In this paper, the exact log-permeability field x exact is defined in a 64 × 64 uniform grid. The generated data by the generative model are usually smooth and blurry compared to the original data due to information compression.
We consider 64 pressure observations that are uniformly located in [0.0625 + 0.125i, 0.0625 + 0.125i], i = 0, 1, 2, . . . , 7. A 5% independent additive Gaussian random noise is considered on these 64 pressure observations which are obtained by the above forward model given the exact log-permeability field.

Test problem 1: Gaussian Random Field (GRF)
For GRF based log-permeability data, one can adopt the KLE method for a reducedorder representation. Bayesian inference often works well on such low-dimensional inversion tasks. However, the KLE expansion requires a-priori knowledge of the length scale in the GRF covariance function. This information is of course not available in most practical applications. Data-driven methods like DGM do not require such restrictive assumptions.
The only prior knowledge available is based on the given training dataset.
In this section, we apply the proposed method on a two-scales scenario and provide a comparison with the reference case of the single-scale method. In the two-scales test case,

Multiscale dataset
As discussed in Section 3.4, we adopt an unsupervised learning method for parameter representation before addressing the solution of the inverse problem. In particular, one can train a generative model for the sampling of log-permeability fields. For such an unsupervised learning problem, we only need to obtain the dataset {x (i) } N i=1 based on historical sample information or the underlying prior distribution π(x). In this example, we assume that the log-permeability field is a Gaussian random field, i.e. x(s) ∼ GP (m(s), k (s 1 , s 2 )), where m(s) and k(s 1 , s 2 ) are the mean and covariance functions, respectively. s 1 = (x 1 , y 1 ) and s 2 = (x 2 , y 2 ) denote two arbitrary spatial locations. The covariance function k(s 1 , s 2 ) in this paper is taken as: where σ 2 log(K) is the variance, and l 1 and l 2 are the length scales along the x and y axes, respectively. In this example, we set m(s) = 1 and σ 2 log(K) = 0.5. As discussed earlier, it is hard to tackle the multiple length scales setup using the KLE method, while the deep generative model has a distinct superiority on the prior assumption since the dataset embodies these assumptions or information in a natural way. We assume that the length scales are not fixed in the prior distribution. Let the length scales be l 1 = l 2 = 0.2 + 0.01i, i = 0, 1, 2, . . . , 9.
The finest-scale parameter x 2 is uniformly discretized into an H × W = 64 × 64 grid. For each length scale, we generate 2500 samples for the training dataset {x where N is 25000. The exact log-permeability field x exact for the Bayesian inversion has never been seen in the training procedure. We set the 16 × 16 grid as the coarsest-scale in the parameter estimation. We adopt the upscaling operator in Eq. (9) with n e = 16 for each realization in The loss function in Eq. (21) is used for training the models q φ 2 (z 2 |x 2 ) and p θ 2 (x 2 |z 1 , z 2 ), where z 2 ∈ R 256 . The training time for the models q φ 2 (z 2 |x 2 ) and p θ 2 (x 2 |z 1 , z 2 ) is about 78 minutes.
Once the generative models are obtained in the different scales, the parameters can be generated in each scale by first sampling the corresponding latent variable, and then by decoding it using the corresponding generative model. Here, we illustrate the MDGM results in Fig. 6. For the pre-trained generative model p θ 1 (x 1 |z 1 ), we can sample the latent variables z 1 from N (0, I), and use these latent variables as input to the model p θ 1 (x 1 |z 1 ).
The first row of Fig. 6 shows 4 realizations with 16 × 16 resolution that are output of this decoder network. It is apparent that all of these images have distinct features in the spatial distribution of high-and low-value regions. However, they are very smooth so that one cannot capture any local information. Fortunately, this is a good choice to highlight the obvious global features, which are of great importance in Bayesian inversion. It is a fundamental and necessary requirement for spatially-varying parameter estimates to be consistent with the exact parameter in capturing global features.
The performance of the 64 × 64 resolution generative model p θ 2 (x 2 |z 1 , z 2 ) in this twoscales MDGM is shown in the second row of Fig. 6. The latent variable z 1 is the mean of the distribution q φ 1 (z 1 |x 1 ), i.e. z 1 = arg max q φ 1 (z 1 |x 1 ), where x 1 ∈ R 16×16 is the first image in this row. The latent variable z 2 is sampled from the Gaussian distribution N (0, I).
We observe that that the three samples generated by the model p θ 2 (x 2 |z 1 , z 2 ) have some particularities. Unlike samples in the first row, they keep similar spatial distribution of highand low-value regions that inherit from the first 16 × 16 resolution image x 1 while they are refined locally in a diverse manner. The model q φ 1 (z 1 |x 1 ) acts as a messenger in the MDGM. The information of low-resolution x 1 that captures global features is encoded by z 1 , which together with the random variables z 2 that is used to supplement local details are decoded by the generative model p θ 2 (x 2 |z 1 , z 2 ). This demonstrates that these two latent variables have different missions in the generative model, in particular z 1 plays an important role in the global feature generation, while z 2 contributes to local details. This disentangled representation [65,53,66,67] for local and global features assists in an accurate and efficient Bayesian multiscale estimation.

The inversion results and discussion
The reference experiment considered is the single-scale method. We use Algorithm 3 with the pre-trained model x = µ θ (z) to estimate directly π(z|D obs ), where x ∈ R 64×64 and z ∈ R 256 . All the implementations using MCMC treat the latent variables as random variables, and use the pre-trained model to generate x for evaluating the likelihood function.
The proposal distribution for the random variables applied in this paper is preconditioned Crank-Nicolson (pCN) that is defined below: where z and z are current and proposed next states, respectively, and ζ ∼ N (0, I). The step size of the random movement from the current state to a new position is controlled by the free parameter γ. We set γ be 0.08 for the first 50% and 0.04 for the last 50% states in the Markov chain for all implementations using Algorithm 3. We run 10000 iterations in MCMC to ensure its convergence. For all implementations of the MCMC algorithm, we collected the last 2000 states as the posterior samples. Fig. 7 shows the estimation results using the reference method. It can be seen that the integration of the deep generative model with MCMC leads to a reasonable estimation of the spatially varying parameter.

Exact Mean Standard deviation
Posterior sample 1 Posterior sample 2 Posterior sample 3 The multiscale parameter estimation is implemented next from coarse-to fine-scale. We use Algorithm 3 with pre-trained model x 1 = µ θ 1 (z 1 ) to estimate π(z 1 |D obs ), where x 1 ∈ R 16×16 and z 1 ∈ R 16 . A Markov chain with length 7000 is constructed for z 1 to estimate the coarse-scale Gaussian log-permeability field. This result explores various global patterns of the log-permeability. Fig. 8 shows that the estimation result is as expected. It is clear To correct and refine the parameter estimation, we use Algorithm 4 with the pre-trained model x 2 = µ θ 2 (z 2 ) to estimate π(z 2 |D obs ), where x 2 ∈ R 64×64 , z 2 = (z 1 , z 2 ), and z 1 ∈ R 16 , and z 2 ∈ R 256 . Since we have obtained the posterior samples in the coarse-scale, the estimation in the fine-scale takes advantage of the coarse estimation. We need to assign two proposal distributions for z 1 and z 2 . The two proposal distributions we used in this paper for Algorithm 4 are the pCN in Eq. (38) with different step sizes. The fixed step size γ for the low-dimensional latent variable z l−1 is 0.01, while the adaptive step size γ for the high-dimensional latent variable z l is 0.08 for the first 50% and 0.04 for the last 50% of the states in the Markov chain. The results in Fig. 9 indicate that the estimated parameter with 64 × 64 grid using the proposed method has even better performance than the reference single-scale results in Fig. 7. The details of the low-value region are closer to the exact log-permeability field.
To compare the estimation result with the exact field, we provide an illustration in the shortcoming of VAE [18] is that it generates blurry samples since the bottleneck layer captures a compressed latent encoding. However, we note that the result obtained from the proposed method is better than the reference method as its mean is much closer to the exact solution.
The main cost of the Bayesian inference using MCMC comes from the forward model evaluation. The forward model's computational cost in Table 2 suggests that the single-scale method takes about 1.9 times the computational cost of the proposed two-scale method. The acceptance rate of MCMC for all implementations is shown in Table 3. We can note that the proposed method has a much higher acceptance rate than the single-scale method in the estimated for the local refinement is assigned with a big step size. The reference method often rejects the proposed samples and gets trapped in local modes. However, in the finescale estimation in the proposed method, the main changes in the local features will lead to small changes in the likelihood function so that most of the proposed samples are accepted. The convergence of the Markov chain used to estimate parameters with a 64 × 64 grid is of concern since the forward evaluation on such a scale is very expensive. To assess the convergence in the desired scale with a 64 × 64 grid, we employ three metrics using variables available in the iterative process. Since the observation data (noisy pressure measurements) is the only basis for parameter estimation, we compute the misfit between the observation data and its corresponding prediction using the sum of squared residuals of observations (SSR obs ) as the iteration proceeds: where N obs = 64 in this paper, and F(µ(z)) (i) and D (i) obs are the i-th predicted pressure value and its corresponding true observation, respectively. Ideally, we expect that the value of SSR obs is close to 0, which suggests there is no discrepancy between the predictions and observations. But this cannot be realized even using the exact input for F due to measurements noise. To remove the impact of noise, we use an enhanced metric i.e. the normalized sum of squared weighted residual (N SSW R) [16,68]: where x exact is the exact log-permeability field, and σ (i) n is the standard deviation of Gaussian random noise imposed in the i-th observation. The SSW R metric is normalized by the SSW R ref . Thus, the inversion process has converged when the N SSW R value is close to 1. Since the target is to estimate the log-permeability field based on the given observations, we also consider the evaluation of the mismatch between the predicted and the exact log-permeability fields using the sum of squared residuals of parameter (SSR para ) as shown below: where M = 4096, since the desired scale is discretized to 64 × 64, and µ(z) (i) and x (i) exact are the i-th value of the predicted and the exact log-permeability fields, respectively. Fig. 11 shows the convergence results in the GRF case. The comparison of the proposed method with the single-scale method is also given in Fig. 11. It can be seen that better performance in all evaluation metrics is produced by the proposed method. A big distinction takes place in the decrease of SSR para illustrated in Fig. 11(c). We show the sampled logpermeability field with 64 × 64 grid at different iterations in the Markov chain in Fig. 12.
For the single-scale method, each dimension of the estimated latent variable z ∈ R 256 has equivalent importance for the generation of the log-permeability. The random walk in such a high-dimensional space has difficulty in efficiently exploring the space and identifying a good estimate and its uncertainty. Together with a random initial state sampled from the prior distribution, the exploration takes a long time to reach the stationary distribution. The second row in Fig. 12 presents the state evolution for the two-scale method. The first-state that inherited the coarse-scale estimation captures well most of the non-local features of the exact log-permeability field. Thus with the iterations shown, we only need to correct the local features to explore the posterior distribution.

Test problem 2: Non-Gaussian Random Field
In the second test problem, we consider a channelized log-permeability field as the exact parameter in the inversion experiment. For such a non-Gaussian permeability field is often difficult to obtain a good parameterization using conventional methods such as sparsegrid approximations [13], wavelets [14], or principal component analysis (PCA) [69,70,71].  Furthermore, the Bayesian inference using random walk MCMC based on recent reported learning-based parameterization methods [16,15,46,22] for high-dimensional non-Gaussian parameters is not a good choice. We show the benefits of the multiscale method in the parameterization of non-Gaussian random fields and Bayesian inference. In this test example, we demonstrate the proposed method with two-and three-scales scenarios and provide comparisons with the reference case of a single-scale method.

Multiscale dataset
Suppose that the prior information for the channel location before any measurement is from the image [22] with size of 2500×2500 shown in Fig. 13 Fig. 13(b). The binary image depicts channels with white regions. We assume that the channelized regions have different log-permeabilities resulting in high-and low-permeability values in white and black regions, respectively. We set the log-permeability values for each region by independently sampling from two Gaussian Random Fields (GRFs). The means of the GRFs for the high-permeability channelized regions and low-permeability regions are 4 and 0, respectively. The covariance function in Eq. (37) with length scales l 1 and l 2 equal to 0.3 is applied. The variance is 0.5 for both GRFs. By imposing two GRFs for the binary image samples, we can generate the training data as shown in Fig. 13(c). We assume the generated samples are i.i.d. sampled from the underlying prior distribution π(x). We take 40000 samples from generated 46208 realizations as training data. The exact log-permeability field x exact for Bayesian inversion is sampled from the remaining samples, which has never been seen in the training procedure. The x exact we used in this test problem is the image shown in Fig. 13(c).  with n e = 4 over the original 64 × 64 grid data, and apply it again with n e = 4 in the generated 32 × 32 grid data to obtain the 16 × 16 grid data. An example with different discretizations has been illustrated in Fig. 2.

Training and results of the MDGM
Using the generated training datasets, we applied the same neural network used in the GRF case to train the MDGM. We will discuss the training procedure and performance below for the two-and three-scales cases. The reference single-scale generative model p θ (x|z) requires 3.15 hours to train 50 epochs using the loss function in Eq. (15) and Algorithm 1, Two scales (16−64). In this example, the x 1 , x 2 denote the parameters with 16×16 grid and 64 × 64 grid, respectively. For the x 1 , the generative model is trained using Algorithm 1 with β = 1 in Eq. (15). It takes about 37 minutes for 30 training epochs. We obtained p θ 1 (x 1 |z 1 ) and q φ 1 (z 1 |x 1 ), where x 1 ∈ R 16×16 and z 1 ∈ R 16 . The model q φ 1 (z 1 |x 1 ) is the input to Algorithm 2 used for training the finer-scale generative model. We use the pre-trained model p θ 1 (x 1 |z 1 ) to reconstruct the parameters to compute the likelihood function when we use Algorithm 3 to estimate π(z 1 |D obs ), where z 1 ∈ R 16 . The generated samples using the model p θ 1 (x 1 |z 1 ) are shown in the first row of Fig. 14 Three scales (16 − 32 − 64). In the second example with three scales, the parameters with 16×16, 32×32 and 64×64 grids are denoted by x 1 , x 2 , x 3 , respectively. The model p θ 1 (x 1 |z 1 ) is the same as the above two-scales example as shown in the first row of Fig. 14(a). The first row of Fig. 14(b) gives the results of the model is encoded from the first image in this row using the model q φ 1 (z 1 |x 1 ), and z 2 ∈ R 64 is sampled from N (0, I). It takes about 1.3 hours to train 50 epochs withβ = 0.7 in Eq. (21).
Compared with the two-scales case, one can notice that the generated samples have highly consistent channels that inherit from the first image but maintain local diversity. In order to sample the desired parameters with 64 × 64 grid, we also adopt Algorithm 2 withβ = 0.7 in Eq. (21) to train the model p θ 3 (x 3 |z 2 , z 3 ), where z 2 ∈ R 80 , z 3 ∈ R 256 , and the training time

The inversion results and discussion
Once all generative models are trained, we can use them in the Bayesian inversion process.
To assess the efficiency and accuracy of the proposed multiscale method, we first evaluate the posterior distribution using the reference single-scale method. Algorithm 3 using the pre-trained model x = µ θ (z) as the input can perform the random walk in the latent space, resolution random samples using the model p θ2 (x 2 |z 1 , z 2 ). We let z 1 = arg max q φ1 (z 1 |x 1 ), and x 1 is the first image in this row, and z 2 ∈ R 64 is randomly sampled from N (0, I). Second row: 64 × 64 resolution random samples using the model p θ3 (x 3 |z 2 , z 3 ). We let z 2 = arg max q φ2 (z 2 |x 2 ), and x 2 is the first image in this row, and z 3 ∈ R 256 is randomly sampled from N (0, I).
where x ∈ R 64×64 and z ∈ R 256 . Although the dimension of the latent variable z is still high, the convergence can be realized by constructing a Markov chain with the length of 30000.
The inferred results are depicted in Fig. 15. The estimated mean and samples mostly match the exact channels and some important local details. Three scales (16 − 32 − 64). To discuss the impact of the number of scales, we also studied the three-scale estimation with pre-trained three-scale generative models. The coarsest scale result was shown in Fig. 16. To refine the estimated results in the scale with 32 × 32 grid, we run 7000 iterations in Algorithm 4 with the pre-trained model x 2 = µ θ 2 (z 2 , z 1 ), where

Exact
x 2 ∈ R 32×32 , z 2 ∈ R 64 , and z 1 ∈ R 16 . Unlike the estimation in the coarsest-scale that can only identify the location of the channels, the refined results with the 32 × 32 grid can also provide a good estimation regarding the spatial distribution of high-and low-values as shown in Fig. 18. Since most salient features were captured in previous scales, the desired scale estimation becomes much easier. We only need to perform 5000 iterations in Algorithm 4 with the pre-trained model x 3 = µ θ 3 (z 3 , z 2 ) to guarantee its convergence, where x 3 ∈ R 64×64 , z 3 ∈ R 256 , and z 2 ∈ R 80 . Fig. 19 summarizes the final three-scales inference results where we can notice that the estimation is more accurate than the previous two experiments in The computational cost of the forward model in each experiment is given in Table 4. among them is the computational cost in the desired scale. We designed the multiscale scheme to reduce the computational cost in the fine-scale with some additional coarse-scale forward evaluations. As a consequence, this leads to significant computational savings during Bayesian inference. One can apply the proposed method to a computationally more intensive model. Correspondingly, more reduced computational time can be expected. The acceptance rate of each MCMC implementation is shown in Table 5. For the reference single-scale experiment, each dimension of the latent variable has equivalent importance in the generation of the log-permeability. At the same time, channels are the most salient features in this Bayesian inversion example. To infer a high-dimensional latent variable one needs to explore a high-dimensional state space, which will lead to a very low acceptance rate. For the multiscale method, the acceptance rate in the coarsest-scale is still very low. Based on the coarse-scale estimation, applying Algorithm 4 to refine and correct the estimation will become much easier since two proposal distributions with different step sizes provide an informed and efficient exploration. Using a small step size for the low-dimensional latent variable to correct the global features and using a big step size for the high-dimensional latent  Figure 19: The three-scale estimation result with the desired 64×64 grid. The results involve inference across scales. Inference was performed with respect to z 2 ∈ R 80 and z 3 ∈ R 256 is performed using Algorithm 4.
variable to explore local features is goal-oriented, which not only improves the acceptance rate but also leads to a better estimation. As with the Gaussian case, we use three evaluation metrics to assess the convergence in the desired scale with the 64 × 64 grid. Fig. 21 indicates a better performance in the channelized log-permeability estimation. It is obvious that the reference single-scale method takes a long exploration to converge to the stationary distribution. As discussed earlier,  the inferred latent variable of the single-scale inference is still high-dimensional, and each dimension keeps equivalent importance in the log-permeability generation. In such scenarios, it is easy to get trapped in local modes. Fig. 22 provides the state evolution of the Markov chain to infer the log-permeability on the 64 × 64 grid. The single-scale method captures well the true solution at about the 15000-th iteration, while the multiscale method only needs to refine the coarse-scale estimation for fine-scale inference. By greatly reducing the fine-scale forward model evaluations, the computational burden in Bayesian inverse problems is reduced. For a parameter like the channelized log-permeability with obvious multiscale characteristics, it can readily be seen that the three-scales experiment performs better than the two-scales experiment with respect to stability, efficiency and accuracy.

Conclusions
In this work, we introduced a novel multiscale parameter estimation framework for   To make the other (l − 1) elements (i.e. z 1 , z 2 , . . . , z l−1 ) in z l to be of consistent size with z l , we use the Upsample operator 3 in the Pytorch library [77] over these elements and then concatenate 4 all of them as input z l ∈ R C×H×W for the l-th scale decoder networks.
The scale factor in the Upsample operator depends on the output and input sizes. Their sizes satisfy the following relationship: