# Seismic Waveform Inversion by Stochastic Optimization

Authors: Tristan van Leeuwen; Aleksandr Y. Aravkin; Felix J. Herrmann
Journal: International Journal of Geophysics
 (2011)

---

## Abstract

We explore the use of stochastic optimization methods for seismic waveform inversion. The basic principle of such methods is to randomly draw a batch of realizations of a given misfit function and goes back to the 1950s.  The ultimate goal of such an approach is to dramatically reduce the computational cost involved in evaluating the misfit.  Following earlier work, we introduce the stochasticity in  waveform inversion problem in a rigorous way via a technique calledrandomized trace estimation. We then review theoretical results that underlie recent developments in the use of stochastic methods for waveform inversion. We present numerical experiments to illustrate the behavior of different types of stochastic optimization methods and investigate the sensitivity to the batch size and the noise level in the data. We find that it is possible to reproduce results that are qualitatively similar to the solution of the full problem with modest batch sizes, even on noisy data. Each iteration of the corresponding stochastic methods requires an order of magnitude fewer PDE solves than a comparable deterministic method applied to the full problem, which may lead to an order of magnitude speedup for waveform inversion in practice.


---

## Body

## 1. Introduction

The use of simultaneous source data in seismic imaging has a long history. So far, simultaneous sources have been used to increase the efficiency of data acquisition [1, 2], migration [3, 4], and simulation [5–7]. Recently, the use of simultaneous source encoding has found its way into waveform inversion. Two key factors play a role in this development: (i) in 3D, one is forced to use modeling engines whose cost is proportional to the number of shots (as opposed to 2D frequency-domain methods where one can reuse the LU factorization to cheaply model any number of shots) and (ii) the curse of dimensionality: the number of shots and the number of gridpoints grows by an order of magnitude.The basic idea of replacing single-shot data by randomly combined “super shots” is intuitively pleasing and has lead to several algorithms [8–11]. All of these aim at reducing the computational costs of full waveform inversion by reducing the number of PDE solves (i.e., the number of simulations). This reduction comes at the cost of introducing random crosstalk between the shots into the problem. It was observed by Krebs et al. [8] that it is beneficial to recombine the shots at every iteration to suppress the random crosstalk and that the approach might be more sensitive to noise in the data. In this paper, we follow Haber et al. [12] and introduce randomized source encoding through a technique called randomized trace estimation [13, 14]. The goal of this technique is to estimate the trace of a matrix efficiently by sampling its action on a small number of randomly chosen vectors. The traditional least-squares optimization problem can now be recast as a stochastic optimization problem. Theoretical developments in this area go back to 1950s, and we review them in this paper. In particular, we discuss two distinct approaches to stochastic optimization. The stochastic approximation (SA) approach consists of a family of algorithms that use a different randomization in each iteration. This idea justifies a key part of the approach described in Krebs et al. [8]. Notably, the idea of averaging the updates over the past is important in this context to suppress the random crosstalk; lack of averaging over the past likely explains the noise sensitivity reported by Krebs et al. [8]. The theory we treat here concerns only first-order optimization methods, though there has been a recent effort to extend similar ideas to methods that exploit curvature information [15].Another approach, called thesample average approximation (SAA), replaces the stochastic optimization problem by an ensemble average over a set of randomizations. The ensemble size should be big enough to suppress the crosstalk. The resulting problem may be treated as a deterministic optimization problem; in particular, one may use any optimization method to solve it.Most theoretical results in SA and SAA assume that the objective function is convex, which is not the case for seismic waveform inversion. However, in practice one starts from a “reasonable” initial model, and we may be able to converge to the closest local minimum. One would expect SA and SAA to be applicable in the same framework. Understanding the theory behind SA and SAA is then very useful in algorithm design, even though the theoretical guarantees derived under the convexity assumption need not apply.As mentioned before, the gain in computational efficiency comes at the cost of introducing random crosstalk between the shots into the problem. Also, the influence of noise in the data may be amplified by randomly combining shots. We can reduce the influence of these two types of noise by increasing the batch size, recombining the shots at every iteration, and averaging over past iterations. We present a detailed numerical study to investigate how these different techniques affect the recovery.The paper is organized as follows. First, we introduce randomized trace estimation in order to cast the canonical waveform inversion problem as a stochastic optimization problem. We describe briefly how SA and SAA can be applied to solve the waveform inversion problem. In Section3, we review relevant theory for these approaches from the field of stochastic optimization. The corresponding algorithms are presented in Section 4. Numerical results on a subset of the Marmousi model are presented in Section 5 to illustrate the characteristics of the different stochastic optimization approaches. Finally, we discuss the results and present the conclusions.


## 2. Waveform Inversion and Trace Estimation

The canonical waveform inversion problem is to find the medium parameters for which the modeled data matches the recorded data in a least-squares sense [16]. We consider the simplest case of constant-density acoustics and model the data in the frequency domain by solving
(1)H[m]u=q,
where H[m] is the discretized Helmholtz operator [ω2m+∇2] for the squared slowness m (with appropriate boundary conditions), u is the discretized wavefield, and q is the discretized source function; both are column vectors. The data are then given by sampling the wavefield at the receiver locations: d=Pu. Note that all the quantities are monochromatic. We hide the dependence on frequency for notational simplicity.We denote the corresponding optimization problem as(2)min⁡m⁡ϕ(m,Q,D)=∑ω‖PH[m]-1Q-D‖F2,
where D=[d1,d2,…,dN] is a frequency slice of the recorded data, and Q=[q1,q2,…,qN] are the corresponding source functions. Note that the dependence of H on ω has been suppressed. ∥·∥F denotes the Frobenius norm, which is defined as ∥A∥F=trace(ATA) (here ·T denotes the complex-conjugate transpose. We will use the same notation for the transpose in case the quantity is real). Note that we assume a fixed-spread acquisition where each receiver sees all the sources.In practice,H-1 is never computed explicitly but involves either an LU decomposition (cf., [17–19]) or an iterative solution strategy (cf., [20, 21]). In the worst case, the matrix has to be inverted separately for each frequency and source position. For 3D full waveform inversion, both the costs for inverting the matrix and the number of sources increase by an order of magnitude. Recently, several authors have proposed to reduce the computational cost by randomly combining sources [8–12].We follow Haber et al. [12] and introduce this encoding in a rigorous manner by using a technique called randomized trace estimation. This technique was introduced by Hutchinson [13] as a technique to efficiently estimate the trace of an implicit matrix. Some recent developments and error estimates can be found in Avron and Toledo [14].This technique is based on the identity(3)trace(ATA)=Ew(wTATAw)=lim⁡K→∞⁡1K∑k=1KwkTATAwk,
where Ew denotes the expectation over w. The random vectors w are chosen such that Ew(wwT)=I (the identity matrix). The identity can be derived easily by using the cyclic permutation rule for the trace (i.e., trace(ABC)=trace(CAB)), the linearity of the expectation, and the aforementioned property of w. At the end of the section, we discuss different choices of the random vectors w. First, we discuss how randomized trace estimation affects the waveform inversion problem.Using the definition of∥A∥F, we have
(4)ϕ(m,Q,D)=Ewϕ(m,Qw,Dw).
This reformulation of (2) is a stochastic optimization problem. We now briefly outline approaches to solve such optimization problems.


### 2.1. Sample Average Approximation

A natural approach to take is to replace the expectation overw by an ensemble average
(5)ϕK(m)=1K∑k=1Kϕ(m,Qwk,Dwk).
This is often referred to in the literature as the sample average approximation (SAA). The random crosstalk can be controlled by picking a “large enough” batch size. As long as the required batch size is smaller than the actual number of sources, we reduce the computational complexity.For a fixedm, it is known that the error |ϕ-ϕK| is of order 1/K (cf., [14]). However, it is not the value of the misfit that we are trying to approximate, but the minimizer. Unfortunately, the difference between the minimizers of ϕ and ϕK is not readily analyzed. Instead, we perform a small numerical experiment to get some idea of the performance of the SAA approach for waveform inversion.We investigate the misfit along the direction of the negative gradientgk (defined below)
(6)fK(α)=ϕK(m-αgK).
The data are generated for the model depicted in Figure 1(a), for 61 colocated, equidistributed sources and receivers along a straight line at 10 m depth and 7 randomly chosen frequencies between 5 and 30 Hz. The source signature is a Ricker wavelet with a peak frequency of 10 Hz. We use a 9-point discretization of the Helmholtz operator with absorbing boundary conditions and solve the system via an (sparse) LU decomposition (cf., [22]). We note that this setup is quite efficient already since the LU decomposition can be reused for each source. Reduction of the number of sources becomes of paramount importance in 3D where one is forced to use iterative methods whose costs grow linearly with the number of sources. The search direction gK is the gradient of ϕK evaluated at the initial model m0, depicted in Figure 1(b). The gradient is computed in the usual way via the adjoint-state method (cf., [23]). The full gradient as well as the gradients for K=1,5,10 are depicted in Figure 2. The error between the full and approximated gradient, caused by the crosstalk, is depicted in Figure 3. As expected, the error decays as 1/K. The misfit as a function of α for various K, as well as the full misfit (no randomization), is depicted in Figure 4. This shows that the minimizer of ϕK is reasonably close to the minimizer of the full misfit ϕ, even for a relatively small batch size K.True (a) and initial (b) squared-slowness models (s2/km2) and the true reflectivity.
(a)(b)(c)The full gradient is depicted in (a). The approximate gradients for variousK are depicted in (b) K=1, (c) K=5, and (d) K=10. For a relatively small batch size, the approximate gradients already show the main features.
(a)(b)(c)(d)Figure 3
Error in the gradient as a function of the batch sizeK. As expected, the error goes down as 1/K (dashed line).Behavior of misfit for variousK. Shown are five different stochastic realizations and the true misfit (dashed line) for (a) K=1, (b) K=5, and (c) K=10. The stochastic misfits approximate the true misfit fairly well for relatively small batch sizes.
(a)(b)(c)


### 2.2. Stochastic Approximation

A second alternative is to apply specialized stochastic optimization methods to problem (4) directly. This is often referred to as the stochastic approximation (SA). The main idea of such algorithms is to pick a new random realization in each iteration and possibly average over past iterations to suppress the resulting stochasticity. In the context of the full waveform inversion problem, this gives an iterative algorithm of the form
(7)mν+1=mν-γν∇ϕK,ν(mν),
where batch size K can be as small as 1, {γν} represent step sizes taken by the algorithm, and the notation ϕK,ν emphasizes that a new randomization is used at every iteration ν (in contrast with the SAA approach).We discuss theoretical performance results and describe SAA and SA in more detail in the next section.


### 2.3. Accuracy and Efficiency of Randomized Trace Estimation

Efficient calculation of the trace of a positive semidefinite matrix lies at the heart of our approach. Factors that determine the performance of this estimation include the random process for thei.i.d.w’s, the size of the source ensemble K, and the properties of the matrix. Hutchinson’s approximation [13], which is based on w’s drawn from a Rademacher distribution (i.e., random ±1), attains the smallest variance for the estimate of the trace. The variance can be used to bound the error via confidence intervals. However, the variance is not the only measure of the error. In particular, Avron and Toledo [14] derive bounds on the batch size in terms of ϵ and δ, defined as follows. A randomized-trace estimator TK=K-1∑wiTBwi is an (ϵ,δ)-approximation of T=trace(B) if
(8)Pr(|TK-T||T|≤ϵ)≥1-δ.
The expressions for the minimum batch size K for which the relative error is smaller than ϵ with probability δ are listed in Table 1 (adapted from Avron and Toledo [14]). Smaller ϵ’s and δ’s lead to larger K, which in turn leads to more accurate trace estimates with increased probability.Table 1
Summary of bounds, adapted from Avron and Toledo [14].
EstimatorDistribution ofwVariance of one sampleBound onK for (ϵ,δ) boundHutchinsonHK=(1/K)∑j=1Kwj⊤AwjPr(wj=±1)=1/22(∥A∥F2-∑i=1NAii2)6ϵ-2ln(2rank(A)/δ)GaussianGK=(1/K)∑j=1Kwj⊤Awjwj∈N(0,1)2∥A∥F220ϵ-2ln(2/δ)Phase encodedLK=(N/K)∑j=1KwjTℱAℱTwjwj drawn uniformly from {e1,…,eN}n/a2ϵ-2ln(4n2/δ)ln(4/δ)Of course, these bounds depend on the choice of the probability distribution of thei.i.d.w’s and the matrix B. Aside from obtaining the lowest value for K, simplicity of computational implementation is also a consideration. In Table 1, we summarize the performance of four different choices for the w’s, namely,(1) the Rademacher distribution, that is,Pr(w[i]=±1)=1/2, (w[i] denotes the ith element in the vector w) yielding E{w[i]}=0  and E{w[i]2}=1 for i=1⋯N. Aside from the fact that this estimator HK (see Table 1) leads to minimum variance, the advantage of this choice is that it leads to a fast implementation with a small memory imprint. The disadvantage of this method is that the lower bound depends on the rank of A and requires larger K compared to w’s defined by the Gaussian (see Table 1);(2) the standard normal distribution, that is,w[i]∈N(0,1) for i=1⋯N. While the variance for this estimator GK (see Table 1) is larger than the variance for Hk, the lower bound for K does not depend on the size or the rank of A and is the smallest of all four methods. This suggests that we can use a fixed value of K for arbitrarily large matrices. However, this method is known to converge slower than Hutchinson’s for matrices A that have significant energy in the off-diagonals. This choice also requires a more complex implementation with a larger memory imprint;(3) the fast phase-encoded method wherew’s selected uniformly from the canonical basis, that is, from {e1,…,eN}. this estimator
(9)LK=NK∑j=1KwjTFAFTwj,
where ℱ is a unitary (i.e, ℱT=ℱ-1) random mixing matrix. The idea is to mix the matrix B such that its diagonal entries are evenly distributed. This is important since the unit vectors only sample the diagonal of the matrix. The flatter the distribution of the diagonal elements, the faster the convergence (if all the diagonal elements were to be the same, we need only one sample to compute the trace exactly).The lower bounds summarized in Table1 tell us that Gaussian w’s theoretically require the smallest K and hence the fewest PDE solves. However, this result comes at the expense of more complex arithmetic, which can be a practical consideration [8]. Aside from the lowest bound, the estimator based on Gaussian w’s has the additional advantage that the bound on K does not depend on the size or rank of the matrix B. Hutchinson’s method, on the other hand, depends logarithmically on the rank of B but has the reported advantage that it performs well for near diagonal matrices [14]. This has important implications for our application because our matrix B is typically full rank and can be considered nearly diagonal only when our optimization procedure is close to convergence. At the beginning of the optimization, we can expect the residual to be large and a B that is not necessarily diagonal dominant.We conduct the following stylized experiment to illustrate the quality of the different trace estimators. We solve the discretized Helmholtz equation at 5 Hz for a realistic acoustic model with 301 colocated sources and receivers located at10m depth. We compute matrix B=ATA for a residue A given by the difference between simulation results for the hard and smooth models shown in Figure 1. As expected, the resulting matrix B, shown in Figure 5, contains significant off-diagonal energy. For the phase-encoded part of the experiment, we use a random mixing matrix based on the DFT, as suggested by Romberg [24]. Such mixing matrices are also commonly found in compressive sensing applications [7, 24–26].Figure 5
Residual matrixA=STS, where S is the data residual corresponding to the smooth model depicted in Figure 1(a) at 5 Hz.We evaluated the different trace estimators 1000 times for batch sizes ranging fromK=1⋯1000. The probability for the error level ϵ is estimated by counting the number of times we were able to achieve that error level for each K. The results for the different trace estimators and error levels are summarized in Figure 6. For this particular example, we see little difference in performance between the different estimators. The corresponding theoretical bounds on the batch size, as given by Table 1, are listed in Table 2. Clearly, these bounds are overly pessimistic in this case. In our experiments, we observed that we get similar reconstruction behavior if we use a finer source/receiver sampling. This suggests that the gain in efficiency will increase with the data size, since we can use larger batch sizes for a fixed downsampling ratio. We also noticed, in this particular example, little or no change in behavior if we change the frequency.Table 2
This table shows the theoretical lower bounds (see Table1) on the batch size K for δ=10-1 for the matrix shown in Figure 5.
ϵ=10-1ϵ=10-2ϵ=10-3Gauss6·1036·1056·107Hutchinson4·1034·1054·107Phase9·1039·1059·107Reconstruction as a function ofK for various methods and error levels: (a) ϵ=10-1, (b) ϵ=10-2, and (c) ϵ=10-3.
(a)(b)(c)


## 2.1. Sample Average Approximation

A natural approach to take is to replace the expectation overw by an ensemble average
(5)ϕK(m)=1K∑k=1Kϕ(m,Qwk,Dwk).
This is often referred to in the literature as the sample average approximation (SAA). The random crosstalk can be controlled by picking a “large enough” batch size. As long as the required batch size is smaller than the actual number of sources, we reduce the computational complexity.For a fixedm, it is known that the error |ϕ-ϕK| is of order 1/K (cf., [14]). However, it is not the value of the misfit that we are trying to approximate, but the minimizer. Unfortunately, the difference between the minimizers of ϕ and ϕK is not readily analyzed. Instead, we perform a small numerical experiment to get some idea of the performance of the SAA approach for waveform inversion.We investigate the misfit along the direction of the negative gradientgk (defined below)
(6)fK(α)=ϕK(m-αgK).
The data are generated for the model depicted in Figure 1(a), for 61 colocated, equidistributed sources and receivers along a straight line at 10 m depth and 7 randomly chosen frequencies between 5 and 30 Hz. The source signature is a Ricker wavelet with a peak frequency of 10 Hz. We use a 9-point discretization of the Helmholtz operator with absorbing boundary conditions and solve the system via an (sparse) LU decomposition (cf., [22]). We note that this setup is quite efficient already since the LU decomposition can be reused for each source. Reduction of the number of sources becomes of paramount importance in 3D where one is forced to use iterative methods whose costs grow linearly with the number of sources. The search direction gK is the gradient of ϕK evaluated at the initial model m0, depicted in Figure 1(b). The gradient is computed in the usual way via the adjoint-state method (cf., [23]). The full gradient as well as the gradients for K=1,5,10 are depicted in Figure 2. The error between the full and approximated gradient, caused by the crosstalk, is depicted in Figure 3. As expected, the error decays as 1/K. The misfit as a function of α for various K, as well as the full misfit (no randomization), is depicted in Figure 4. This shows that the minimizer of ϕK is reasonably close to the minimizer of the full misfit ϕ, even for a relatively small batch size K.True (a) and initial (b) squared-slowness models (s2/km2) and the true reflectivity.
(a)(b)(c)The full gradient is depicted in (a). The approximate gradients for variousK are depicted in (b) K=1, (c) K=5, and (d) K=10. For a relatively small batch size, the approximate gradients already show the main features.
(a)(b)(c)(d)Figure 3
Error in the gradient as a function of the batch sizeK. As expected, the error goes down as 1/K (dashed line).Behavior of misfit for variousK. Shown are five different stochastic realizations and the true misfit (dashed line) for (a) K=1, (b) K=5, and (c) K=10. The stochastic misfits approximate the true misfit fairly well for relatively small batch sizes.
(a)(b)(c)


## 2.2. Stochastic Approximation

A second alternative is to apply specialized stochastic optimization methods to problem (4) directly. This is often referred to as the stochastic approximation (SA). The main idea of such algorithms is to pick a new random realization in each iteration and possibly average over past iterations to suppress the resulting stochasticity. In the context of the full waveform inversion problem, this gives an iterative algorithm of the form
(7)mν+1=mν-γν∇ϕK,ν(mν),
where batch size K can be as small as 1, {γν} represent step sizes taken by the algorithm, and the notation ϕK,ν emphasizes that a new randomization is used at every iteration ν (in contrast with the SAA approach).We discuss theoretical performance results and describe SAA and SA in more detail in the next section.


## 2.3. Accuracy and Efficiency of Randomized Trace Estimation

Efficient calculation of the trace of a positive semidefinite matrix lies at the heart of our approach. Factors that determine the performance of this estimation include the random process for thei.i.d.w’s, the size of the source ensemble K, and the properties of the matrix. Hutchinson’s approximation [13], which is based on w’s drawn from a Rademacher distribution (i.e., random ±1), attains the smallest variance for the estimate of the trace. The variance can be used to bound the error via confidence intervals. However, the variance is not the only measure of the error. In particular, Avron and Toledo [14] derive bounds on the batch size in terms of ϵ and δ, defined as follows. A randomized-trace estimator TK=K-1∑wiTBwi is an (ϵ,δ)-approximation of T=trace(B) if
(8)Pr(|TK-T||T|≤ϵ)≥1-δ.
The expressions for the minimum batch size K for which the relative error is smaller than ϵ with probability δ are listed in Table 1 (adapted from Avron and Toledo [14]). Smaller ϵ’s and δ’s lead to larger K, which in turn leads to more accurate trace estimates with increased probability.Table 1
Summary of bounds, adapted from Avron and Toledo [14].
EstimatorDistribution ofwVariance of one sampleBound onK for (ϵ,δ) boundHutchinsonHK=(1/K)∑j=1Kwj⊤AwjPr(wj=±1)=1/22(∥A∥F2-∑i=1NAii2)6ϵ-2ln(2rank(A)/δ)GaussianGK=(1/K)∑j=1Kwj⊤Awjwj∈N(0,1)2∥A∥F220ϵ-2ln(2/δ)Phase encodedLK=(N/K)∑j=1KwjTℱAℱTwjwj drawn uniformly from {e1,…,eN}n/a2ϵ-2ln(4n2/δ)ln(4/δ)Of course, these bounds depend on the choice of the probability distribution of thei.i.d.w’s and the matrix B. Aside from obtaining the lowest value for K, simplicity of computational implementation is also a consideration. In Table 1, we summarize the performance of four different choices for the w’s, namely,(1) the Rademacher distribution, that is,Pr(w[i]=±1)=1/2, (w[i] denotes the ith element in the vector w) yielding E{w[i]}=0  and E{w[i]2}=1 for i=1⋯N. Aside from the fact that this estimator HK (see Table 1) leads to minimum variance, the advantage of this choice is that it leads to a fast implementation with a small memory imprint. The disadvantage of this method is that the lower bound depends on the rank of A and requires larger K compared to w’s defined by the Gaussian (see Table 1);(2) the standard normal distribution, that is,w[i]∈N(0,1) for i=1⋯N. While the variance for this estimator GK (see Table 1) is larger than the variance for Hk, the lower bound for K does not depend on the size or the rank of A and is the smallest of all four methods. This suggests that we can use a fixed value of K for arbitrarily large matrices. However, this method is known to converge slower than Hutchinson’s for matrices A that have significant energy in the off-diagonals. This choice also requires a more complex implementation with a larger memory imprint;(3) the fast phase-encoded method wherew’s selected uniformly from the canonical basis, that is, from {e1,…,eN}. this estimator
(9)LK=NK∑j=1KwjTFAFTwj,
where ℱ is a unitary (i.e, ℱT=ℱ-1) random mixing matrix. The idea is to mix the matrix B such that its diagonal entries are evenly distributed. This is important since the unit vectors only sample the diagonal of the matrix. The flatter the distribution of the diagonal elements, the faster the convergence (if all the diagonal elements were to be the same, we need only one sample to compute the trace exactly).The lower bounds summarized in Table1 tell us that Gaussian w’s theoretically require the smallest K and hence the fewest PDE solves. However, this result comes at the expense of more complex arithmetic, which can be a practical consideration [8]. Aside from the lowest bound, the estimator based on Gaussian w’s has the additional advantage that the bound on K does not depend on the size or rank of the matrix B. Hutchinson’s method, on the other hand, depends logarithmically on the rank of B but has the reported advantage that it performs well for near diagonal matrices [14]. This has important implications for our application because our matrix B is typically full rank and can be considered nearly diagonal only when our optimization procedure is close to convergence. At the beginning of the optimization, we can expect the residual to be large and a B that is not necessarily diagonal dominant.We conduct the following stylized experiment to illustrate the quality of the different trace estimators. We solve the discretized Helmholtz equation at 5 Hz for a realistic acoustic model with 301 colocated sources and receivers located at10m depth. We compute matrix B=ATA for a residue A given by the difference between simulation results for the hard and smooth models shown in Figure 1. As expected, the resulting matrix B, shown in Figure 5, contains significant off-diagonal energy. For the phase-encoded part of the experiment, we use a random mixing matrix based on the DFT, as suggested by Romberg [24]. Such mixing matrices are also commonly found in compressive sensing applications [7, 24–26].Figure 5
Residual matrixA=STS, where S is the data residual corresponding to the smooth model depicted in Figure 1(a) at 5 Hz.We evaluated the different trace estimators 1000 times for batch sizes ranging fromK=1⋯1000. The probability for the error level ϵ is estimated by counting the number of times we were able to achieve that error level for each K. The results for the different trace estimators and error levels are summarized in Figure 6. For this particular example, we see little difference in performance between the different estimators. The corresponding theoretical bounds on the batch size, as given by Table 1, are listed in Table 2. Clearly, these bounds are overly pessimistic in this case. In our experiments, we observed that we get similar reconstruction behavior if we use a finer source/receiver sampling. This suggests that the gain in efficiency will increase with the data size, since we can use larger batch sizes for a fixed downsampling ratio. We also noticed, in this particular example, little or no change in behavior if we change the frequency.Table 2
This table shows the theoretical lower bounds (see Table1) on the batch size K for δ=10-1 for the matrix shown in Figure 5.
ϵ=10-1ϵ=10-2ϵ=10-3Gauss6·1036·1056·107Hutchinson4·1034·1054·107Phase9·1039·1059·107Reconstruction as a function ofK for various methods and error levels: (a) ϵ=10-1, (b) ϵ=10-2, and (c) ϵ=10-3.
(a)(b)(c)


## 3. Optimization

### 3.1. Sample Average Approximation

The sample average approximation (SAA) is used to solve the following class of stochastic optimization problems:(10)min⁡x∈X⁡{f(x)=Ew{F(x,w)}},
where X⊂ℝn is the set of admissible models (assumed to be a compact convex set, for example, box constraints xmin≤x≤xmax), w is a random vector with distribution supported on W⊂ℝd, F:X×W→ℝ, and the function f(x) is convex [27]. The last important assumption is the law of large numbers (LLN), that is, f̂K(x)→f(x) with probability 1 as K→∞. These assumptions are required for most of the known theoretical results about convergence of SAA methods. The convexity assumption and LLN assumption can be relaxed in the case when F(·,w) is continuous on X for almost every w∈W and F(x,w) is dominated by an integrable function G(w), so that |f(x)|≤Ew{G(w)} for every x∈X [28]. Given an optimization problem of type (10), the SAA approach [27] is to generate a random sample w1,…,wK and solve the approximate (or sample average) problem
(11)min⁡x∈X⁡{f̂K(x)=1K∑j=1KF(x,wj)}.
When these assumptions are satisfied, the optimal value of (11) converges to the optimal value of the full problem (10) with probability 1. Moreover, under more technical assumptions on the distribution of the random variable w, conservative bounds have been derived on the batch size K necessary to obtain a particular accuracy level ϵ [29, equation (22)]. These bounds do not require the convexity assumptions but instead require assumptions on local behavior of F(·,w). It is worth underscoring that “accuracy” here of solution x¯ with respect to the optimal solution x* is defined with respect to the function value difference f(x)-f(x*), rather than in terms of ∥x-x*∥ or other measure in the space of model parameters. From a practical point of view, the SAA approach is appealing because it allows flexibility in the choice of algorithm for the solution of (11). This works on two levels. First, if a faster algorithm becomes available for the solution of (11), it can immediately impact (10). Second, having fixed a large K¯ and f̂K¯ to obtain reasonable accuracy in the solution of (10), one is free to approximately solve a sequence of smaller problems (Ki≪K¯) with warm starts on the way to solving f̂K¯ [12]. In other words, SAA theory guarantees the existence of a K large enough for which the approximate problem is close to the full problem; however, the algorithm for solving the approximate problem (11) is left completely to the practitioner and in particular may require the evaluation of very few samples at early iterations.


### 3.2. Stochastic Approximation

Stochastic approximation (SA) methods go back to Robbins and Monro [30], who considered the root-finding problem
(12)g(x)=g0,
in the case where g(x) cannot be evaluated directly. Rather, one has access to a function G(x,w) for which Ew{G(x,w)}=g(x). The approach can be translated to optimization problems of the form
(13)min⁡⁡f(x)
by considering g to be the gradient of f and setting g0=0. Again, we cannot evaluate f(x) directly, but we have access to F(x,w) for which Ew{F(x,w)}=f(x). More generally, for problems of type (10), Bertsekas and Tsitsiklis [31] and Bertsekas and Tsitsiklis [32] consider iterative algorithms of the form
(14)xν+1=xν-γνsν,
where γν are a sequence of step sizes determined a priori that satisfy certain properties, and sν can be thought of as noisy unbiased estimates of the gradient (i.e., Ewsν=∇f(xν)). Note that right away we are forced into an algorithmic framework, which never appears in the SAA discussion. The positive step sizes γν are chosen to satisfy
(15)∑ν=0∞γν=∞,∑ν=0∞γν2<∞.
The main idea is that the step sizes go to zero, but not too fast. A commonly used example of such a sequence of step sizes is
(16)γν∝1ν.
The main result of Bertsekas and Tsitsiklis [31] is that if ∇f satisfies the Lipshitz condition with constant L(17)‖∇f(x)-∇f(y)‖≤L‖x-y‖,
that is, the changes in the gradient are bounded in norm by changes in the parameter space, and if the directions sν on average point “close to” the gradient and are not too noisy, then the sequence f(xν) converges, and every limit point x¯ of {xν} is a stationary point of f (i.e., ∇f(x¯)=0). Under stronger assumptions that the level sets of f are bounded and the minimum is unique, this guarantees that the algorithms described above will find it. A similar family of algorithms was studied by Polyak and Juditsky [33], who considered larger step sizes γν but included averaging model estimates into their algorithm. In the context discussed above, the step size rule 10 is replaced by
(18)γν-γν+1γν=o(γν).
A particular example of such a sequence cited by the paper is
(19)γν∝ν-β,0<β<1.
The iterative scheme is then given by
(20)xν+1=xν-γνsν,x¯ν=1ν∑i=0ν-1xi.
Under assumptions similar in spirit to the ones in Bertsekas and Tsitsiklis [31], there is a result for the convergence of the iterates xν to the true estimate x*, namely, x¯ν→x* almost surely and
(21)ν(x¯ν-x*)⟶DN(0,V),
where the convergence is in distribution, and the matrix V is in some sense optimal and is related to the Hessian of f at the solution x*. A more recent report [34] also considers averaging of model iterates in the context of optimizing (not necessarily smooth) convex functions of the form
(22)f(x)=Ew{F(x,w)}
over a convex set X. When f is smooth, this situation reduces to the previous discussion. Nesterov and Vial [34] choose a finite sequence of N step sizes a prior and consider the error in the expected value function
(23)Ex¯ν{f(x¯ν)}-f(x*)
after N iterations. This is similar to the SAA analysis but is much easier to interpret, because now the desired accuracy in the objective value directly translates to the number of iterations of a particular algorithm
(24)xν+1=πX(xν-γνsν),x¯=∑ν=0N-1γνxν∑ν=0N-1γν,
where πX is projection onto the convex set of admissible models X. Unfortunately, the error is O(L2(∑γν2/∑γν)+R2(1/∑γν)), where R is the diameter of the set X (related to the bounds on x from our earlier example) and L is a uniform bound on ∥∇f∥, and so the estimate may be overly conservative. If all the γν are chosen to be uniform, the optimal size is γ=R/LN, and then the result is simply
(25)Ex¯ν{f(x¯ν)}-f(x*)≤LRN.
For a recent survey of stochastic optimization and new robust SA methods, please see Nemirovski et al. [27].Note that the error rate in the objective values isO(1/N), where the constant depends in a straightforward way on the size of the set X and the behavior of ∥∇f∥. Compare this to the O(1/K)error bound for the SAA approach. In contrast to the SAA, the SA approach translates directly into a particular algorithm. This makes it easier to implement for full waveform inversion, but also leaves less freedom for algorithm design than in SAA, where any algorithm can be used to solve the deterministic ensemble average problem.


## 3.1. Sample Average Approximation

The sample average approximation (SAA) is used to solve the following class of stochastic optimization problems:(10)min⁡x∈X⁡{f(x)=Ew{F(x,w)}},
where X⊂ℝn is the set of admissible models (assumed to be a compact convex set, for example, box constraints xmin≤x≤xmax), w is a random vector with distribution supported on W⊂ℝd, F:X×W→ℝ, and the function f(x) is convex [27]. The last important assumption is the law of large numbers (LLN), that is, f̂K(x)→f(x) with probability 1 as K→∞. These assumptions are required for most of the known theoretical results about convergence of SAA methods. The convexity assumption and LLN assumption can be relaxed in the case when F(·,w) is continuous on X for almost every w∈W and F(x,w) is dominated by an integrable function G(w), so that |f(x)|≤Ew{G(w)} for every x∈X [28]. Given an optimization problem of type (10), the SAA approach [27] is to generate a random sample w1,…,wK and solve the approximate (or sample average) problem
(11)min⁡x∈X⁡{f̂K(x)=1K∑j=1KF(x,wj)}.
When these assumptions are satisfied, the optimal value of (11) converges to the optimal value of the full problem (10) with probability 1. Moreover, under more technical assumptions on the distribution of the random variable w, conservative bounds have been derived on the batch size K necessary to obtain a particular accuracy level ϵ [29, equation (22)]. These bounds do not require the convexity assumptions but instead require assumptions on local behavior of F(·,w). It is worth underscoring that “accuracy” here of solution x¯ with respect to the optimal solution x* is defined with respect to the function value difference f(x)-f(x*), rather than in terms of ∥x-x*∥ or other measure in the space of model parameters. From a practical point of view, the SAA approach is appealing because it allows flexibility in the choice of algorithm for the solution of (11). This works on two levels. First, if a faster algorithm becomes available for the solution of (11), it can immediately impact (10). Second, having fixed a large K¯ and f̂K¯ to obtain reasonable accuracy in the solution of (10), one is free to approximately solve a sequence of smaller problems (Ki≪K¯) with warm starts on the way to solving f̂K¯ [12]. In other words, SAA theory guarantees the existence of a K large enough for which the approximate problem is close to the full problem; however, the algorithm for solving the approximate problem (11) is left completely to the practitioner and in particular may require the evaluation of very few samples at early iterations.


## 3.2. Stochastic Approximation

Stochastic approximation (SA) methods go back to Robbins and Monro [30], who considered the root-finding problem
(12)g(x)=g0,
in the case where g(x) cannot be evaluated directly. Rather, one has access to a function G(x,w) for which Ew{G(x,w)}=g(x). The approach can be translated to optimization problems of the form
(13)min⁡⁡f(x)
by considering g to be the gradient of f and setting g0=0. Again, we cannot evaluate f(x) directly, but we have access to F(x,w) for which Ew{F(x,w)}=f(x). More generally, for problems of type (10), Bertsekas and Tsitsiklis [31] and Bertsekas and Tsitsiklis [32] consider iterative algorithms of the form
(14)xν+1=xν-γνsν,
where γν are a sequence of step sizes determined a priori that satisfy certain properties, and sν can be thought of as noisy unbiased estimates of the gradient (i.e., Ewsν=∇f(xν)). Note that right away we are forced into an algorithmic framework, which never appears in the SAA discussion. The positive step sizes γν are chosen to satisfy
(15)∑ν=0∞γν=∞,∑ν=0∞γν2<∞.
The main idea is that the step sizes go to zero, but not too fast. A commonly used example of such a sequence of step sizes is
(16)γν∝1ν.
The main result of Bertsekas and Tsitsiklis [31] is that if ∇f satisfies the Lipshitz condition with constant L(17)‖∇f(x)-∇f(y)‖≤L‖x-y‖,
that is, the changes in the gradient are bounded in norm by changes in the parameter space, and if the directions sν on average point “close to” the gradient and are not too noisy, then the sequence f(xν) converges, and every limit point x¯ of {xν} is a stationary point of f (i.e., ∇f(x¯)=0). Under stronger assumptions that the level sets of f are bounded and the minimum is unique, this guarantees that the algorithms described above will find it. A similar family of algorithms was studied by Polyak and Juditsky [33], who considered larger step sizes γν but included averaging model estimates into their algorithm. In the context discussed above, the step size rule 10 is replaced by
(18)γν-γν+1γν=o(γν).
A particular example of such a sequence cited by the paper is
(19)γν∝ν-β,0<β<1.
The iterative scheme is then given by
(20)xν+1=xν-γνsν,x¯ν=1ν∑i=0ν-1xi.
Under assumptions similar in spirit to the ones in Bertsekas and Tsitsiklis [31], there is a result for the convergence of the iterates xν to the true estimate x*, namely, x¯ν→x* almost surely and
(21)ν(x¯ν-x*)⟶DN(0,V),
where the convergence is in distribution, and the matrix V is in some sense optimal and is related to the Hessian of f at the solution x*. A more recent report [34] also considers averaging of model iterates in the context of optimizing (not necessarily smooth) convex functions of the form
(22)f(x)=Ew{F(x,w)}
over a convex set X. When f is smooth, this situation reduces to the previous discussion. Nesterov and Vial [34] choose a finite sequence of N step sizes a prior and consider the error in the expected value function
(23)Ex¯ν{f(x¯ν)}-f(x*)
after N iterations. This is similar to the SAA analysis but is much easier to interpret, because now the desired accuracy in the objective value directly translates to the number of iterations of a particular algorithm
(24)xν+1=πX(xν-γνsν),x¯=∑ν=0N-1γνxν∑ν=0N-1γν,
where πX is projection onto the convex set of admissible models X. Unfortunately, the error is O(L2(∑γν2/∑γν)+R2(1/∑γν)), where R is the diameter of the set X (related to the bounds on x from our earlier example) and L is a uniform bound on ∥∇f∥, and so the estimate may be overly conservative. If all the γν are chosen to be uniform, the optimal size is γ=R/LN, and then the result is simply
(25)Ex¯ν{f(x¯ν)}-f(x*)≤LRN.
For a recent survey of stochastic optimization and new robust SA methods, please see Nemirovski et al. [27].Note that the error rate in the objective values isO(1/N), where the constant depends in a straightforward way on the size of the set X and the behavior of ∥∇f∥. Compare this to the O(1/K)error bound for the SAA approach. In contrast to the SAA, the SA approach translates directly into a particular algorithm. This makes it easier to implement for full waveform inversion, but also leaves less freedom for algorithm design than in SAA, where any algorithm can be used to solve the deterministic ensemble average problem.


## 4. Algorithms

To test the performance of the SAA approach, we chose to use a steepest descent method with an Armijo line search (cf., [35]). Although one could in principle use a second-order method (such as L-BFGS), we chose to use a first-order method to allow for better comparison to the SA results. The pseudocode is presented in Algorithm 1.Algorithm 1:Steepest descent.
While not converged dos←-∇ϕ[mi]/∥∇ϕ[mi]∥2findλ s.t. ϕ[mi+λs]≤ϕ[mi]+cλ∇ϕ[mi]Tsmi+1←mi+λsi←i+1end whileThe SA methods are closely related to the steepest descent method. The main difference is that for each iteration a new random realization is drawn from a prescribed distribution and that the result is averaged over past iterations. We chose to implement a few modifications to the standard SA algorithms. First, we use an Armijo line search to determine the step size instead of using a prescribed sequence such as that discussed in the previous section. This assures some descent at each iteration with respect to the current realization ofϕK, and we found that this greatly improved the convergence. Second, we allow for averaging over the past n iterations instead of the full history. This prevents the method from stalling. The pseudo-code is presented in Algorithm 2.Algorithm 2:Stochastic descent.
While not converged dodraww from a pre-scribed distributions←-∇ϕ[mi,w]/∥∇ϕ[mi,w]∥2findλ s.t. ϕ[mi+λs,w]≤ϕ[mi,w]+cλ∇ϕ[mi,w]Tsmi+1←(1/(n+1))(∑i-nimi+λs)i←i+1end while


## 5. Results

For the numerical experiments, we use the true and initial squared-slowness models depicted in Figure1. The data are generated for 61 equispaced, colocated sources and receivers at 10 m depth and 7 randomly chosen (but fixed) frequencies between 5 and 30 Hz. The latter strategy is inspired by results from compressive sensing (cf., [7, 36, 37]). The basic idea is to turn aliases that are introduced by sub-Nyquist sampling into random noise.The Helmholtz operator is discretized on a grid with 10 m spacing, using a 9-point finite difference stencil and absorbing boundary conditions. The point sources are represented as narrow Gaussians. As a source signature, we use a Ricker wavelet with a peak frequency of 10 Hz. The noise is Gaussian with a prescribed SNR.We run each of the optimization methods for 500 iterations and compare the performance for various batch sizes and noise levels to the result of steepest descent on the full problem. Remember that by using small batch sizes, the iterations are very cheap, so we can afford to do more. The random vectors are drawn from a Gaussian distribution with zero mean and unit variance. We chose to use the Gaussian because the theoretical bounds onK do not depend on properties of the residual matrix. Although the matrix will change constantly during the optimization, we can at least expect a uniform quality of the approximation.In a realistic application, one might want to add a regularization term. In particular, this would prevent the overfitting that we observe in the noisy case. Note that limiting the amount of iterations also serves as a form of regularization [38].


### 5.1. Sample Average Approximation

We choose a set ofK Gaussian random vectors with zero mean and unit variance and run the steepest descent algorithm presented previously on the resulting deterministic optimization problem. The results after 500 iterations on data without noise are shown in the first column of Figure 7. The error between the recovered and true model is shown in Figure 8(a). As reference, the error between the true and recovered model for the inversion with all the sequential sources is also shown. As expected, the recovery is better for larger batch sizes. The recovered models for data with noise are shown in the second column (SNR = 20 dB) and third (SNR = 10 dB) columns of Figure 7. The corresponding recovery error is shown in Figures 8(b) and 8(c), respectively. It shows that the SAA approach starts overfitting in an earlier stage than the full inversion. Also, we are not able to reach the same model error as the full inversion.Inversion result for the SAA approach with various batch sizes and noise levels. The rows represent different batch sizesK=1,5,10,20, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. The reconstruction with K=20 for noiseless data (j) is qualitatively comparable to the full reconstruction. The quality deteriorates quickly for small batch sizes and noisy data.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Error between the inverted and true model for the SAA approach with various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). On noiseless data, we achieve a qualitatively comparable result withK=20, as can be seen from (a). For noisy data, however, the largest batch size is not enough to prevent overfitting.
(a)(b)(c)


### 5.2. Stochastic Approximation

We run the stochastic descent algorithm for varying batch sizes (K=1,5,10) and history sizes (n=0,10,500).The results obtained without averaging are shown in Figure9. The columns represent different batch sizes, while the rows represent different noise levels. The recovery errors for the different batch sizes and noise levels are shown in Figure 10. In the noiseless case, we are able to achieve the same recovery error as the full inversion with only one simultaneous source. When noise is present in the data, one simultaneous source is not enough, however. Still, we can achieve the same recovery error as the full problem with only 10 simultaneous sources. This yields an order of magnitude improvement in our computation, since the total number of iterations needed by the stochastic method to achieve a given level of accuracy is roughly the same as required by a deterministic first-order method used on the full system, but each stochastic iteration requires ten times fewer PDE solves than a deterministic iteration on the full system.Inversion result for the SA approach without averaging for various batch sizes and noise levels. The rows represent different batch sizesK=1,5,10, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. We obtain good results with K=1, and the quality does not improve dramatically for larger batch sizes, except for the highest noise level.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Error between the inverted and true model for the SA approach without averaging for various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). We get qualitatively similar results, compared to the full inversion, withK=1 for noiseless data and data with 10 dB of noise. For very noisy data (20 dB), we need a larger batch size. Although the SA approach requires roughly the same number of iterations as the full inversion, the iterations are much cheaper. For K=1, we model the data for only one simultaneous source per iteration, compared to 61 for the full inversion.
(a)(b)(c)Results obtained with averaging over the past 10 iterations are shown in Figure11. The rows represent different batch sizes, while the columns represent different noise levels. The corresponding recovery errors are shown in Figure 12. It shows that averaging helps to overcome some of the noise sensitivity, and we are now able to achieve a good reconstruction with only 5 simultaneous sources. Also, the averaging damps the irregularity of the convergence somewhat.Inversion result for the SA approach with limited averaging (n=10) for various batch sizes and noise levels. The rows represent different batch sizes K=1,5,10, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. We obtain good results with K=1, and the quality does not improve dramatically for larger batch sizes, except for the highest noise level.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Error between the inverted and true model for the SA approach with limited averaging for various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). The convergence is smoother than that of SA without averaging, especially when the data is very noisy (10 dB). The averaging seems to slow down the convergence slightly, however, and we need a batch sizeK=5 for the best results.
(a)(b)(c)Finally, we show the result obtained by averaging over the full history in Figure13. The corresponding recovery error is shown in Figure 14. It shows that too much averaging slows down the convergence.Inversion result for the SA approach with full averaging (n=500) for various batch sizes and noise levels. The rows represent different batch sizes K=1,5,10, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. Averaging over the full past dramatically deteriorates the reconstruction.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Error between the inverted and true model for the SA approach with full averaging for various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). Averaging over the full past slows down the convergence dramatically.
(a)(b)(c)


## 5.1. Sample Average Approximation

We choose a set ofK Gaussian random vectors with zero mean and unit variance and run the steepest descent algorithm presented previously on the resulting deterministic optimization problem. The results after 500 iterations on data without noise are shown in the first column of Figure 7. The error between the recovered and true model is shown in Figure 8(a). As reference, the error between the true and recovered model for the inversion with all the sequential sources is also shown. As expected, the recovery is better for larger batch sizes. The recovered models for data with noise are shown in the second column (SNR = 20 dB) and third (SNR = 10 dB) columns of Figure 7. The corresponding recovery error is shown in Figures 8(b) and 8(c), respectively. It shows that the SAA approach starts overfitting in an earlier stage than the full inversion. Also, we are not able to reach the same model error as the full inversion.Inversion result for the SAA approach with various batch sizes and noise levels. The rows represent different batch sizesK=1,5,10,20, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. The reconstruction with K=20 for noiseless data (j) is qualitatively comparable to the full reconstruction. The quality deteriorates quickly for small batch sizes and noisy data.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)Error between the inverted and true model for the SAA approach with various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). On noiseless data, we achieve a qualitatively comparable result withK=20, as can be seen from (a). For noisy data, however, the largest batch size is not enough to prevent overfitting.
(a)(b)(c)


## 5.2. Stochastic Approximation

We run the stochastic descent algorithm for varying batch sizes (K=1,5,10) and history sizes (n=0,10,500).The results obtained without averaging are shown in Figure9. The columns represent different batch sizes, while the rows represent different noise levels. The recovery errors for the different batch sizes and noise levels are shown in Figure 10. In the noiseless case, we are able to achieve the same recovery error as the full inversion with only one simultaneous source. When noise is present in the data, one simultaneous source is not enough, however. Still, we can achieve the same recovery error as the full problem with only 10 simultaneous sources. This yields an order of magnitude improvement in our computation, since the total number of iterations needed by the stochastic method to achieve a given level of accuracy is roughly the same as required by a deterministic first-order method used on the full system, but each stochastic iteration requires ten times fewer PDE solves than a deterministic iteration on the full system.Inversion result for the SA approach without averaging for various batch sizes and noise levels. The rows represent different batch sizesK=1,5,10, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. We obtain good results with K=1, and the quality does not improve dramatically for larger batch sizes, except for the highest noise level.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Error between the inverted and true model for the SA approach without averaging for various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). We get qualitatively similar results, compared to the full inversion, withK=1 for noiseless data and data with 10 dB of noise. For very noisy data (20 dB), we need a larger batch size. Although the SA approach requires roughly the same number of iterations as the full inversion, the iterations are much cheaper. For K=1, we model the data for only one simultaneous source per iteration, compared to 61 for the full inversion.
(a)(b)(c)Results obtained with averaging over the past 10 iterations are shown in Figure11. The rows represent different batch sizes, while the columns represent different noise levels. The corresponding recovery errors are shown in Figure 12. It shows that averaging helps to overcome some of the noise sensitivity, and we are now able to achieve a good reconstruction with only 5 simultaneous sources. Also, the averaging damps the irregularity of the convergence somewhat.Inversion result for the SA approach with limited averaging (n=10) for various batch sizes and noise levels. The rows represent different batch sizes K=1,5,10, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. We obtain good results with K=1, and the quality does not improve dramatically for larger batch sizes, except for the highest noise level.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Error between the inverted and true model for the SA approach with limited averaging for various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). The convergence is smoother than that of SA without averaging, especially when the data is very noisy (10 dB). The averaging seems to slow down the convergence slightly, however, and we need a batch sizeK=5 for the best results.
(a)(b)(c)Finally, we show the result obtained by averaging over the full history in Figure13. The corresponding recovery error is shown in Figure 14. It shows that too much averaging slows down the convergence.Inversion result for the SA approach with full averaging (n=500) for various batch sizes and noise levels. The rows represent different batch sizes K=1,5,10, while the columns represent different noise levels: no noise, SNR = 20 dB, and SNR = 10 dB. Averaging over the full past dramatically deteriorates the reconstruction.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Error between the inverted and true model for the SA approach with full averaging for various batch sizes and the full problem, (a) without noise, (b) with noise (SNR = 20 dB), and (c) with noise (SNR = 10 dB). Averaging over the full past slows down the convergence dramatically.
(a)(b)(c)


## 6. Conclusions and Discussion

Following Haber et al. [39], we reduce the dimensionality of full waveform inversion via randomized trace estimation. This reduction comes at the cost of introducing random crosstalk between the sources into the updates. The resulting optimization problem can be treated as a stochastic optimization problem. Theory for such methods goes back to the 1950s and justifies the approach presented by Krebs et al. [8]. In particular, we use theoretical results by Avron and Toledo [14] on randomized trace estimation to get bounds for the batch size needed to approximate the misfit to a given accuracy level with a given probability. Numerical tests show, however, that these bounds may be overly pessimistic and that we get reasonable approximations for modest batch sizes.Theory from the field of stochastic optimization suggests several approaches to tackle the optimization problem and reduce the influence of the crosstalk introduced by the randomization. The first approach, thesample average approximation, dictates the use of a fixed set of random sources and relies solely on increasing the batch size to get rid of the crosstalk. The stochastic approximation, on the other hand, dictates that we redraw the randomization each iteration and average over the past in order to suppress the stochasticity of the gradients.We note that, as opposed torandomized dimensionality reduction, several authors have proposed methods for deterministic dimensionality reduction [39, 40]. These techniques are related to optimal experimental design and try to determine the source combination that somehow optimally illuminates the target. It is not quite clear how such methods compare to the randomized approach discussed here. It is clear, however, that by using random superpositions we have access to powerful results from the field of compressive sensing to further improve the reconstruction. Li and Herrmann [11] use sparse recovery techniques instead of Monte Carlo sampling to get rid of the crosstalk.In our experiments, we were able to obtain results that are comparable to the full optimization with a small fraction of the number of sources. In the noiseless case, we needed onlyone simultaneous source for the SA approach. Even with noisy data, five simultaneous sources proved sufficient. This is a very promising result, since using five simultaneous sources for the SA method means that every iteration requires 20 times fewer PDE solves, which directly translates to a 20x computational speedup compared to a first-order deterministic method. The key point is that both SA and the full deterministic approach require roughly the same number of iterations to achieve the same accuracy.Averaging over a limited number of past iterations improved the results for a fixed batch size and allows for the use of fewer simultaneous sources. However, too much averaging slows down the convergence.The results of the SA approach, where a new realization of the random vectors is drawn at every iteration, are superior to the SAA results, where the random vectors are fixed. However, one could use a more sophisticated (possibly black box) optimization method for the SAA approach to get a similar result with fewer iterations. The tradeoff between using a smaller batch size and first-order methods (i.e., more iterations) versus using a larger batch size and second-order methods (i.e., less iterations) needs to be investigated further. Random superposition of shots only makes sense if those shots are sampled by the same receivers. In particular, this hampers straightforward application to marine seismic data. One way to get around this is to partition the data into blocks that are fully sampled. However, this would not give the same amount of reduction in the number of shots because only shots that are relatively close to each other can be combined without losing too much data.The type of encoding used will most likely affect the behavior of both SA and SAA methods. It remains to be investigated which encoding is most suitable for waveform inversion.



---

*Source: 689041-2011-06-26.xml*