title: Deep generative models in inversion: a review and development of a new approach based on a variational autoencoder
authors: Jorge Lopez-Alvis Eric Laloy Frédéric Nguyen Thomas Hermans
text: 
  
  
  
  [t1]Jorge Lopez-Alvis conceived the main idea of the paper, developed and tested the proposed method (including its software implementation) and wrote a first draft of the manuscript. Eric Laloy conceived the main idea of the paper, performed some software tests, suggested improvements to the proposed method and edited the first draft. Frédéric Nguyen and Thomas Hermans are supervisors of the first author, conceived the main idea of the paper, suggested improvements to the proposed method, obtained funding for the research and edited the first draft.
  
  
  address1,address2]Jorge Lopez-Alvismycorrespondingauthor
  [mycorrespondingauthor]Corresponding author
  jlopez@uliege.be
  
  address3]Eric Laloy
  
  address1]Frédéric Nguyen
  
  address2]Thomas Hermans
  
  [address1]Urban and Environmental Engineering, Applied Geophysics, University of Liege, Belgium
  [address2]Department of Geology, Ghent University, Belgium
  [address3]Engineered and Geosystems Analysis, Institute for Environment, Health and Safety, Belgian Nuclear Research Center, Belgium
  
  
  When solving inverse problems in geophysical imaging, deep generative models (DGMs) may be used to enforce the solution to display highly structured spatial patterns which are supported by independent information (e.g. the geological setting) of the subsurface. In such case, inversion may be formulated in a latent space where a low-dimensional parameterization of the patterns is defined and where Markov chain Monte Carlo or gradient-based methods may be applied. However, the generative mapping between the latent and the original (pixel) representations is usually highly nonlinear which may cause some difficulties for inversion, especially for gradient-based methods. In this contribution we review the conceptual framework of inversion with DGMs and study the principal causes of the nonlinearity of the generative mapping. As a result, we identify a conflict between two goals: the accuracy of the generated patterns and the feasibility of gradient-based inversion. In addition, we show how some of the training parameters of a variational autoencoder, which is a particular instance of a DGM, may be chosen so that a tradeoff between these two goals is achieved and acceptable inversion results are obtained with a stochastic gradient-descent scheme. A test case using truth models with channel patterns of different complexity and cross-borehole traveltime tomographic data involving both a linear and a nonlinear forward operator is used to assess the performance of the proposed approach.
  
  
  
  deep generative model geophysical inversion geological prior information variational autoencoder stochastic gradient descent deep learning
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  A common task in the geosciences is to solve an inverse problem in order to obtain a model (or image) from a set of measurements sensing a heterogeneous spatial domain. When measurements are sparse, the inverse problem is usually ill-posed and its solution non-unique. In such cases it is possible to constrain the solution to be found only among models with certain spatial patterns. In practice, such patterns are supported by independent (prior) information of the sensed domain (e.g. knowledge of the geological setting) and used with the aim of appropriately reconstructing heterogeneity. Classical regularization may be used to impose certain structures to the solution <cit.> but these are generally too simple to be realistic. Recently, the use of deep generative models (DGMs) to constrain the solution space of inverse problems has been proposed so that resulting models have specific spatial patterns <cit.>. DGMs can deal with realistic (natural) patterns which are not captured by classical regularization or random processes defined by second-order statistics <cit.>. In this way, inversion with DGMs provides an alternative to inversion with either multiple-point geostatistics (MPS) <cit.> or example-based texture synthesis (ETS) <cit.>.
  
  All the previously mentioned methods generally rely on gridded representations for the models (i.e. by dividing the spatial domain in cells or pixels). However, a key difference between them is that MPS and ETS directly extract spatial patterns from a training image (or exemplar) <cit.> whereas DGMs require an initial learning phase in which samples of the patterns are used as a training set—such samples may be obtained by different means, e.g. cropped from a training image or taken directly from a large set of images. While any of the methods (DGMs, MPS or ETS) may be used with inversion, DGMs rely on a so called latent space where a low-dimensional parametric representation (referred to as latent vector or code) is defined and where Markov Chain Monte Carlo (MCMC) or gradient-based methods may be applied. Note that when using gridded representations each model may be seen as a vector in "pixel" space (a space where each pixel is one dimension). Since only highly structured spatial patterns are allowed, vectors will only occupy a subset of the pixel space. This subset defines a manifold of lower dimensionality than the pixel space <cit.> and the latent space is simply a low-dimensional space where such manifold is mapped. Most inversion methods require a perturbation step to search for models that fit the data but such a step is not straightforward to compute for highly structured patterns <cit.>. The latent space of DGMs provides a useful frame to compute a perturbation step <cit.> or even a local gradient-descent direction <cit.> which generally results in better exploration of the posterior distribution and/or faster convergence compared to inversion with MPS or ETS.
  
  So far, inversion with DGMs has been done successfully with regular MCMC sampling methods <cit.>. However, when applicable, gradient-based methods may be preferred given their lower computational demand. Gradient-based deterministic inversion with DGMs has been pursued with encouraging results (Richardson, 2018, Laloy et al., 2019), however, convergence to the true model was shown to be dependent on the initial model. In the framework of probabilistic inversion, MCMC methods that use the gradient to guide the sampling in the latent space have shown to be less prone to get trapped in local minima than gradient-based deterministic methods while they are also expected to reach convergence faster than regular MCMC <cit.>. A different inversion strategy that has also been applied successfully with DGMs and has a relatively low computational cost is the Ensemble Smoother <cit.>.
  
  Recently, <cit.> studied the difficulties of performing gradient-based deterministic inversion with a specific DGM. They concluded that the non-linearity of their generative function or "generator" (i.e. the mapping from the latent space to the pixel space) was high enough to hinder gradient-based optimization, causing the latter to often fail in finding the global minimum even when the objective function was known to be convex (in pixel space). Such high non-linearity is expected since DGMs approximate highly complex patterns by using low-dimensional inputs usually generated by simple probability distributions (e.g. normal or uniform distributions). The complexity of realistic patterns means that the corresponding manifold will generally have both a curvature and a topology that is radically different from the region (or subset) defined indirectly in the latent space by the chosen probability distribution. Then, the generative function has to deform this region (when mapping it to the pixel space) in such a way as to approximate (or cover) the manifold as close as possible. The combined deformation needed to curve the region and to approximate its topology causes the generative function to be highly nonlinear. In order to approximate manifolds of realistic patterns, most common DGMs involve (artificial) neural networks with several layers and non-linear (activation) functions. Recently, it has been shown that deep neural networks with a ReLU activation function are able to change topology of the input <cit.> so, besides nonlinearity, the generator of DGMs may also induce changes in topology when mapping from the latent space. When the sole purpose of the DGM is for generating new samples, high nonlinearity and induced changes in topology are not important but they might be an issue when the DGM is used for additional tasks, such as inversion.
  
  For a specific subsurface pattern, the degree of non-linearity and the changes induced in topology by the generative function may be controlled mainly by its architecture and the way it is trained <cit.>. Regarding difference in training, two common types of DGMs can be distinguished: generative adversarial networks (GANs) <cit.> and variational autoencoders (VAEs) <cit.>—in both cases training generally takes the form of optimizing a loss function. GANs and VAEs require specification of a probability distribution in the latent space and an architecture for the discriminator or encoder (respectively) in addition to the one for their generators. They might also require other parameters to be specified such as the weights on the different terms of the loss function. Frequently, some of these choices use default values, but generally all of them may affect the degree of non-linearity of the generator <cit.>.
  
  Given all these possible controls for learning the generator it is interesting to investigate whether there exist some combinations of such controls that allow both for a good reproduction of the patterns and a good performance of less computationally demanding gradient-based inversion. In this work, we review some of the difficulties of performing inversion with DGMs and show how to obtain a well-balanced tradeoff between accuracy in patterns and applicability of gradient-based methods. In particular, we propose to use the training choice of a VAE as DGM and to select some of its parameters in order to achieve good results with gradient-based inversion. Then, we compare this to the training choice of a GAN that has been tested with gradient-based inversion in prior studies <cit.>. Furthermore, we show that since the resulting VAE inversion is only mildly nonlinear, modified stochastic gradient-descent (SGD) methods are generally sufficient to avoid getting trapped in local minima and provide a better alternative than regular gradient-based methods while also retaining a low computational cost.
  
  The remainder of this paper is structured as follows. Section <ref> explains DGMs and their conceptualization as approximating the real (pattern) manifold. In Section <ref> the use of DGMs to represent prior information in inversion and the difficulties of performing gradient-based inversion are reviewed. Sections <ref> and <ref> show how to use a VAE and SGD to cope with some of the mentioned difficulties. Then, Section <ref> shows some results of the proposed approach. Section <ref> discusses the obtained results and points to some remaining challenges. Finally, Section <ref> presents the conclusions of this work.
  
  
  
  § METHODS
  
  
  
  
   §.§ Deep generative models (DGM) to represent realistic patterns.
  
  
  
  The term "deep learning" generally refers to machine learning methods that involve several layers of multidimensional functions. This general "deep" setting has been shown to allow for complex mappings to be accurately approximated by building a succession of intermediate (simpler) representations or concepts <cit.>. Consider, for instance, deep neural networks (DNNs) which are mappings defined by a composition of a set of (multidimensional) functions ϕ_k as:
  
  
      𝐠(𝐱) = (ϕ_L ∘…∘ϕ_2 ∘ϕ_1)(𝐱)
   
  
  where 𝐱 is a multidimensional (vector) input, k = {1, …, L} denotes the function (layer) index and composition follows the order from right to left. Furthermore, each ϕ_k is defined as:
  
  
      ϕ_k(ξ)=ψ_k(𝐌_kξ+𝐛_k)
  
  
  in which ψ_k is a (nonlinear) activation function, 𝐌_k and 𝐛_k are vectors of weights and biases, respectively, and ξ denotes the output of the previous function (layer) ϕ_k-1 for k>1 or the initial input 𝐱 for k = 1. Then, training the DNN involves estimating the values for the all the parameters θ={𝐌_k,𝐛_k | 1 ≤ k ≤ L} where each 𝐌_k or 𝐛_k may be of different dimensionality depending on the layer. In practice, the number of parameters θ for such models may reach the order of 10^6, therefore training is achieved by relying on autodifferentiation <cit.> and fast optimization techniques based on stochastic gradient descent (SGD) <cit.>, both usually implemented for and run in highly parallel (GPU) computing architectures.
  
  A deep generative model (DGM) is a particular application of such deep methods <cit.>. In a DGM a set of training examples 𝐗={𝐱^(i)| 1 ≤ i ≤ N} and a simple low-dimensional probability distribution p(𝐳) are used to learn a model 𝐠(𝐳) that is capable of generating new samples of 𝐱 (which are consistent with the training set) by using as input samples from p(𝐳). This can be written as:
  
  
      𝐱=𝐠(𝐳), 𝐳∼ p(𝐳)
  
  
  where 𝐠(𝐳) is referred to as the "generator" and 𝐳 denotes a vector of latent variables or "code". While the training (and generated) samples 𝐱 are usually represented in a high-dimensional space ℝ^D, the probability distribution p(𝐳) is defined in a low-dimensional space ℝ^d. The space ℝ^D is often referred to as "ambient space" while the space ℝ^d is called the "latent space". Fig. <ref>a,c shows a schematic representation where the dimensionality of the ambient space is D=3  and the one of the latent space is d=2. A typical application of DGMs is the generation of images <cit.> for which the ambient space is just the pixel space. Gridded representations of subsurface models may be seen as two- or three-dimensional images of the subsurface.
  
  
  
  The underlying assumption in DGMs is that real-world data are generally structured in their high-dimensional ambient space ℝ^D and therefore have an intrinsic lower dimensionality—such assumption is known in machine learning literature as the manifold hypothesis <cit.> because it states that high-dimensional data usually lie on (or lie close to) a lower-dimensional manifold ℳ⊂ℝ^D. For instance, when studying a subsurface region it is usually assumed that geological processes gave it certain degree of structure then, to allow for a flexible base on which to represent the distribution of the different subsurface materials, the region is usually divided in pixels (or cells) within each of which the material is assumed to be homogeneous. Such gridded representation "lives" in the high-dimensional pixel space (the ambient space) but since it has some structure there should be a lower dimensional space (the latent space) where the same distribution of subsurface materials might be represented. Technically, while both the latent space ℝ^d and the manifold ℳ are usually low-dimensional, they may differ in dimensionality and/or the manifold may only occupy a certain portion of the latent space (e.g. the shaded region in Fig. <ref>c).
  
  Considering the manifold assumption described above, a DGM may be regarded as a model to implicitly approximate the "real" manifold ℳ by generating samples that closely follow such manifold, i.e. that lie on an approximate manifold ℳ' (Fig. <ref>b). Samples of this approximate manifold are generated by sampling first from a simple probability distribution p(𝐳) in latent space (e.g. a normal or uniform distribution) and then passing them through the generator 𝐠(𝐳). Since the probability distribution p(𝐳) defines indirectly a region (or subset) in latent space that generally has a different curvature and topology than the real manifold, the generator 𝐠(𝐳) must be able to approximate both curvature and topology when mapping the samples of p(𝐳) to ambient space. This generally requires the generator to be a highly nonlinear function. As an instance, consider the case of certain spatial patterns whose real manifold is a highly curved surface with "holes" in ambient space and the (input) region defined by a uniform p(𝐳) is a (flat) plane in a two-dimensional latent space. Regarding their topological properties, one technically says that this plane is simply connected while the real manifold is not <cit.>. Then, the generative function has to deform this plane in such a way as to approximate (or cover) the real manifold as close as possible. An important property of DGMs is that since a probability distribution in latent space is used, the sample "density" of such plane (and its mapping) also plays an important role. For instance, the generative function may approximate the "holes" of the real manifold by creating regions of very low density of samples when mapping to ambient space (to picture this one can imagine locally stretching a flexible material without changing its curvature). The combined deformation needed to curve the plane and to "make" the holes causes the generative function to be highly nonlinear. Note that when considering a DGM that uses a DNN with ReLU activation functions as generator 𝐠(𝐳), it is also possible for 𝐠(𝐳) to change topology of the input by "folding" transformations <cit.>.
  
  While one should always strive to accurately approximate the real manifold, since a finite set of training samples is used a tradeoff between accuracy and diversity in the generated samples may be a better objective. Indeed, the use of the prescribed probability distributions is done to continuously "fill" the space between the samples and therefore generate samples of a continuous manifold. Recent success—in terms of accuracy and diversity of generated samples—has been achieved with two DGMs that are based on deep neural networks (DNNs): generative adversarial networks (GANs) <cit.> and variational autoencoders (VAEs) <cit.>. The generator g(𝐳) on both strategies is a mapping from low-dimensional input 𝐳∈ℝ^d to high-dimensional output 𝐱∈ℝ^D. In contrast, the mappings corresponding to the discriminator and the encoder take high-dimensional inputs 𝐱 and return low-dimensional outputs.
  
  
  
   §.§ Gradient-based inversion with DGMs
  
  
  
  Consider a survey or experiment for which a vector of noisy measurements 𝐲 = (y_1,…,y_Q)^T ∈ℝ^Q of a physical process is available. A simplified description of the process may be expressed by a (mathematical) forward operator 𝐟:ℝ^D →ℝ^Q that takes as input a subsurface model vector 𝐱 = (x_1,…,x_D)^T ∈ℝ^D (obtained by discretizing the spatial distribution of physical properties) and outputs a simulated response 𝐟(𝐱). Commonly, this operator is in the form of a (numerical) discretization of a set of PDEs describing the process under study and is an approximation of the real process. As a result of such approximation and the use of noisy data, a noise term η is added to the simulation to represent total uncertainty. Then, the relation between the operator and the measurements may be written as <cit.>:
  
  
      𝐲=𝐟(𝐱)+η
  
  
  The corresponding inverse problem or inversion of Eq. (<ref>), aims to obtain an estimation of the vector 𝐱 from the (noisy) data 𝐲. Deterministic inversion does so by optimizing a misfit function γ(𝐱) that is usually given in the form of a distance function between simulated response 𝐟(𝐱) and data 𝐲, e.g. by the l_2 norm:
  
  
      γ(𝐱) = 𝐟(𝐱)-𝐲^2
  
  
  In this way, traditional gradient-based inversion requires the gradient ∇_𝐱γ(𝐱) whose elements are:
  
  
      [∇_𝐱γ(𝐱)]_i = ∂γ(𝐱)/∂ x_i
  
  
  and are computed by considering Eq. (<ref>) together with the chosen misfit.
  
  DGMs may be used with inversion of subsurface data 𝐲 to obtain geologically realistic spatial distributions of physical properties 𝐱 <cit.>. While this is also possible with traditional deterministic inversion where a regularization term is added directly in Eq. (<ref>) (i.e. in ambient space) to obtain models with the imposed structures that minimize the misfit <cit.>, DGMs are more flexible because they can enforce different structures provided they are appropriately trained. In the DGM setting, the low-dimensional samples 𝐳 that input to the generator 𝐠(𝐳) may be seen as defining a low-dimensional parameterization (or encoding) of realistic patterns 𝐱 and therefore exploration of the set of feasible models may be done in the latent space ℝ^d, as long as the search is done within the region where the approximated manifold ℳ' is defined (depicted by shading in Fig. <ref>c).
  
  Since the misfit γ(𝐱) is typically defined in ambient space ℝ^D (e.g.  in Fig. <ref>a), gradient-based inversion with DGMs may be seen as optimizing the intersection of γ(𝐱) with the approximate manifold ℳ' (Fig. <ref>e). Such intersected misfit is mapped into the latent space (Fig. <ref>f) and may be expressed as γ(𝐠(𝐳)). Also note that when probability distributions p(𝐳) with infinite support are used (e.g. a normal distribution), one can guide the search in the latent space by adding controlling (regularization) terms to the mapped misfit <cit.> and the resulting objective function may be written as:
  
  
      ζ(𝐳)    = γ(𝐠(𝐳)) + λ R(𝐳) 
         = 𝐟(𝐠(𝐳))-𝐲^2 + λ R(𝐳)
  
  
  where R(z) is a regularization term defined in the latent space and λ is the corresponding regularization factor. 
  
  In practice, no exhaustive mapping has to be done and the gradient ∇_𝐳ζ(𝐳) is only computed for the points in latent space where optimization lands in each iteration (in Fig. <ref>f the gradient is represented for one iteration). The gradient ∇_𝐳ζ(𝐳) is computed by adding a derivative layer corresponding to ∇_𝐱γ(𝐱) to the autodifferentation that was set up for 𝐠(𝐳) while training the DGM <cit.>. Such autodifferentiation setup may be seen as implicitly obtaining the jacobian 𝐉(𝐳) of size D × d whose elements are:
  
  
      [𝐉(𝐳)]_i,j = ∂ g_i(𝐳)/∂ z_j
  
  
  Then, the gradient ∇_𝐳ζ(𝐳) is obtained from Eq. (<ref>) by using the chain rule given by the product of Eqs. (<ref>) and (<ref>):
  
  
      ∇_𝐳ζ(𝐳)    = ∇_𝐳γ(𝐠(𝐳)) + λ∇_𝐳 R(𝐳) 
         = 𝐉(𝐳)^T ∇_𝐱γ(𝐱) + λ∇_𝐳 R(𝐳)
  
  
  The latter may also be done implicitly by incorporating directly in the autodifferentiation framework, e.g. putting it on top of the so called computational graph <cit.>.
  
  Assuming the considered misfit function γ(𝐱) is convex in ambient space ℝ^D (as depicted by concentric contours in Fig. <ref>a), difficulties to perform gradient-based deterministic inversion may arise due to the generator 𝐠(𝐳) <cit.>. We propose that such difficulties arise because the generator (1) is highly nonlinear and (2) changes the topology of the input region defined by p(𝐳). Both of these properties often cause distances (between samples) in latent space to be significantly different than distances in ambient space. Consider again the example of a real manifold that is a highly curved surface with "holes" in it and a uniform distribution p(𝐳) is used as input to the generator, then the latter might be able to approximate both the curvature and the holes at the cost of increasing nonlinearity and/or changing topology. When considering this backwards—e.g. when mapping the misfit function γ(𝐱) in the latent space—the approximation of both high curvature and differences in topology often translate in discontinuities or high nonlinearities because a continuous mapping onto the uniform distribution is enforced. This results in high curvature being effectively "flattened" and holes effectively "glued", both of which cause distances to be highly distorted. In this work, we will call a generator "well-behaved" when it is only mildly nonlinear and preserves topology.
  
  Both the generator's nonlinearity and its ability to change topology, may be controlled by two factors: (1) the generator architecture (type and size of each layer and total number of layers) and (2) the way it is trained (including training parameters). If the goal is to perform gradient-based inversion with DGMs, one should try to preserve convexity of γ(𝐱) as much as possible when mapping it to the latent space as γ(𝐠(𝐳)) while not degrading the generator's ability to reproduce the desired patterns. To aid in preserving such convexity, we propose to enforce the generator 𝐠(𝐳) to be well-behaved. This means that the generator will approximate the real manifold ℳ with a manifold ℳ' with a moderate curvature and whose topology is the same as the region defined in latent space by p(𝐳). By enforcing a moderate curvature manifold, local oscillations that may give rise to local minima (as those shown in Fig. <ref>d) but only have minimum impact in pattern accuracy are avoided in the approximate manifold ℳ' (the local minima are no longer present in Fig. <ref>e). In turn, when the generator is encouraged to preserve topology no more local minima should arise in ℝ^d than the ones resulting from intersecting γ(𝐱) with the approximate manifold ℳ' in ℝ^D (note e.g. there is one local minima in both Fig. <ref>e,f). The latter is in line with the proposal of <cit.>, where they argue that for the purpose of representation learning (which basically means learning encodings that are useful for other tasks than just generative modeling) the mapping should preserve topology.
  
  GANs often produce highly nonlinear generators that do not preserve topology, which may result in challenging inversion in the latent space. <cit.> provide an example of how architecture of a GAN is set to obtain a relatively well-behaved generator 𝐠(𝐳). They propose to use a model called spatial generative adversarial network (SGAN) <cit.> that enforces different latent variables to affect different local regions in the ambient space. Their architecture results in a high compression (lower dimensionality of the latent space) and controls nonlinearity which allowed them to successfully perform MCMC-based inversion in the latent space. However, gradient-based deterministic inversion performed with the same DGM was shown to be highly dependent on the initial model <cit.> pointing towards the existence of local minima. In this work we aim for robust gradient-based inversion in latent space by considering a VAE, the other predominant type of DGM, and its ability to produce a well-behaved generator.
  
  
  
   §.§ Variational autoencoder (VAE) as DGM for inversion
  
  
  
  A variational autoencoder (VAE) is the model resulting from using a reparameterized gradient estimator for the evidence lower bound while applying (amortized) variational inference to an autoencoder, i.e. an architecture involving an encoder and a decoder which are both (possibly deep) neural networks <cit.>. To train a VAE one uses a dataset 𝐗={𝐱^(i)| 1 ≤ i ≤ N} where each 𝐱^(i) is a sample (e.g. an image) with the desired patterns and then maximizes the sum of the evidence (or marginal likelihood) lower bound of each individual sample. The evidence lower bound for each sample can be written as <cit.>
  
  
      ℒ(θ,ϑ;𝐱^(i))=ℒ^x + ℒ^z
  
  
  with
  
  
      ℒ^x = 𝔼_q_ϑ(𝐳|𝐱^(i))[log(p_θ(𝐱^(i)|𝐳)]
  
  
  and
  
  
      ℒ^z=-D_KL(q_ϑ(𝐳|𝐱^(i))||p(𝐳))
  
  
  where 𝐳 refers to the codes or latent vectors, p_θ(𝐱|𝐳) is the (probabilistic) decoder, q_ϑ(𝐳|𝐱) is the (probabilistic) encoder, 𝔼 denotes the expectation operator, D_KL denotes the Kullback-Leibler distance and, θ and ϑ are the parameters (weights and biases) of the DNNs for the decoder and encoder, respectively.
  
  In order to maximize the evidence lower bound in Eq. (<ref>), its gradient with respect to both θ and ϑ is required, however, this is generally intractable and therefore an estimator is used.  This estimator is based on a so called reparameterization trick of the random variable 𝐳∼ q_ϑ(𝐳|𝐱) which uses an auxiliary noise ϵ. In the case of a VAE, the encoder is defined as a multivariate Gaussian with diagonal covariance:
  
  
      q_ϑ(𝐳|𝐱) = 𝒩(𝐡_ϑ(𝐱),𝐮_ϑ(𝐱) · I_d)
  
  
  where 𝐡_ϑ(𝐱) and log𝐮_ϑ(𝐱) are modeled with DNNs and I_d is a d × d diagonal matrix. Then, the encoder and the auxiliary noise ϵ are used in the following way during training <cit.>
  
  
      𝐳 = 𝐡_ϑ(𝐱)+𝐮_ϑ(𝐱) ⊙ϵ, ϵ∼ p(ϵ)
  
  
  where ⊙ denotes an element-wise product. Often Eq. (<ref>) has an analytical solution, then only Eq. (<ref>) is approximated with the estimator as <cit.>
  
  
      ℒ^x = 1/M∑_j=1^M log(p_θ(𝐱^(i)|𝐳^(i,j)))
  
  
  where 𝐳^(i,j) = 𝐡_ϑ(𝐱^(i))+𝐮_ϑ(𝐱^(i)) ⊙ϵ^(j), ϵ^(j)∼ p(ϵ) and M is the number of samples used for the estimator. Further, if we set the decoder p_θ(𝐱|𝐳) as a multivariate Gaussian with diagonal covariance structure, then
  
  
      p_θ(𝐱|𝐳) = 𝒩(𝐠_θ(𝐳),𝐯_θ(𝐳) · I_D)
  
  
  where 𝐠_θ(𝐳) and log𝐯_θ(𝐳) are modeled with DNNs and I_D is a D × D diagonal matrix. In this work, we consider only the mean of the decoder p_θ(𝐱|𝐳) which is just the (deterministic) generator 𝐠_θ(𝐳). Then, the corresponding (mean-square error) loss function may be written as
  
  
      ℒ^x = 1/M∑_j=1^M 𝐠_θ(𝐳^(i,j))-𝐱^(i)^2
  
  
  The described setting allows for the gradient to be computed with respect to both θ and ϑ and then stochastic gradient descent is used to maximize the lower bound in Eq. (<ref>). In the rest of this work, we drop the subindex θ in 𝐠(𝐳) to simplify notation and also because once the DGM is trained, the parameters θ do not change, i.e. they are fixed for the subsequent inversion.
  
  As previously mentioned, it is often possible to analytically integrate the Kullback-Leibler distance in Eq. (<ref>). In this work, we consider that p(𝐳) and q_ϑ(𝐳|𝐱) are both Gaussian therefore Eq. (<ref>) may be rewritten as <cit.>:
  
  
      ℒ^z=1/2∑_i=1^d (1+log((u_i)^2)-(h_i)^2-(u_i)^2)
  
  
  where the sum is done for the d output dimensions of the encoder.
  
  Note that the term in Eqs. (<ref>), (<ref>) and (<ref>) may be interpreted as a reconstruction term that causes the outputs of the encode-decode operation to look similar to the training samples, while the term in Eqs. (<ref>) and (<ref>) may be considered a regularization term that enforces the encoder q_ϑ(𝐳|𝐱) to be close to a prescribed distribution p(𝐳). In practice, one may add a weight to the second term <cit.> of the lower bound as:
  
  
      ℒ(θ,ϑ;𝐱^(i))=ℒ^x + βℒ^z
  
  
  to prevent samples to be encoded far from each other in the latent space, which may cause overfitting of the reconstruction term and degrade the VAE's generative performance. The overall process of training and generation for a VAE is depicted in Fig. <ref>.
  
  
  
  Note that in setting up the VAE one has to choose: (1) the architectures of the encoder and decoder, (2) the probability distribution p(𝐳), (3) the noise distribution p(ϵ) and (4) the regularization weight β. As mentioned in Section <ref>, these choices may impact the nonlinearity of the generator and its ability to preserve topology, which in turn affect the mapping of the data misfit function γ(𝐱) in latent space and possibly diminish the performance of inversion methods. In this work, we assume that an architecture for the generator (decoder) 𝐠(𝐳) is chosen so that it performs sufficiently good in terms of reproducing the patterns. For instance, when gridded spatial distributions (images) are considered, a typical choice is a deep-convolutional neural network <cit.>. While the choice of the probability distribution p(𝐳) may aid in obtaining a well-behaved generator, e.g. by selecting a probability distribution with the same (or similar) topology as the real manifold <cit.>, we expect such a choice to be highly problem (pattern) dependent. Therefore herein we focus on the other two possible controls: the noise distribution p(ϵ) and the regularization weight β. 
  
  The effect of the regularization weight β is such that when increased the encoded training samples tend to lie closer to the prescribed probability distribution p(𝐳). Then, one may picture the transformation of the encoder as taking the low-dimensional approximate manifold in the ambient space and charting it (e.g. by bending, stretching and even folding) into the region defined by p(𝐳) in the latent space and the generator as the transformation undoing such charting.
  
  While the effect of β in a VAE is relatively easy to understand, the effect of the noise distribution p(ϵ) is not so straightforward. First, note that the typical choice of a diagonal noise as p(ϵ) = 𝒩(0, α· I_d) where α denotes a constant variance (frequently set to α=1.0) and I_d is a d × d diagonal matrix is usually done for tractability or computational convenience <cit.>. However, it has been proposed recently that the choice of a diagonal noise has an impact on a property called disentanglement <cit.>. Such disentanglement basically means that different latent directions control different independent characteristics of the training (or generated) samples. They explain that a diagonal p(ϵ) might induce an encoding that preserves local orthogonality of the ambient space. In this work, we argue that the choice of a diagonal p(ϵ) (which is usually done only for computational convenience) might be useful in producing a well-behaved generator. 
  
  In order to visualize the joint effect of α and β, Fig. <ref> shows a synthetic example where samples in a two-dimensional ambient space lie close to a rotated "eight-shaped" manifold (Fig. <ref>a). In addition, to study the impact on inversion, a convex data misfit function γ(𝐱) in the same space (created synthetically with a negative isotropic Gaussian function) is shown in Fig. <ref>b. The latent space is also chosen two-dimensional for visualization purposes but recall that for a real case the dimensionality of the latent space is usually much lower than the one of the ambient space. Then, Fig. <ref> considers nine different combinations for the values of α and β to show how the (nonlinear) generator 𝐠(𝐳) maps a region of the latent space (denoted by the z-axes in the first three rows) into the ambient space (denoted by the x-axes in the last three rows) in order to approximate the manifold in Fig. <ref>a. To visualize the deformation caused by the generator, an orthogonal grid in the z-axes and its mapping into the x-axes (a deformed grid) are shown (both on the left of each inset). The corresponding encoded training samples are shown in red in the z-axes (left of each inset) and their reconstruction (resulting from the operation of encode-decode) is shown also in red in the the x-axes (right of each inset), where also the original training samples are shown (in blue) to assess the accuracy of reconstruction. Samples obtained from a Gaussian distribution with a unitary diagonal covariance p(𝐳) are shown in the z-axes in orange (left of each inset), while their generator-mapped values are shown also in orange in the x-axes (right of each inset). Finally, the mapping of the data misfit function in Fig. <ref>b into the latent space is shown in the z-axes (right of each inset).
  
  
  
  
  
  It is worth mentioning a few effects visible in the illustrative example of Figs. <ref> and <ref>. First, note that increasing the variance of α seems to cause the grid to be more "rigid" locally (grid lines tend to intersect more at right angles) while going through the generator which may in turn help in preserving topology and controlling non-linearity (e.g. compare the deformation of the grids for different values of α for β=0.01), and more importantly, in preserving the convexity of the data misfit function in the latent space (the mapped misfit function using α=0.1 and β=0.01 has a single global minimum, while the misfit function for α=0.01 and β=0.01 has two minima in latent space). Also note that both α and β should be set in order to not cause a significant degradation in: (1) the reconstruction of the patterns, e.g. the cases of α=1.0 with both β=0.1 and β=0.01 show that the "eight-shape" is not completely reconstructed (seen in red samples not fully overlaying the blue samples in x-axes), or (2) the similarity of the encoded samples to the prescribed distribution p(𝐳), e.g. the case of α=0.01 and β=0.1 shows that encoded samples (in red) are too concentrated (lower variance) and therefore far from the prescribed normal distribution with unit variance. In this case, the intermediate values (α=0.1 and β=0.01) seem to provide the best choice in terms of reconstruction of the patterns, generative accuracy and convexity of the misfit function in latent space.
  
  In summary, a generator 𝐠(𝐳) that preserves topology and contains nonlinearity is the best choice for gradient-based inversion in the latent space because it preserves convexity of the objective function. Note, however, that if the topology of the probability distribution p(𝐳) is different to the one of the real manifold ℳ, this strategy may result in approximate manifolds ℳ' that do not account for all topological differences—e.g. that partially cover holes of the real one (see e.g. Fig. <ref>b)—and therefore might produce models that have non-accurate patterns when sampling from p(𝐳). We argue that the two training parameters α and β of a VAE may be chosen in order for the latter issue to not be severe, i.e. the generated patterns do not deviate too much from the training patterns, while still approximately preserving convexity of the objective function in the latent space.
  
  To test our proposed method we implement a VAE in PyTorch <cit.> and use training samples cropped from a "training image" which is large enough to have many repetitions of the patterns at the cropping size—a requirement similar in MPS. For our synthetic case, we use the training image of 2500 × 2500 pixels from <cit.> and the cropping size is chosen to fit the setting of our synthetic experiment (explained in detail in Sec. <ref>). Fig. <ref>a shows a patch of the training image and the position of the three (cropped) training samples shown Fig. <ref>b. Three generated samples from our proposed VAE trained with such croppings are shown Fig. <ref>c. For comparison, Fig. <ref>d shows three samples generated with the SGAN proposed by <cit.>. Patterns of generated samples in Fig. <ref>c are not completely accurate comparing to those of the training image or the SGAN—they might display e.g. some breaking channels and smoothed edges. As mentioned above, this is expected for our proposed VAE because the approximate manifold fills some holes of the real manifold and may have less curvature. However, we argue that such inaccuracies may not cause significant error or bias while performing inversion in practice because an informative dataset will generally make the inversion land in appropriate models (given the prescribed patterns were selected correctly). More importantly, in contrast to the SGAN, a modified gradient-based inversion (such as that presented in Sec. <ref>) will generally find a consistent minimum when applied with our proposed VAE regardless of the initial model.
  
  
  
  
  
   §.§ Stochastic gradient descent with decreasing step size
  
  
  
  Note that even when topology is preserved and nonlinearity is contained, the data misfit function in the latent space might still present some local minima. Using our proposed VAE approach in the synthetic case study, the resulting misfit function seems to have the shape of a global basin of attraction with some local minima of less amplitude. To deal with such remaining local minima we propose to use a stochastic gradient descent (SGD) method instead of regular gradient-based optimization.
  
  SGD methods are commonly used in training machine learning models to cope with large datasets <cit.> and it has also been shown they are able to find minima that are useful in terms of generalization <cit.>. They essentially use an estimator for the gradient of the objective function computed only with a batch of the data. Such estimator is used in each gradient descent iteration and may be written for the case of inversion in the latent space as:
  
  
      𝐳_k+1 = 𝐳_k - ℓ·∇_𝐳ζ(𝐳)_k
  
  
  where k denotes the iteration index, ℓ is the step size (or learning rate) and the gradient estimator ∇_𝐳ζ(𝐳)_k is computed by using Eq. (<ref>) for a data batch (i.e. a subset of 𝐲) which is different for each k-th iteration but of constant size b. Relying on such estimator makes SGD methods less likely to get trapped in local minima if a sufficiently large step is chosen. However, if such a step is too large the optimization will have issues when it is close to the global minimum, usually seen in the form of high misfit and oscillations in the value of the objective function. Note that similar to other stochastic optimization methods, SGD only guarantees convergence to the global minimum with a certain probability, however if modified in the right way for the type of problems to be solved and its parameters chosen appropriately such probability could be very close to one.
  
  Recently, it has been proposed that using SGD may be seen as optimizing a smoothed version of the objective function obtained by convolving it with the gradient "noise" resulting from batching <cit.>. The degree of noise (and therefore the degree of smoothness) is controlled by the ratio of the learning rate to the batch size ℓ / b <cit.>. Therefore if we choose to decrease the value of ℓ (while keeping b constant) as the optimization progresses we might be able to achieve lower misfit values i.e. get sufficiently close the global minimum. This may be implemented by using:
  
  
      ℓ_k+1 = c_ℓ·ℓ_k
  
  
  where a constant value of c_ℓ < 1.0 and a starting value ℓ_0 must be chosen. In practice, the method may be further improved by also decreasing the controlling (regularization) term in Eqs. (<ref>) and (<ref>) in order to prevent that large initial steps diverge from the region of the latent space where the manifold is defined <cit.>. Then, similarly to ℓ this may be done as:
  
  
      λ_k+1 = c_λ·λ_k
  
  
  again a constant c_λ < 1.0 and a starting value λ_0 must be selected.
  
  The combined effect of simultaneously decreasing ℓ and λ is illustrated in Fig. (<ref>) for a simple synthetic problem in a two-dimensional (d=2) latent space ℛ^d. The misfit term (i.e. first term of Eq. (<ref>)) of the synthetic problem is shown in Fig. <ref>a. Assuming that p(𝐳) is a normal distribution 𝒩(0,I_d) where I_d is a d × d identity matrix, we propose a specific regularization term R(𝐳) that will preferentially stay in the regions of higher mass (where most samples are located). This is done by radially constraining the search space by means of a χ-distribution, i.e. the regularization term is written as:
  
  
      R(𝐳)= (𝐳-μ_χ)^2
  
  
  where μ_χ is the mean for a χ-distribution with d degrees of freedom. Dashed lines in Fig. <ref>a denote this mean together with the 16- and 84-th percentiles. In general, this is especially useful for higher dimensionalities where most of the mass of a normal distribution is far from its center <cit.>. Then, Eq. (<ref>) may be rewritten as:
  
  
      ζ(𝐳) = 𝐟(𝐠(𝐳))-𝐲^2 + λ (𝐳-μ_χ)^2
  
  
  and correspondingly Eq. (<ref>) may be expressed as:
  
  
      ∇_𝐳ζ(𝐳) =  𝐉(𝐳)^T ∇_𝐱γ(𝐱) + 2λ𝐳(1-μ_χ/𝐳)
  
  
  As mentioned above, this gradient is often computed simply by adding a layer to the autodifferentiation of the generator. One optimization instance for a random initial model is shown in Fig. <ref>b, while the behavior of the misfit and 𝐳 is shown in Fig. <ref>c,d. Notice the rather "noisy" inversion trajectory, but also its ability to escape local minima. The effect of decreasing ℓ is seen in Fig. <ref>c by the decreasing of the oscillations amplitude as the optimization progresses, while the effect of decreasing λ is noticeable in Fig. <ref>d by the progressive shifting of 𝐳 away from μ_χ.
  
  
  
  The strategy described above and stated by Eq. (<ref>) is generally applicable to DGMs that use an independent normal distribution as its probability distribution p(𝐳) and whose generator is well-behaved. In this work, we consider a VAE whose training parameters β and p(ϵ) are chosen so that it results in a mildly nonlinear inversion for which such SGD strategy is generally useful.
  
  
  
   §.§ Inverse problem: traveltime tomography
  
  
  
  To test our proposed method and compare it with a previous instance of inversion with a DGM, we consider an identical setting to that used in <cit.>. Such setting considers a dataset of borehole ground penetrating radar (GPR) traveltime tomography. To obtain a subsurface model 𝐱∈ℝ^D this method relies in contrasts of electromagnetic wave velocity which is related to moisture content and therefore to porosity for saturated media. The tomographic array considers a transmitter antenna in one borehole and a receiver antenna in the other, each of which is moved to different positions and a vector of measurements 𝐲∈ℝ^Q is obtained by taking the traveltime of the wave's first arrival for each transmitter-receiver combination. We assume that the sensed physical domain is a 6.5 × 12.9 m plane (i.e. the two-dimensional region between the boreholes) and is discretized in 0.1 × 0.1 m cells of constant velocity to represent spatial heterogeneity (i.e. a representation of D= 65 × 129 = 8385 cells is obtained). We consider a binary subsurface (e.g. composed of two materials with different porosity) with respective wave velocities of 0.06 and 0.08 m ns-1. Measurements are taken every 0.5 m in depth (the first being at 0.5 m and the last at 12.5 m) resulting in a dataset of Q= 625 traveltimes. For one instance of our synthetic case, we add normal independent noise η∼𝒩(0,σ^2 · I_Q) where σ^2 is the noise variance and I_Q is a 625 × 625 diagonal matrix.
  
  Similarly to <cit.>, we first consider a fully linear forward operator 𝐟 for which raypaths are always straight, i.e. independent of the velocity spatial distribution. For this case Eq. (<ref>) may be rewritten as:
  
  
      𝐲 = 𝐅𝐱 + η
  
  
  where 𝐅 is a matrix of dimension Q × D in which a certain row contains the length of the raypath in each cell of the model for a certain transmitter-receiver combination. The corresponding gradient of the misfit ∇_𝐱γ(𝐱) to be used in Eq. (<ref>) for the solution of the inversion is:
  
  
      ∇_𝐱γ(𝐱) = -2𝐅^T (𝐲-𝐅𝐱)
  
  
  We also consider the case of a more physically realistic nonlinear forward operator 𝐟 (see Eq. (<ref>)) for which raypaths are not straight. In particular, we consider a shortest path (graph) method which uses secondary nodes to improve the accuracy of the simulated traveltimes as proposed by <cit.> and implemented in PyGIMLi <cit.>. For this case, when inversion with Eq. (<ref>) is pursued, we linearize the forward operator 𝐟 in order to compute the gradient:
  
  
      ∇_𝐱γ(𝐱) = -𝐒(𝐱)^T(𝐲-𝐟(𝐱))
  
  
  where is 𝐒(𝐱) is the Q × D jacobian matrix of the forward operator whose elements are:
  
  
      [𝐒(𝐱)]_i,j = ∂ f_i(𝐱)/∂ x_j
  
  
  The elements of the jacobian 𝐒(𝐱) are computed by the shortest path method and also represent lengths of raypaths. In contrast to the linear case, these have to be recomputed in every iteration. Both the nonlinear forward operator and the need for recomputing the jacobian result in higher computational cost compared to the linear operator.
  
  The method proposed in Sec. <ref> to perform gradient-based inversion with a VAE should work for the linear forward operator because the nonlinearity in the inverse problem arises only due to the generator 𝐠(𝐳) which is moderate when the latter is well-behaved. However, since the considered nonlinear forward operator in Eq. (<ref>) is only mildly nonlinear (when contrast in velocities is not extreme), the same method may also provide good inversion results for this operator.
  
  
  
  § RESULTS
  
  
  
  
  
   §.§ Training of VAE
  
  
  As previously mentioned, our proposed method relies on a VAE whose training parameters are selected in order to improve gradient-based inversion. The training samples are the croppings detailed in Sec. <ref> whose dimensionality is D= 8325 and we consider a latent code of dimensionality d= 20. The probability distribution p(𝐳) is an independent multinormal distribution 𝒩(0,I_d) with I_d an identity matrix of size 20 × 20. The architecture of the encoder and the decoder includes 4 convolutional layers, 2 fully-connected layers and instance normalization is used between each layer (details may be consulted in the associated code). The training parameters relevant to our proposed  method are chosen in the following way: (1) β in Eq. (<ref>) is given a value of 1000, chosen by visually assessing the generated samples which also coincided with a moderate visual deviation from the prescribed p(𝐳); (2) the distribution p(ϵ) in Eq. (<ref>) is given the typical value of a diagonal unit variance (α=1.0), which enforces local orthogonality while passing through the generator and therefore aids in preserving topology and controlling non-linearity. The value of β is rather high compared to previous studies, e.g. <cit.> used a value of 20 for similar two-dimensional patterns, but seems to cause slightly higher compression ratios (in our work the compression ratio is 420 compared to 200 in the mentioned study). The VAE is trained by maximizing the lower bound in Eq. (<ref>) using 105 iterations and batches of 100 random croppings in each iteration (a GPU was used in order to reduce training time). In the following, we test the performance of this VAE when used for our proposed SGD-based inversion with a linear and a mildly nonlinear forward operator.
  
  
  
   §.§ Case with a linear forward model
  
  
  In this section, we consider the linear operator in Eq. (<ref>) and assess the performance of our proposed DGM inversion approach: using a VAE trained as above (to have a well-behaved generator) and SGD with both decreasing step size and regularization to optimize Eq. (<ref>). We aim to show that such approach is robust regarding its convergence to the global minimum and therefore assess its performance by using 100 different initial models. Compared to previous studies, our proposed approach involves changes in both the DGM and the optimization, therefore we compare with the base cases listed in Table <ref> to show the impact of each change proposed. As denoted by the columns of this table, the different cases consider: (1) VAE and SGAN as DGMs, (2) SGD and Adam <cit.> as stochastic optimizers, (3) data batching for computing the gradient ∇_𝐳ζ(𝐳), which basically means using SGD when batching and using (regular) gradient-descent when not batching, (4) regularization in the latent space, with "origin" being the one proposed in <cit.> and "ring" the one proposed herein, and (5) decreasing of the step size (or learning rate). Our proposed approach is then labeled as "VSbrd". We also show the chosen values for the step size ℓ and its decreasing factor c_ℓ when applicable—for these cases the values of λ = 10.0 and c_λ = 0.999 are used. The number of iterations for inversion is set to 3000 for all cases. When data batching is used, the batch size b is 25 (of a total of 625) and is sampled with no replacement, then the whole dataset is used every 25 iterations (i.e. the number of epochs is 120 for a total of 120 × 25 = 3000 iterations).  Note that we compare against the approach in <cit.>, where SGAN is used as DGM and Adam (gradient-descent with adaptive moments) are used to optimize the resulting objective function—this case is labeled "SAnnn" in Table <ref>. We also consider the case where we apply our proposed SGD to the same SGAN (labeled as "SSbnd"). For both of these cases instead of regularization we use stochastic clipping in the latent space (Laloy et al. 2018, 2019) because a uniform p(𝐳) with finite support is used.
  
  
  
  We consider 6 different truth subsurface models to assess our method and compare with the base cases: (1) a set of three models cropped directly from the training image and (2) a set of three models obtained by generating from the trained VAE. Both sets include models with three different degrees of complexity. These truth models are shown in the first row of Fig. <ref> where "mc" refers to the first set, "mv" refers to the second set and the degree of complexity is denoted by a subscript, where "1" denotes least complex and "3" most complex. For the first set (mc), to avoid "memorizing" the croppings we exclude them (and any overlapping cropping) from the samples used to trained the VAE. The second set (mv) is similar to the one used by <cit.> to test the performance of their setup, only in their case the models were generated from a SGAN instead of a VAE. For each one of these truths, we generate synthetic data applying the forward operator 𝐅 and use these data to perform gradient-based inversion for each case in Table <ref>. 
  
  
  
  We first consider no added noise to the synthetic dataset, hence after inversion the data misfit should be close to zero for inverted models that are sufficiently close to the global minimum. To define a threshold for this data misfit beyond which inverted models are "accepted", we use the RMSE between these synthetic data and data obtained by applying the forward operator on models resulting from passing the truth models through a VAE's encoding-decoding (these models are shown in the second row of Fig. <ref> and the corresponding values for the threshold are shown in Table <ref>). This is done because we found the encode-decode reconstructed models to be visually very similar to the truth models (compare first and second rows of Fig. <ref>) and also show a low model RMSE when compared to them (computed just as the difference of pixel values between truth model and the encode-decode model and shown Table <ref>). Once such threshold is defined for each truth model, gradient-based inversion is run for the same 100 initial models for all cases in Table <ref>. Note that no convergence criteria were set in order to compare to all base cases (some cases such as "SAnnn" do not allow for easily defining such criteria) but in practice it is possible to set them for our proposed approach (VSbrd) in terms of a minimal change in either step size and/or data misfit. This also means that for some cases (including our proposed VSbrd) the 3000 iterations may not be necessary for all truths and all initial models. Results for the number of accepted inverted models are shown in Table <ref> while the corresponding mean of the misfit (expressed as RMSE) for the 100 inversions is shown in Table <ref>.
  
  
  
  
  
  
  
  As seen in Table <ref>, given our defined threshold: (1) the cases where VAE and SGD with decreasing step were used (VSbod and VSbrd) resulted in all inverted models being accepted, (2) the cases where SGAN was used (SAnnn and SSbnd) resulted in almost all models being rejected (only two models accepted for SAnnn with truth mv1), and (3) the cases where VAE and non-decreasing step size SGD was used (VSnnn and VSbnn) resulted in some inverted models being accepted. Note also that using SGD (data batching) without a decreasing step size (VSbnn) results in less accepted models compared to GD (VSnnn), highlighting the importance of our proposed decreasing step size and regularization. As shown in Table <ref> a higher mean RMSE is related to a lower number of accepted models. Furthermore, two things worth highlighting in Table <ref> are (1) the general improvement for inversion with SGAN caused by our proposed SGD compared to Adam for most truth models (compare SSbnd and SAnnn), which means that our proposed SGD is advantageous regardless of the DGM, and (2) the slight improvement caused by our proposed regularization compared to the one from <cit.>.
  
  Examples of inverted models obtained for the different cases in Table <ref> using the cropped truth with moderate complexity (mc2) are shown in Fig. <ref>. Here, truth models are shown in Fig. <ref>a while Fig. <ref>b shows one example of an accepted model for cases that have at least one (VSnnn, VSbnn, VSbod and VSbrd). Similarly, Fig. <ref>c shows one example of a rejected model for applicable cases (VSnnn, VSbnn, SAnnn and SSbnd). Finally, the corresponding data RMSE vs. iteration number plots are shown in Fig. <ref>d (in blue for accepted models and red for rejected ones) and corresponding model RMSE plots are shown in Fig. <ref>e. Note both the higher similarity with the truth model (i.e. note the low model RMSE and compare models in Fig. <ref>b,c with those in Fig. <ref>a) and the lower RMSE for accepted models. Also, examples of inverted models for our proposed approach (VSbrd) using all the truths are shown in Fig. <ref>b, together with plots of RMSE vs. iteration number (Fig. <ref>d) and norm of 𝐳 vs. iteration number (Fig. <ref>e). For cropped truths (mc) it seems that visual similarity decreases and final data RMSE of inverted models increases as complexity increases, whereas for generated truths they seem independent of complexity. Notice the overshoot in 𝐳 in the initial iterations and its eventual convergence close to μ_χ as defined in Eq. (<ref>).
  
  
  
  
  
  To study the effect of noise for our proposed approach (VSbrd), we added noise with a standard deviation σ = 0.25 ns to the synthetic traveltime data. Corresponding results are shown in the rightmost column of Tables <ref> and <ref> and in Fig. <ref>c (with corresponding data RMSE and 𝐳 norm plots in Fig. <ref>d,e). The threshold in this case is set equal to the one for the noise-free case plus σ and when using it all inverted models with our proposed approach are accepted. It is also worth noticing the relative robustness of the method to noise, as shown by the corresponding mean misfit values in Table <ref> that indicate no significant overfitting, i.e. the mean misfit values are close to the noise-free threshold plus σ even if no traditional regularization was used. The latter means that optimizing in the latent space of the DGM is effectively constraining the inverted models to display the prescribed patterns.
  
  
  
   §.§ Case with a nonlinear forward model
  
  
  After showing that our proposed method works with the linear forward operator for the synthetic case considered, we now test its performance with a nonlinear forward operator. For inversion, the general form of Eq. (<ref>) is used and the gradient in the latent space given in Eq. (<ref>) is computed using Eq. (<ref>). As mentioned in Sec. <ref>, we consider a shortest path method to solve for the traveltime for which we use 3 secondary nodes added to the edges of the velocity grid. Note that the jacobian 𝐒(𝐱) in Eq. (<ref>) has to be recomputed at every iteration. Given the higher computational demand for inversion with the nonlinear forward operator and since it was already shown to be the best performing approach for the linear forward operator, we only test our proposed approach VSbrd with all the truths and for a single initial model (Fig. <ref>). This was done both without noise and with noise added using the same standard deviation σ = 0.25 as in the linear operator scenario. We select the following values for the required inversion parameters: ℓ= 0.1, c_ℓ= 0.8, λ= 1.0 and c_λ= 0.99. The total number of iterations is 750 with data batching of size 25 similar to the linear case. Note that to further reduce the number of iterations required for inversion we use a lower c_ℓ compared to the linear case, but the decreasing in Eq. (<ref>) is only done every 5 iterations. This may cause the method to converge to the global minimum with lower probability, however it seems to still be high enough since all of the inversions with no added noise are very similar to the truth models. Also, using the threshold obtained by encoding-decoding the truth models (now computed with the nonlinear forward operator) all inverted models are accepted (these models are shown in Fig. <ref>b). When considering added noise, results are similar but inversion seems to converge to the global minimum with slightly lower probability (6 out of 8 inversions are accepted) and accepted models are shown in in Fig. <ref>c. The behavior of the misfit during optimization (Fig. <ref>d) is similar to the linear case, although oscillations of a slightly higher amplitude are still visible in the last iterations (mainly due to the lower number of iterations). To partially solve the latter issue, we take as inverted model the model with lowest misfit and not the one for the final iteration (these are the models shown in Fig. <ref>b,c). The plot of the norm of 𝐳 vs. iterations in Fig. <ref>e shows a similar behavior to the linear case, although there seems to be more oscillations in 𝐳 during initial iterations.
  
  
  
  
  
  § DISCUSSION
  
  
  
  When training the VAE for our considered synthetic case, only the selection of β was done while the variance of p(ϵ) was not changed (a unity variance α=1.0 was used). However, as noted in Sec. <ref> and in Fig. <ref>, the variance α also has an impact on (gradient-based) inversion as it might be used to control the generator's nonlinearity and its induced changes in topology.
  
  In order to select SGD parameters in our proposed approach, we suggest looking jointly at the behavior of the misfit and norm of 𝐳. For instance, if a certain number of iterations is desired for computational reasons, we suggest choosing first ℓ and c_ℓ that produce a behavior of the misfit similar to that in Fig. <ref>d, i.e. oscillations of high amplitude at the beginning and then progressive attenuation of the oscillations in such a way that at the end they are negligible. Note, however, that inversion may have to be run a few times because divergence may occur during initial iterations (this is easily seen in the value of 𝐳 taking values far from μ_χ). Once ℓ and c_ℓ are chosen, the selection of λ and c_λ is done only to prevent divergence, this may be achieved by looking for a behavior similar to that in Fig. <ref>e. An initial overshoot in 𝐳 is normal (and even necessary) since the method is exploring more rapidly the latent space, however, it should eventually converge to a value close to μ_χ.
  
  The results for gradient-based inversion using our proposed approach point to a (possible) conflict between the accuracy of the reproduced patterns and the feasibility of gradient-based inversion with DGMs. As mentioned above, this is due to a non-convex objective function in latent space resulting from the generator's nonlinearity and its induced changes in topology. In this work, we argue that nonlinearity and changes in topology might be safely controlled by selecting certain values of α and β while training a VAE in order to improve performance of gradient-based inversion. We empirically prove the validity of this statement for our case study with specific values α and β. In general (for inversion with DGMs), this implies that a tradeoff between generative accuracy and a well-behaved generator may be found. The latter statement also supports our assumption regarding the "holes" of the real manifold for the case of channel patterns (as mentioned in Section <ref>): when approximating the real manifold using a VAE with a well-behaved generator, the approximate manifold will tend to fill the holes and therefore produce breaking channels. While the generator's nonlinearity was already identified by <cit.> as a potential factor for hindering gradient-based inversion, its causes (curvature and topology of the real manifold) and the possible induced changes in topology have not been previously explained as factors in degrading the performance of gradient-based inversion in the latent space (to the authors' knowledge).
  
  In general, good performance of DNNs for some tasks is usually associated with their ability to change topology <cit.>. However, when one wants to use the latent variables or codes of DGMs for further tasks and not just for generation, these changes in topology might become an issue. For instance, we interpret the misfit "jumps" seen in gradient-based inversion with SGAN (as seen in Fig. <ref>c for case SAnnn) as resulting from the "gluing" or "collapsing" in latent space of holes in the real manifold—either caused by an induced change in topology or a high nonlinearity in the SGAN generator. Some studies have even suggested that if one wants to obtain useful geometric interpretations in the latent space (e.g. to perform interpolation), the activation functions should be restricted to ones that are smooth <cit.>, that means e.g. not using the ReLU activation function that is generally recognized to result in faster learning. In contrast, in this work we do consider ReLU activation functions but control the changes in topology by means of a combination of α and β, whether this might nullify the advantages of ReLU is still an open question. Note  however that, in general, control of induced changes in topology and high nonlinearities (as in our proposed approach) might be useful for any inversion method that relies in the concept of a neighborhood (e.g. MCMC and ensemble smoothers).
  
  Besides its good performance for gradient-based inversion, a further advantage of our approach when compared to the previous approaches is that when the data used for inversion is not sufficiently informative, regularization in the latent space might be used to constrain to the most common patterns with our regularization term in Eq. (<ref>). This statement provides an interesting paradigm where regularization in latent space might be seen as a flexible way to incorporate complex regularization. In contrast, a disadvantage of our proposed approach is that GANs in general result in higher generative accuracy (all generated patterns look more similar to those in the training image), however, as previously mentioned this may be in conflict with performance of certain inversion methods. Also, as may be noticed in the relation between the data misfit and the degree of complexity for cropped truths, a limiting factor in using our VAE is its inability to produce new highly complex patterns. However, this lack of innovation (or sample diversity) is generally present in other methods and may be even more severe for regular GANs, where the phenomenon is known as mode collapse. Recently, different ways to control such mode collapse in VAEs and GANs have been proposed <cit.>.
  
  Finally, regarding our proposed SGD we must note that similar results might be obtained with a MCMC method where information about the gradient is taken into account. For example, <cit.> use a Metropolis-adjusted Langevin method which basically follows a gradient-descent and adds some noise to the step. However, the noise added to the gradient step in our approach is different—SGD noise has been shown to be approximately constant but anisotropic <cit.>. Another possible alternative to our method is to use Riemannian optimization, which is possible when the DGM approximate manifold is smooth. Although it is possible to compute the direction of the gradient by using the pullback Riemannian metric, which may be obtained as suggested by e.g. <cit.>, it is not straightforward to compute the step because it would have to be along a geodesic curve instead of a straight path and such geodesics are computationally demanding to compute.
  
  
  
  § CONCLUSIONS
  
  
  
  In this work the principal difficulties of performing inversion with deep generative models (DGMs) are reviewed and a conflict between generated pattern accuracy and feasibility of gradient-based inversion is identified. Also, an approach based on a variational autoencoder (VAE) as DGM and a modified stochastic gradient descent method for optimization is proposed to address such conflict. We show that two training parameters of the VAE (the weight factor β and the variance α of the encoder's noise distribution p(ϵ)) may be chosen in order to obtain a well-behaved generator 𝐠(𝐳), i.e. one that is mildly nonlinear and approximately preserves topology when mapping from latent space to ambient space. This helps in maintaining the convexity of the misfit function in the latent space and therefore improves the behavior of gradient-based inversion. We highlight changes in topology which have not been previously identified as impacting the convexity of the inversion objective function. In contrast to prior studies where gradient-based inversion was used, our approach converges to the neighborhood of the global minimum with very high probability for both a linear forward operator and a mildly nonlinear forward operator with and without noise. We argue that when using DGMs in inversion, a tradeoff may be found where inverted models are close enough to the prescribed patterns while low cost gradient-based inversion is still applicable—our proposed approach relies on such tradeoff and produces inverted models with significant similarity to the training patterns and a low data misfit.
  
  
  
  § ACKNOWLEDGEMENTS
  
  
  This work has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement number 722028 (ENIGMA ITN).
  
  
  
  § COMPUTER CODE AVAILABILITY
  
  
  All code necessary to reproduce the test case is available at: https://github.com/jlalvis/VAE_SGDhttps://github.com/jlalvis/VAE_SGD.
  
  
  
  
  
  
  
  
  
  
  