# Full waveform inversion using extended and simultaneous sources††thanks:
Corresponding author: Eran Treister. This research was partially supported by
The Israel Science Foundation (grant No. 1589/19).

Sagi Buchatsky22footnotemark: 2 and Eran Treister Department of Computer
Science, Ben-Gurion University of the Negev, Beer Sheva, Israel.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

PDE-constrained optimization problems are often treated using the reduced
formulation where the PDE constraints are eliminated. This approach is known
to be more computationally feasible than other alternatives at large scales.
However, the elimination of the constraints forces the optimization process to
fulfill the constraints at all times. In some problems this may lead to a
highly non-linear objective, which is hard to solve. An example to such a
problem, which we focus on in this work, is Full Waveform Inversion (FWI),
which appears in seismic exploration of oil and gas reservoirs, and medical
imaging. In an attempt to relieve the non-linearity of FWI, several approaches
suggested to expand the optimization search space and relax the PDE
constraints. This comes, however, with severe memory and computational costs,
which we aim to reduce. In this work we adopt the expanded search space
approach, and suggest a new formulation of FWI using extended source
functions. To make the source-extended problem more feasible in memory and
computations, we couple the source extensions in the form of a low-rank
matrix. This way, we have a large-but-manageable additional parameter space,
which has a rather low memory footprint, and is much more suitable for solving
large scale instances of the problem than the full rank additional space. In
addition, we show how our source-extended approach is applied together with
the popular simultaneous sources technique—a stochastic optimization technique
that significantly reduces the computations needed for FWI inversions. We
demonstrate our approaches for solving FWI problems using 2D and 3D models
with high frequency data only.

###### keywords:

Inverse problems, PDE-constrained optimization, Gauss-Newton, Full Waveform
Inversion, Trace estimation, Extended Sources, Low rank minimization,
Alternating minimization.

86A22, 86A15, 65M32, 65N22, 35Q86, 35R30

##  1 Introduction

Many computational science applications involve parameter estimation of
partial differential equations (PDEs), known as inverse problems [15, 47, 38].
One challenging inverse problem is _Full Waveform Inversion_ (FWI), where we
aim to determine the wave velocity, density, and possibly other parameters of
a heterogeneous medium, given sources and waveform observations at receiver
locations on its boundary [25]. Inverse problems of this type arise in seismic
exploration of oil and gas reservoirs, earth sub-surface mapping, ultrasound
imaging [8], optical diffraction tomography [39], brain imaging [14] and more.

The estimation of the velocity model is usually performed by fitting
numerically simulated data to observed field data. This results in an
optimization problem that needs to be solved by an iterative descent
algorithm, which gradually reduces the misfit between the simulated and field
data. This data fitting problem is usually accompanied by a suitable
regularization that aims to introduce prior information on the estimated
coefficients [40]. In this paper we focus on FWI in the context of seismic
exploration, where the wave velocity of the earth subsurface is estimated.
This problem has gained popularity in the last decade with the advances in
data acquisition techniques, computing power and numerical algorithms [27, 12,
21, 37, 10, 45, 44, 25].

Solving the FWI problem is challenging in two main aspects, both of which may
be relevant to other inverse problems. First, it is an ill-posed and highly
non-convex problem, especially if the sources and receivers are placed on the
same surface. Then the problem typically has multiple minima, and convergence
to a local minimum leads to a non-plausible estimated model. Many approaches
have been proposed for solving FWI, however, it is still considered difficult
to solve, as solutions techniques can be unstable, converging to local minima
in many scenarios. Some recent works suggest to complement FWI with other more
stable inverse problems such as travel time tomography [24, 41] or
electromagnetic inversion [19]. Other approaches can be found in [10, 22].
Such techniques are useful in some cases, but still, FWI is considered
challenging and not robust enough in real life.

In addition to being ill-posed, FWI is highly computationally expensive, both
in memory usage and computations. The problem typically involves a large
number of sources and several frequencies, and requires multiple solutions of
the forward problem to simulate the data, each of which is challenging by
itself. These numerical simulations are either obtained by propagating the
wave equation (in time domain), or by solving the Helmholtz equation (in
frequency domain). Both options introduce computational challenges due to the
scale and properties of the problem. The aim of this work is to help dealing
with both of the aforementioned challenges—non-linearity and computational
cost—and make the FWI solution process more robust, but also more
computationally feasible at the same time, involving less forward simulations.

_Relieving the non-linearity by enlarging the search-space_ Our first
objective in this work is to relieve some of non-linearity in FWI. Recent
works [2, 20, 46] have demonstrated that relaxing the PDE constraints and
enlarging the search-space may allow us to avoid local minima, and find
plausible solution of FWI from rather arbitrary initial guesses. This can be
obtained by introducing the full waveforms as variables, or extending the
point-wise sources to the full domain. However, this comes with a severe
memory footprint in 3D, coming from the introduction of the full waveforms or
extended sources as variables of the optimization problem. Both options
introduce several hundreds (or more) vectors of unknowns to keep track of and
iterate on as part of the solution process, each at the size of the 3D domain.

In our first contribution in this paper we suggest a new formulation of FWI
that has the additional search parameters in the form of _low-rank and sparse_
extended source functions. This way, we relax the PDE constraints, while at
the same time we keep the inversion manageable in memory (comparable to the
reduced FWI version), by using low-rank source extensions and sparsity
promoting penalties.

_Reducing computational complexity_ To further reduce the computational
complexity of the extended inversion problem, we adopt a stochastic
optimization technique called _simultaneous sources_ [18]. In this approach we
reduce the number of sources involved in each iteration of the optimization,
by applying a randomized trace estimation technique [7] to the data misfit
term. While standard stochastic techniques are based on random subsets of the
misfit terms (sources), in the simultaneous technique we project all the
sources onto a smaller dimension using a random matrix. This way, the sources
are randomly mixed, and the inversion is not biased towards any geometry of
the sources’ locations as in the subset-based methods [18]. Therefore the
simultaneous sources technique was found to be more effective than the subset-
of-sources technique in general PDE-constrained optimization like DC-
resistivity [32, 30], and FWI in particular [43].

However, as it is, the simultaneous sources technique is not suitable for the
experimental setting of typical real-life FWI scenarios. That includes common
settings where, for example, each source only affects part of the receivers,
or sporadically missing data due to malfunctioning receivers, or cases where
the measurement noise parameters (the co-variance matrices) are not similar
for all the sources. The particular case of sporadically missing data was
addressed in [31] using a smooth data completion regularization for the DC-
resistivity problem. This solution, however, is not suitable for FWI, for
which the data is not smooth. Recently, the work [23] adjusted the
simultaneous sources technique in the case of missing data by splitting the
data misfit function into two terms, so that one is suitable to the
simultaneous sources method and the other is easy to handle computationally
(we elaborate on this later). This was accompanied by a low-rank
regularization to the data term, which is suitable for FWI, in particular.

Our second contribution in this paper is to show the combination of extended
and simultaneous sources, since, similarly to the scenarios mentioned above,
the simultaneous sources is not entirely suitable to accelerate FWI when
extended sources are used. In particular, we adjust the method to the low-rank
source extensions that we suggest in this paper.

The rest of this paper is organized as follows. In Section 2 we introduce some
background and notation used throughout this paper. In Section 3 we introduce
our main contribution: low-rank and sparse extended sources. Then, in Section
4 we show how to merge the extended sources method with simultaneous sources.
Following that, we present numerical results in Section 5.

##  2 Preliminaries and background

To model the waveforms for FWI in frequency domain, we consider the discrete
acoustic Helmholtz equation as a forward problem, assuming a constant density
media

(2.1)
$\displaystyle\mathcal{H}(\mathbf{m},\omega)\mathbf{u}=\Delta_{h}\mathbf{u}+\omega^{2}\mathbf{m}\odot\mathbf{u}=\mathbf{q}_{s}.$

The symbol $\Delta_{h}$ represents the discretized Laplacian on a nodal
regular grid, $\mathbf{u}=\mathbf{u}(\mathbf{m},\omega,x_{s})$ is the discrete
wavefield, $\mathbf{m}>0$ is the model for the inverse of the squared wave
velocity, $\omega$ is the frequency, and the symbol $\odot$ denotes the
Hadamard product. The source, $\mathbf{q}_{s}$ is assumed to be a
discretization of a delta function, $\delta(x-x_{s})$, that is located at
$x_{s}$. The equation is discretized on a finite domain and is accompanied
with absorbing boundary conditions on all sides of the domain, and possible
attenuation term [27]. At high wave-numbers, the Helmholtz linear system (2.1)
is very challenging to solve numerically—it is indefinite, highly ill-
conditioned, and it requires a very fine mesh to accurately capture the wave
behavior. To solve FWI, we require multiple solutions of (2.1), which makes
the solution of the inverse problem highly expensive and cumbersome [41] in
addition to the other difficulties mentioned earlier.

In FWI, sources are located at many locations on the top part of the grid, and
the waveform that is generated by each source is recorded at locations where
receivers are placed. In our formulation, each observed data sample
corresponding to a source and frequency $\omega$, is given by

(2.2)
$\mathbf{d}^{{\sf{{}_{obs}}}}(\omega,x_{s})=\mathbf{P}_{s}^{\top}\mathbf{u}(\mathbf{m}_{{\sf{{}_{true}}}},\omega,x_{s})+\epsilon$

where $\mathbf{P}_{s}$ is a sampling matrix that measures the wave field
$\mathbf{u}$ that is generated by the source at $x_{s}$ at the locations of
receivers. The data contain noise $\epsilon$, which we assume to be i.i.d,
Gaussian, and with zero mean and (co-)variance $\Sigma$. Given data for many
sources and several frequencies, we aim to estimate the true model,
$\mathbf{m}_{{\sf{{}_{true}}}}$, by minimizing the difference between the
measured and simulated data, obtained by solving (2.1) across all the sources
and frequencies. This may be done by solving the PDE-constrained optimization
problem

(2.3)
$\displaystyle\min_{\begin{subarray}{c}\mathbf{m}_{L}\leq\mathbf{m}\leq\mathbf{m}_{H}\\\
\\{\mathbf{u}_{sj}\\}_{j=1,s=1}^{n_{f}\;\;n_{s}}\end{subarray}}\Phi_{\sf{{}_{Constrained}}}(\mathbf{m},\\{\mathbf{u}_{sj}\\})$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n_{f}}{\sum_{s=1}^{n_{s}}{\left\|\mathbf{P}_{s}^{\top}\mathbf{u}_{sj}-\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}\right\|^{2}_{\Sigma_{sj}^{-1}}}}+\alpha
R(\mathbf{m})$ $\displaystyle\mbox{s.t.
}\mathcal{H}(\mathbf{m},\omega_{j})\mathbf{u}_{sj}$ $\displaystyle=$
$\displaystyle\mathbf{q}_{s},\quad\quad\quad\quad\quad s=1,\ldots,n_{s},\ \
j=1,\ldots,n_{f}$

where $\mathbf{u}_{sj}$ is the waveform for source $s$ and frequency
$\omega_{j}$, which is predicted for a given model $\mathbf{m}$, according to
the forward problem (2.1). The data terms
$\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}=\mathbf{d}^{{\sf{{}_{obs}}}}(\omega_{j},x_{s})$
are the corresponding observed data as in (2.2).

We use the upper and lower bounds $\mathbf{m}_{H}>\mathbf{m}_{L}>0$, and a
regularization term $R(\mathbf{m})$, which is accompanied by a parameter
$\alpha>0$. These promote prior information in the inversion, and help us find
a reasonable solution to the otherwise ill-posed problem [47, 38]. To this
end, we assume $\mathbf{m}$ to be a layered model, and choose $R$ to promote
piece-wise smooth functions like the total variation regularization term [33].
We elaborate on these choices later in the results section.

Since we have many sources and several frequencies, the solution of the
inverse problem at large scales requires parallel/distributed software and
resources. Typically, the solution will be distributed to workers according to
frequencies, sources, and even sub-domains, and some components (e.g., the
gradient for $\mathbf{m}$) will be assembled on a master worker or a small
group of workers. For more information, and an open-source code, see [34] and
references therein.

###  2.1 The method of frequency continuation

As stated before, reaching local minima is a major problem in solving FWI. The
first step in trying to avoid local minima is adopting a method called
frequency continuation, which has proven to be effective in [28, 27].
Frequency continuation is obtained by first approximately solving the problem
(2.3) (or one of the equivalent formulations that we show later) using the
lower frequencies only in its cost function, to build a smooth approximation
of the velocity model. We then add more and more frequencies and approximately
solve the problems, each time starting from the previously obtained model,
until we cover all of the frequencies. Algorithm 1 summarizes the “window-
wise” frequency continuation approach that we use in this work, where we
consider a window of frequencies at each time, instead of all of them.
Furthermore, in this work we use the Gauss-Newton (GN) method to approximately
minimize the corresponding $\Phi$ in every frequency continuation iteration
[28]. We briefly present the projected GN later, and refer the reader to [34,
41] for a more detailed description of GN using similar notation. Projected
Conjugate Gradients (CG) is used for the inner Newton problems.

# Assume $\omega_{1}<...<\omega_{n_{f}}$.

# $ws$: Window size of frequencies that we work on each time.

# $i_{start}$, $i_{end}$ \- Initial and final frequencies to consider.
$i_{start}\geq 1$, $i_{end}\leq n_{f}$.

Initialize $\mathbf{m}$ by some reference model (or have it from a previous
cycle).

for  _$i=i_{start}:i_{end}$_ do

Approximately solve $\Phi$ using data for
$\omega_{max\\{i-ws,1\\}},...,\omega_{i}$, starting from a previous model
$\mathbf{m}$.

end for

Algorithm 1 Frequency continuation

While having good results, the frequency continuation has a major obstacle -
starting with a good initial guess $\mathbf{m}_{ref}$ is critical for an
accurate estimation. Finding such a $\mathbf{m}_{ref}$ is not a trivial task,
since our data does not correspond to low enough frequencies because of
acquisition limitations. Some of the previously mentioned methods aim exactly
at finding the initial guess, while in this work we offer a more general
approach to try to escape local minima efficiently.

###  2.2 All-at-once methods

One of the main approaches to solve (2.3) considers $\mathbf{m}$ and all the
fields $\mathbf{u}_{sj}$ as variables, and handles the PDE constraints using
the Lagrange multipliers approach, adding a Lagrange variable to each pair
$s,j$ [17, 25]. However, in real life 3D scenarios, solving the problem is
computationally demanding both in terms of memory and computations. The
forward problem (2.1) requires a very fine mesh of hundreds of millions of
grid points for each $\mathbf{u}_{sj}$. The inverse problem (2.3) typically
includes several hundreds of sources $\mathbf{q}_{s}$ and tens of frequencies
$\omega_{j}$. Handling the additional $\mathbf{u}_{sj}$ variables imposes
severe memory issues at large scales, since these vectors are huge, and there
are many of them. Manipulating them iteratively in the optimization process is
quite cumbersome at such large scales.

###  2.3 The reduced and penalized formulations

Another common approach to solve (2.3) is by eliminating the PDE constraints,
setting

(2.4)
$\mathbf{u}_{sj}=\mathbf{u}_{sj}(\mathbf{m})=\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{q}_{s}.$

This results in the reduced unconstrained (with respect to the PDEs)
optimization problem

(2.5)
$\displaystyle{\min_{\mathbf{m}_{L}\leq\mathbf{m}\leq\mathbf{m}_{H}}}\Phi_{{\sf{{}_{Reduced}}}}(\mathbf{m})=\sum_{j=1}^{n_{f}}{\sum_{s=1}^{n_{s}}{\left\|\mathbf{P}_{s}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{q}_{s}-\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}\right\|^{2}_{\Sigma_{sj}^{-1}}}}+\alpha
R(\mathbf{m})\\\ $

for the model $\mathbf{m}$ only. This problem is equivalent to (2.3), but it
imposes different solution techniques. In particular, using traditional
methods like GN [28], it is possible to handle large instances of (2.5) by
occasional use of the disk [41], instead of frequently manipulating the many
large vectors $\mathbf{u}_{sj}$.

However, while the unconstrained (2.5) is more feasible to handle than (2.3)
at large scales, it is highly non-linear. This non-linearity partially stems
from the elimination of the PDE constraints—the optimization process is forced
to fulfill those constraints at all times. Consequently, [46] suggested to
relax these constraints. That is, to allow the fields $\mathbf{u}_{sj}$ not to
satisfy the PDE (2.1), which is weakly enforced by a penalty term using the
$L_{2}$ norm. The resulting constrained problem is given by

$\displaystyle\min_{\begin{subarray}{c}\mathbf{m}_{L}\leq\mathbf{m}\leq\mathbf{m}_{H}\\\
\\{\mathbf{u}_{sj}\\}_{j=1,s=1}^{n_{f}\;\;n_{s}}\end{subarray}}\Phi_{\sf{{}_{Penalized}}}(\mathbf{m},\\{\mathbf{u}_{sj}\\})$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{n_{f}}{\sum_{s=1}^{n_{s}}{\left\|\mathbf{P}_{s}^{\top}\mathbf{u}_{sj}-\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}\right\|^{2}_{\Sigma_{sj}^{-1}}}}+$
$\displaystyle\beta\sum_{j=1}^{n_{f}}{\sum_{s=1}^{n_{s}}{\left\|\mathcal{H}(\mathbf{m},\omega_{j})\mathbf{u}_{sj}-\mathbf{q}_{s}\right\|_{2}^{2}}}+\alpha
R(\mathbf{m}).$

This problem has penalty terms which essentially replace the PDE constraints,
and $\beta>0$ is a penalty parameter that controls how accurately the
constraints are fulfilled. This $L_{2}$ penalty method was further developed
(along with suitable handling of bound-constraints and regularization) in [1,
2, 3] using the alternating direction method of multipliers (ADMM), which
relieves the user from choosing a large $\beta$ to fulfil the constraints in
high accuracy. These approaches, like the all-at-once method to solve (2.3),
include the iterative updates of the vectors $\\{\mathbf{u}_{sj}\\}$, which
require a lot of memory and are almost impractical in certain scenarios.

###  2.4 The formulation of extended sources

A different approach, albeit in the same spirit as the penalty approach in the
context of this work, enlarges the search space of the problem by extending
the point-sources $\mathbf{q}_{s}$ [20]. Even though this work suggested the
approach in a time-domain formulation, the equivalent formulation in frequency
domain may result in the problem

$\displaystyle\min_{\begin{subarray}{c}\mathbf{m}_{L}\leq\mathbf{m}\leq\mathbf{m}_{H}\\\
\\{\mathbf{z}_{s}\\}_{s=1}^{n_{s}}\end{subarray}}\Phi_{{\sf{{}_{ExtSrc}}}}(\mathbf{m},\\{\mathbf{z}_{s}\\}_{s=1}^{n_{s}})$
$\displaystyle=$
$\displaystyle\sum_{s=1}^{n_{s}}{\sum_{j=1}^{n_{f}}{\left\|\mathbf{P}_{s}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{z}_{s}-\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}\right\|^{2}_{\Sigma_{sj}^{-1}}}}+$
$\displaystyle\beta\sum_{s=1}^{n_{s}}\|\mathbf{z}_{s}\|_{W_{s}}^{2}+\alpha
R(\mathbf{m}).$

Here, the extended sources $\mathbf{z}_{s}$ are introduced as variables, and
are penalized by a weighted $\ell_{2}$ norm. The weight matrices $W_{s}$ do
not penalize $\mathbf{z}_{s}$ at the sources locations $x_{s}$, hence setting
the extended sources to be the original ones, i.e.,
$\mathbf{z}_{s}=\mathbf{q}_{s}$, does not penalize the objective. This is
another form of enlarging the search space and relaxing the PDE constraints.
Similarly to before, the variables $\mathbf{z}_{s}$ are large and typically we
have too many of them in 3D. We also note that here, we need to choose $\beta$
sufficiently large so that the extension of the sources vanishes at the end of
the optimization. This method showed very promising results, but is expensive.
In this work we wish to make it more applicable at large scales. We note that
compared to (2.3) or [46], we have only $n_{s}$ unknown vectors
$\mathbf{z}_{s}$ in addition to $\mathbf{m}$, instead of $n_{s}\cdot n_{f}$
vectors in the fields $\mathbf{u}_{sj}$. The downside of the extended sources
approach compared to (2.3) or [46], is that the forward problem (2.1) needs to
be solved repeatedly for the source variables $\mathbf{z}_{s}$.

###  2.5 Simultaneous sources via split formulation

As mentioned before, one of the most effective ways to reduce the
computational cost of solving inverse problems is by simultaneous sources,
which is obtained by trace estimation [18]. That is, given a matrix
$\mathbf{A}\in\mathbb{R}^{m\times n}$, we can approximately calculate its
Frobenius norm by [7]:

(2.8)
$\|\mathbf{A}\|_{F}^{2}=trace(\mathbf{A}^{\top}\mathbf{A})=\mathbb{E}_{x}\|\mathbf{A}\mathbf{x}\|_{2}^{2}\approx\frac{1}{p}\sum_{i=1}^{p}{\|\mathbf{A}\mathbf{x}_{i}\|_{2}^{2}}$

where each $\mathbf{x}_{i}\in\mathbb{R}^{n}$ is chosen from Radermacher
distribution—each element is randomly chosen from the set $\\{-1,1\\}$ with
equal probability. We can reformulate in matrix notation:
$\frac{1}{p}\sum_{i=1}^{p}{\|\mathbf{A}\mathbf{x}_{i}\|_{2}}=\frac{1}{p}\|\mathbf{A}\mathbf{X}\|_{F}^{2}$
where $\mathbf{X}$ is a $n\times p$ matrix whose column $i$ is
$\mathbf{x}_{i}$. In practice we calculate a norm of an $m\times p$ matrix
instead of an $m\times n$ matrix, which can save significant computations when
$p\ll n$.

To apply the simultaneous sources technique, we wish to approximate the data
misfit terms in formulations like (2.3) or (2.5), by applying (2.8). However,
this is possible only when the operators $\mathbf{P}_{s}$ in these
formulations do not depend on $s$, i.e., when all the receivers record the
waveforms from all the sources. The same goes for the error (co-)variance
matrices $\Sigma_{sj}$. Such requirements are not met in many realistic
scenarios. The work of [23] handled this issue using a quite general approach:
by splitting the data term. In this approach a new set of data variables
$\hat{\mathbf{d}}_{sj}$ is introduced, and is defined on the union of the
supports of $\mathbf{P}_{s}$ for all $s$—that is the union of the locations of
all the receivers for all the sources. Next, $\hat{\mathbf{d}}_{sj}$ should be
similar to the observed data $\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}$ at the
locations of the receivers that actually record the waveform for each source
$s$. To this end, (2.5), for example, is reformulated as

$\displaystyle\displaystyle{\min_{\begin{subarray}{c}\mathbf{m}_{L}\leq\mathbf{m}\leq\mathbf{m}_{H}\\\
\\{\hat{\mathbf{d}}_{sj}\\}_{j=1,s=1}^{n_{f}\;\;n_{s}}\end{subarray}}}\Phi_{\sf{{}_{ReducedSplit}}}(\mathbf{m},\\{\hat{\mathbf{d}}_{sj}\\})=$
$\displaystyle\quad\quad\quad\sum_{j=1}^{n_{f}}{\sum_{s=1}^{n_{s}}{\left\|\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{q}_{s}-\hat{\mathbf{d}}_{sj}\right\|^{2}_{\Sigma_{j}^{-1}}+\eta\|\hat{\mathbf{P}}_{s}^{\top}\hat{\mathbf{d}}_{sj}-\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}\|^{2}_{\Sigma_{sj}^{-1}}}}+\alpha
R(\mathbf{m}),$

where $\mathbf{P}$ is an operator that projects a vector $\mathbf{u}$ onto the
union of the receivers’ locations, and $\hat{\mathbf{P}}_{s}$ is defined to
choose the subset of receivers for source $s$ out of that union (that is, such
that $\mathbf{P}_{s}^{\top}=\hat{\mathbf{P}}_{s}^{\top}\mathbf{P}^{\top}$).
$\eta>0$ is a data fitting parameter that needs to be chosen relatively high
to ensure $\hat{\mathbf{d}}_{sj}$ will be similar to
$\mathbf{d}^{{\sf{{}_{obs}}}}_{sj}$ at the support of the receivers of each
source $s$. Minimizing (2.5) with respect to $\hat{\mathbf{d}}_{sj}$ is
trivial, and minimizing it with respect to $\mathbf{m}$ is similar to (2.5),
but can be applied by stochastic trace estimation more efficiently. Hence,
techniques such as alternating minimization or variable projection [6] are
highly favorable here. Adding a regularization to $\hat{\mathbf{d}}_{sj}$ like
in [23] might improve this formulation and is also a subject of research.

To apply the simultaneous sources technique to (2.5), we first select the size
of the subspace we want (denoted by $p$), then define a random matrix
$\mathbf{X}$ of size $n_{s}\times p$ from Rademacher distribution. By (2.8) we
obtain

$\displaystyle\sum_{s}\|\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{q}_{s}-\hat{\mathbf{d}}_{s,j}\|_{\Sigma_{j}^{-1}}^{2}$
$\displaystyle=$
$\displaystyle\|\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{Q}-\hat{\mathbf{D}}_{j}\|_{\Sigma_{j}^{-1}}^{2}$
$\displaystyle\approx$
$\displaystyle\frac{1}{p}\|\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{Q}\mathbf{X}-\hat{\mathbf{D}}_{j}\mathbf{X}\|_{\Sigma_{j}^{-1}}^{2}$

where $\mathbf{Q}$ is the matrix whose columns are all the sources
$\mathbf{q}_{s}$, and $\hat{\mathbf{D}}_{j}$ is the matrix of all data
unknowns $\hat{\mathbf{d}}_{sj}$. This effectively reduces the number of
sources to $p$ in every iteration. Following this splitting, we henceforth
present our methods assuming that $\mathbf{P}_{s}=\mathbf{P}$, and
$\Sigma_{sj}=\Sigma_{j}$ for all the sources $s$, and ignore the extra data
term in (2.5).

##  3 Robust FWI with low-rank and sparse extended sources

In this work we aim to relax the PDE-constraints without introducing full
(memory consuming) variables such as $\mathbf{u}_{sj}$, aiming for a problem
that is similar to the more memory-friendly (2.5) rather than to (2.3) or
(2.4). To this end, we consider source extensions $\mathbf{z}_{s}$ as
variables. In our formulation, the new sources are
$\mathbf{q}_{s}+\mathbf{z}_{s}$ (or $\mathbf{Q}+\mathbf{Z}$ in matrix
notation), and we aim that at the end of the reconstruction, $\mathbf{z}_{s}$
will vanish. However, since we have many sources, this again introduces a
significant amount of extra variables, like (2.4), that we wish to prevent. To
achieve that, we couple between all the extended sources, and define all of
them as a low-rank matrix. Also, we use the $\ell_{1}$ norm as penalty, and
aim that the source extensions are sparse to further reduce the memory
consumption. Because the vectors $\mathbf{q}_{s}$ are sparse (they are a
discrete $\delta$ function), they do not require a lot of memory—we wish the
same for the extensions $\mathbf{z}_{s}$. The $\ell_{1}$ norm gives a higher
penalty for small non-zero numbers compared to the $\ell_{2}$ norm, and
encourages sparse results with a high amount of zeros.

More explicitly, we use the matrix notation of (2.5), and define the source
extension matrix as $\mathbf{Z}=\mathbf{Z}_{1}\mathbf{Z}_{2}$, where
$\mathbf{Z}_{1}\in\mathbb{C}^{N\times n_{es}}$ and
$\mathbf{Z}_{2}\in\mathbb{C}^{n_{es}\times n_{s}}$ are the matrices that form
the low-rank decomposition of $\mathbf{Z}$. $N$ is the total number of grid
nodes in the domain, $n_{s}$ is the number of original sources and
$n_{es}<n_{s}$ is the maximal rank of the extended sources matrix $\mathbf{Z}$
(essentially the number of extended sources). $n_{es}$ is chosen to be small
enough such that $\mathbf{Z}_{1}$ can reasonably fit in memory, and the
computations involving $\mathbf{Z}_{1}$ are reasonable. On the other hand,
$\mathbf{Z}_{1}$ needs to be rich (wide) enough to allow the relaxation of the
PDE-constraints. This is a balance that we need to manage. Our optimization
problem (in matrix notation) becomes:

(3.1)
$\displaystyle\quad\min_{\begin{subarray}{c}\mathbf{m}_{L}\leq\mathbf{m}\leq\mathbf{m}_{H}\\\
\mathbf{Z}_{1},\mathbf{Z}_{2}\end{subarray}}\Phi_{{\sf{{}_{LowRankExtSrc}}}}(\mathbf{m},\mathbf{Z}_{1},\mathbf{Z}_{2})=$
$\displaystyle\quad\quad\quad\sum_{j=1}^{n_{f}}{\left\|\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}(\mathbf{Q}+\mathbf{Z}_{1}\mathbf{Z}_{2})-\mathbf{D}^{{\sf{{}_{obs}}}}_{j}\right\|^{2}_{\Sigma_{j}^{-1}}}+\beta_{1}\|\mathbf{Z}_{1}\|_{1}+\frac{\beta_{2}}{2}\|\mathbf{Z}_{2}\|_{F}^{2}+\alpha
R(\mathbf{m}).$

$\beta_{1},\beta_{2}>0$ are regularization parameters that controls how much
we allow the $\mathbf{Z}$ matrices to change and reduce the objective. The
penalty term for $\mathbf{Z}_{2}$ is the Frobenius norm, and the penalty for
$\mathbf{Z}_{1}$ is the element-wise $\ell_{1}$ norm to promote sparsity in
$\mathbf{Z}_{1}$. This stems from memory considerations to help cases where
the required rank of $\mathbf{Z}$ needs to be high. This way, the PDE-
constraints are relaxed, since we remove the demand to strictly satisfy the
PDE constraints for the original sources $\mathbf{Q}$. The advantage is that
we achieve that without increasing the required memory by a lot. The closely
related extended sources approach of [20], presented in (2.4) uses a full rank
$\mathbf{Z}$ and a weighted $\ell_{2}$ norm instead of an $\ell_{1}$ norm.
These differences have a significant practical implication in 3D, as the low-
rank structure and sparsity of $\mathbf{Z}$ are beneficial because of
computational and memory considerations.

###  3.1 Iterative solution by alternating minimization

We solve the objective (3.1) by alternating minimization (ALM), where we
alternate between the approximate minimization of $\mathbf{m}$,
$\mathbf{Z}_{1}$, and $\mathbf{Z}_{2}$ in turns. Minimizing for $\mathbf{m}$
is equivalent to minimizing (2.5), and can be obtained by the LBFGS [29] or GN
methods. Minimizing for $\mathbf{Z}_{2}$ can be obtained directly, as it is a
minimization of a rather small quadratic function (details are given later).
Minimizing for $\mathbf{Z}_{1}$ can be done efficiently using iterative re-
weighted least squares (IRLS) or proximal methods (in particular, proximal CG
[48, 42]). Alg. 2 presents the ALM algorithm which is integrated into the
standard frequency continuation procedure in Alg. 1.

#Assume $\omega_{1}<...<\omega_{n_{f}}$.

# $ws$: Window size of frequencies that we work on each time.

# $i_{start}$, $i_{end}$ \- Initial and final frequencies to consider.
$i_{start}\geq 1$, $i_{end}\leq n_{f}$.

Initialize $\mathbf{m}^{(0)}$ by some reference model.

Initialize $\mathbf{Z}_{1}$ randomly (if there is no previous solution).

for  _$i=i_{start}:i_{end}$_ do

for  _$j=1:iter_{ALM}$_ do

  1. 1.

Solve for $\mathbf{Z}_{2}$ given $\mathbf{Z}_{1}$.

  2. 2.

Apply CG iterations for $\mathbf{Z}_{1}$.

  3. 3.

Solve for $\mathbf{Z}_{2}$ given $\mathbf{Z}_{1}$.

  4. 4.

Apply a GN iteration for $\Phi(\mathbf{m})$ using data for
$\omega_{max\\{i-ws,1\\}},...,\omega_{i}$,
starting from a previous model $\mathbf{m}$.

  5. 5.

Update $\beta_{1},\beta_{2}$ as in Section 3.2.

end for

end for

Algorithm 2 Frequency continuation with alternating minimization.

####  3.1.1 Solving for $\mathbf{m}$ using Gauss-Newton

The solution of (3.1) with respect to $\mathbf{m}$ is hardly influenced by the
extended sources, and is similar to the solution of (2.5). At each iteration
$k$ of GN we obtain the linear approximation

(3.2)
$\mathbf{u}_{s,j}(\mathbf{m}^{(k)}+\delta\mathbf{m})\approx\mathbf{u}_{s,j}(\mathbf{m}^{(k)})+\mathbf{J}_{s,j}(\mathbf{m}^{(k)})\delta\mathbf{m},$

where
$\mathbf{u}_{sj}(\mathbf{m})=\mathcal{H}(\mathbf{m},\omega_{j})^{-1}(\mathbf{q}_{s}+\mathbf{z}_{s})$,
and $\mathbf{J}_{s,j}=\nabla_{\mathbf{m}}\mathbf{u}_{s,j}$ is the Jacobian
matrix of the data for every source $s$ and frequency $\omega_{j}$.
$\mathbf{z}_{s}$ is the $s$-th column of $\mathbf{Z}_{1}\mathbf{Z}_{2}$. In
real-life scales, the Jacobian matrix cannot be stored in memory [16], but we
can apply matrix-vector products with this matrix and its conjugate transpose.

At each step, we place (3.2) in (3.1) and solve a quadratic minimization by an
iterative solver where only matrix vector products of the Jacobians are
computed. To obtain the GN step $\delta\mathbf{m}$, we first compute the
gradient of (3.1), which is given by

(3.3)
$\nabla_{\mathbf{m}}\Phi(\mathbf{m}^{(k)})=\sum_{s,j}\mathbf{J}_{sj}(\mathbf{m}^{(k)})^{\top}\mathbf{P}(\mathbf{P}^{\top}\mathbf{u}_{sj}(\mathbf{m}^{(k)})-\mathbf{d}_{sj}^{{\sf{{}_{obs}}}})+\alpha\nabla_{\mathbf{m}}R(\mathbf{m}^{(k)}),$

and then we approximately solve the linear system
$\mathbf{H}\delta\mathbf{m}=-\nabla_{\mathbf{m}}\Phi(\mathbf{m}^{(k)}),$ where
the Gauss-Newton Hessian is defined by

(3.4)
$\mathbf{H}=\sum_{s,j}\mathbf{J}_{sj}(\mathbf{m}^{(k)})^{\top}\mathbf{P}\mathbf{P}^{\top}\mathbf{J}_{sj}(\mathbf{m}^{(k)})+\alpha\nabla^{2}R(\mathbf{m}^{(k)}).$

Once the linear system is approximately solved, the model is updated,
$\mathbf{m}\leftarrow\mathbf{m}+\mu\delta\mathbf{m}$ where $\mu$ is a line
search parameter that is chosen such that the objective function is
sufficiently decreased at each iteration (the Armijo rule).

The Jacobian matrix required in (3.3)-(3.4), is given by

(3.5)
$\displaystyle\mathbf{J}_{sj}(\mathbf{m})=-\omega^{2}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\,\mbox{diag}(\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{z}_{s}).$

Assuming that the fields $\mathbf{u}_{sj}$ are stored in memory (or the disk),
the multiplication of this matrix with a vector requires one forward solution
(for each pair of source and frequency). The fields can be stored in a rather
low precision [41].

####  3.1.2 Solving for $\mathbf{Z}_{2}$

The minimization of (3.1) with respect to $\mathbf{Z}_{2}$ can be obtained
directly, as this is a small-sized quadratic minimization problem. Denote the
residual $\mathbf{R}_{j}$ and temporary matrix $\mathbf{T}_{j}$ by

(3.6)
$\mathbf{R}_{j}=\mathbf{D}_{j}^{{\sf{{}_{obs}}}}-\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\mathbf{Q},\quad\quad\mathbf{T}_{j}=\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}\in\mathbb{C}^{n_{r}\times
N},$

where $n_{r}$ is the number of receivers. The solution for $\mathbf{Z}_{2}$ is
given by solving $n_{s}$ small linear systems of size $n_{es}\times n_{es}$
with the same symmetric and positive definite matrix:

(3.7)
$\mathbf{Z}_{2}=\left(\sum_{j}(\mathbf{T}_{j}\mathbf{Z}_{1})^{*}\Sigma_{j}^{-1}(\mathbf{T}_{j}\mathbf{Z}_{1})+\beta_{2}\mathbf{I}\right)^{-1}\sum_{j}(\mathbf{T}_{j}\mathbf{Z}_{1})^{*}\Sigma_{j}^{-1}\mathbf{R}_{j}.$

The residuals $\mathbf{R}_{j}$ and matrices $(\mathbf{T}_{j}\mathbf{Z}_{1})$
are computed as part of the update for $\mathbf{Z}_{1}$, which are both kept
fixed in this part. We note that $\mathbf{Z}_{2}$ is easily computed given
$\mathbf{Z}_{1}$, hence, there is no need to keep its iterative state or
initialize it.

####  3.1.3 Solving for $\mathbf{Z}_{1}$

When limited to $\mathbf{Z}_{1}$ only, problem (3.1) is a quadratic
minimization with an $\ell_{1}$ penalty, which is also called the least
absolute shrinkage and selection operator (LASSO) regression. This problem was
originally suggested in [35] for seismic inversion like here, but went on to
be highly popular in other applications, mainly signal processing [48]. It is
a well understood problem with a variety of available solvers, like IRLS,
proximal CG [42] or SEquential Subspace OPtimization (SESOP) [48]. These three
are the most suitable methods in our context, because they allow us to apply
the relevant matrices as “black-box” operators, some of which include the
solution of (2.1).

Again, using the residual $\mathbf{R}_{j}$ and $\mathbf{T}_{j}$ in (3.6), the
problem for $\mathbf{Z}_{1}$ is given by:

(3.8)
$\displaystyle\quad\quad\sum_{j=1}^{n_{f}}{\left\|\mathbf{T}_{j}\mathbf{Z}_{1}\mathbf{Z}_{2}-\mathbf{R}_{j}\right\|^{2}_{\Sigma_{j}^{-1}}}+\beta_{1}\|\mathbf{Z}_{1}\|_{1}.$

Because of the coupling, we get a single linear system inside the $\ell_{2}$
norm for all the unknowns in $\mathbf{Z}_{1}$. That is, this is not a block
linear system with multiple right-hand-sides, like (3.7). Using the $\otimes$
symbol for the Kronecker product, the matrix equation in (3.8) can also be
vectorized as:

(3.9)
$\displaystyle\quad\quad\sum_{j=1}^{n_{f}}{\left\|(\mathbf{Z}_{2}^{*}\otimes\mathbf{T}_{j}){\sf{vec}}(\mathbf{Z}_{1})-{\sf{vec}}(\mathbf{R}_{j})\right\|^{2}_{\Sigma_{j}^{-1}}}+\beta_{1}\|{\sf{vec}}(\mathbf{Z}_{1})\|_{1},$

where ${\sf{vec}}()$ denotes the column-stacking of any matrix into a vector.
In this work we solve the LASSO minimization using IRLS, where at each ALM
iteration for $\mathbf{Z}_{1}$ in Alg. (2) we replace the $\ell_{1}$ norm with
a weighted $\ell_{2}$ norm

(3.10)
$\|{\sf{vec}}(\mathbf{Z}_{1})\|_{1}\rightarrow\frac{1}{2}\|{\sf{vec}}(\mathbf{Z}_{1})\|_{\mathbf{W}}^{2},\mbox{
where
}\mathbf{W}=\mbox{diag}\\{(|{\sf{vec}}(\mathbf{Z}_{1})|+\varepsilon)^{-1}\\},$

where $|{\sf{vec}}(\mathbf{Z}_{1})|$ denotes the stacked vector of absolute
entries of $\mathbf{Z}_{1}$. This way, the gradient of the temporary IRLS
objective smoothly approximates the gradient of (3.9). Then, the IRLS
objective is approximately minimized by a few iterations of standard
preconditioned CG for the normal equations, which directly minimizes the IRLS
approximation of (3.9) in each of its iterations. For efficiency, the
operators for CG are computed in matrix form as follows

(3.11)
$\mathcal{OP}(\mathbf{Z}_{1})=\mathbf{T}_{j}\mathbf{Z}_{1}\mathbf{Z}_{2},\quad\mathcal{OP}^{*}(\mathbf{R})=\mathbf{T}_{j}^{*}\mathbf{R}\mathbf{Z}_{2}^{*},$

where the conjugate transposed operator $\mathcal{OP}^{*}(\mathbf{R})$ is
equivalent to

(3.12)
$(\mathbf{Z}_{2}^{*}\otimes\mathbf{T}_{j})^{*}{\sf{vec}}(\mathbf{R})=(\mathbf{Z}_{2}\otimes\mathbf{T}_{j}^{*}){\sf{vec}}(\mathbf{R})={\sf{vec}}(\mathbf{T}_{j}^{*}\mathbf{R}\mathbf{Z}_{2}^{*})$

in vectorized form. We apply a few such CG iterations (specifically, about 5
each time), as an approximate minimization for $\mathbf{Z}_{1}$ inside the
most inner ALM algorithm. As a preconditioner, we use the matrix $\mathbf{W}$
in (3.10).

The solution for $\mathbf{Z}_{1}$ and $\mathbf{Z}_{2}$ and update for
$\mathbf{m}$ are repeated $iter_{ALM}$ times in Algorithm (2). In this
process, the LU factorization or preconditioner of $\mathcal{H}$ needed to
multiply the matrices $\mathbf{T}_{j}$ with vectors, is obtained as part of
the iteration for $\mathbf{m}$.

###  3.2 The choice of the regularization parameters

The regularization parameters $\beta_{1}$ and $\beta_{2}$ play a significant
role in the extended sources framework. They control how much we keep the
problem (3.1) close to the original one (2.5). On the one hand, we wish to
allow the matrix $\mathbf{Z}=\mathbf{Z}_{1}\mathbf{Z}_{2}$ to have significant
enough values to influence the optimization process. On the other hand, we
wish to keep this process close to the reality, where the extensions do not
exist, since they are only artificial. The work [20] chose the parameter
$\beta$ in (2.4) to be fixed, so initially, when the misfit is high with
respect to $\mathbf{m}$ this parameter is relatively low, and as $\mathbf{m}$
improves with the iterations, $\beta$ becomes more significant with respect to
$\mathbf{Z}$. As one can expect, we observed a similar behavior in our
formulation (3.1) with $\beta_{1}$ and $\beta_{2}$. As a rule of thumb, the
work [20] chose $\beta$ such that the misfit with the extended sources is
about half of the misfit without the extended sources (if setting
$\mathbf{Z}=0$). This is obviously a quite strong penalty, as the full
$\mathbf{Z}$ can easily zero out the misfit in (3.1) almost for any
$\mathbf{m}$, if not penalized. That is part of the motivation for this work.
By setting a high-enough $\beta$, the work [20] essentially limits
$\mathbf{Z}$. Here, we limit it by using a low-rank structure that is much
more favorable computationally, even though it does introduce some algorithmic
complications, as discussed above. In our framework, we keep the ratio between
$\beta_{1}$ and $\beta_{2}$ fixed—specifically, we choose
$\beta_{2}=100\beta_{1}$, which was chosen based on trial and error.
Throughout the iterations, we change $\beta_{i}$ together to keep the ratio
between misfits (with and without extended sources) to be approximately 0.5.
More explicitly, to keep the ratio in the section $[r_{1},r_{2}]$ we apply the
following rule:

(3.13) If
$\displaystyle\frac{\mbox{misfit}(\mathbf{Z}_{1}\mathbf{Z}_{2})}{\mbox{misfit}(0)}>r_{2}\mbox{
then }\beta_{1},\beta_{2}\leftarrow\beta_{1}/\gamma,\beta_{2}/\gamma,$ Else if
$\displaystyle\frac{\mbox{misfit}(\mathbf{Z}_{1}\mathbf{Z}_{2})}{\mbox{misfit}(0)}<r_{1}\mbox{
then }\beta_{1},\beta_{2}\leftarrow\beta_{1}\cdot\gamma,\beta_{2}\cdot\gamma.$

In this work we choose $r_{1}=0.3$, $r_{1}=0.5$, and $\gamma=1.5$.

###  3.3 Computational costs

The cost of the entire inversion is dominated by 1) the cost of the GN
iterations for $\mathbf{m}$, and 2) the cost of the IRLS-CG iterations for
$\mathbf{Z}_{1}$. The cost of the minimization for $\mathbf{Z}_{2}$ is quite
negligible compared to the cost of the other two. Below we provide details
regarding each of these components.

####  3.3.1 The cost of GN (standard reduced-formulated FWI)

The cost of each GN iteration is entirely dominated by the costs of the
forward solvers, and that is also the dominant cost in standard reduced FWI in
(2.5). At each gradient and Hessian-vector multiplication, we need to solve
the forward problem (2.1) twice: once to compute the Jacobian in (3.5), and
once for its adjoint, where the adjoint Helmholtz equation is solved (that is
assuming that the fields $\mathbf{u}_{sj}$ are stored). Given
$\mathbf{m}^{(k)}$, if possible (e.g., in 2D) we can factorize the matrices
$\mathcal{H}(\omega_{j})$ for all frequencies, using a direct solver [4, 36],
and then only need to apply the forward-backward substitutions for each
source. Solutions in 3D typically require using iterative solvers with
effective preconditioners (e.g., [13, 26]) to be computationally efficient.
The preconditioner often dominates the computational cost. The preconditioner
setup, like the LU factorization, is applied only once per iteration. It is
clear that the cost of GN is controlled by the number of frequencies we
consider in a frequency continuation window, and more importantly, the number
of sources. To summarize:

(3.14)
$\mbox{cost(GN)}=\sum_{frequencies}\mbox{cost(LinSetup)}+\\#iter(\mbox{GN.CG})\cdot
2\cdot n_{s}\cdot\mbox{cost(LinSolve)},$

where $\\#iter(\mbox{GN.CG})$ is the number of inner CG iterations in GN.
Setting up the sources, i.e., multiply $\mathbf{Z}_{1}\mathbf{Z}_{2}$, may be
costly if the rank of source extension $n_{es}$ is not low as we choose here.

In terms of memory, the footprint of the GN iterations is given by

(3.15)
$\mbox{mem(GN)}=\sum_{frequencies}\mbox{mem(LinSetup)}+n_{s}\cdot\mathcal{O}(N),$

where mem(LinSetup) is the memory requirements of the Helmholtz solver, and
$N$ is the forward mesh size. The $\mathcal{O}(N)$ storage mainly reflects the
storage of the fields $\mathbf{u}_{sj}$ that are needed in the sensitivity
computations in Eq. (3.5). If the fields are not stored, the multiplication of
the sensitivities with vector requires twice the number of forward
simulations.

####  3.3.2 The cost of the minimization for $\mathbf{Z}_{1}$

The cost of the minimization for $\mathbf{Z}_{1}$ _per extended source_ is
similar to that of GN. The matrices $\mathbf{T}_{j}$ that are defined in (3.6)
and needed in (3.11) typically cannot be stored in memory. Hence, for example,
to apply $\mathcal{OP}$ to $\mathbf{Z}_{1}$ in (3.11), we need to apply a
forward simulation for each pair of extend source and frequency. This shows
that the low-rank structure of the source extension
$\mathbf{Z}_{1}\mathbf{Z}_{2}$ is crucial for the method to be computationally
attractive. Note that the solution for $\mathbf{Z}_{1}$ does not require the
solver or preconditioner setups for the Helmholtz operators, and can reuse the
ones computed as part of GN. To summarize:

(3.16) $\mbox{cost(}\mathbf{Z}_{1}\mbox{
solve)}=\sum_{frequencies}\\#iter(\mbox{IRLS.CG})\cdot 2\cdot
n_{es}\cdot\mbox{cost(LinSolve)},$

where $\\#iter(\mbox{IRLS.CG})$ is the number of inner CG iterations in IRLS.
In terms of memory, the footprint of the $\mathbf{Z}_{1}$ minimization is
given by

(3.17) $\mbox{mem(}\mathbf{Z}_{1}\mbox{
solve)}=\sum_{frequencies}\mbox{mem(LinSetup)}+n_{es}\cdot\mathcal{O}(N),$

which can be expensive if $n_{es}$ is large. In particular, in [20] the number
of extended sources equals to the number of sources $n_{s}$ as in (2.4), and
the cost is proportional to $n_{s}$ instead of $n_{es}$. In this analysis we
neglect the memory saving that we can exploit from the sparsity of
$\mathbf{Z}_{1}$, since the peak memory of the sparse solvers can reach a high
percentage of the unknowns, and should be carefully controlled. In any case,
this is only crucial if $n_{es}$ is high.

Comparison: As one can observe, the cost of the two dominant components of the
algorithm are controlled by similar factors: number of CG iterations, and
number of sources involved. If, as we expect, the total number of sources
$n_{s}$ is significantly larger than the rank of extended sources $n_{es}$,
then the additional computations and memory for the inversion following the
low-rank source extension is low.

##  4 FWI using both extended and simultaneous sources

To further ease the computational cost of the inversion we will effectively
reduce the number of sources in the misfit at each GN iteration using the
simultaneous sources technique presented in Section 2.5. Here we describe how
to combine it with the low-rank extended sources objective (3.1). We note
that, as far as we know, no work describes the combination of simultaneous
sources with the standard extended sources (2.4). The task is not
straightforward, since in some sense, the simultaneous sources technique
compactly “summarizes” the many sources $\mathbf{Z}$ into a few by
$\mathbf{Z}\mathbf{X}$. However, if all those many sources in $\mathbf{Z}$ are
part of the inversion unknowns, then it is not clear how to update them based
on their compactly estimated version, without investing the computations for
all of them. A standard update for a full $\mathbf{Z}$ costs proportionally to
$n_{s}$ forward linear solves—that is the type of computation that we wish to
prevent.

Basically, we wish to combine equations (2.5) and (3.1). Given a random matrix
$\mathbf{X}\in\mathbb{R}^{n_{s}\times p}$, the combination leads to the
objective

(4.1)
$\displaystyle\quad\;\;\quad\min_{\begin{subarray}{c}\mathbf{m}_{L}\leq\mathbf{m}\leq\mathbf{m}_{H}\\\
\mathbf{Z}_{1},\mathbf{Z}_{2}\end{subarray}}\Phi_{{\sf{{}_{LowRankExtSimSrc}}}}(\mathbf{m},\mathbf{Z}_{1},\mathbf{Z}_{2};\mathbf{X})=$
$\displaystyle\frac{1}{p}\sum_{j=1}^{n_{f}}{\left\|\mathbf{P}^{\top}\mathcal{H}(\mathbf{m},\omega_{j})^{-1}(\mathbf{Q}+\mathbf{Z}_{1}\mathbf{Z}_{2})\mathbf{X}-\mathbf{D}^{{\sf{{}_{obs}}}}_{j}\mathbf{X}\right\|^{2}_{\Sigma_{j}^{-1}}}+\beta_{1}\|\mathbf{Z}_{1}\|_{1}+\frac{\beta_{2}}{2p}\|\mathbf{Z}_{2}\mathbf{X}\|_{F}^{2}+\alpha
R(\mathbf{m}).$

Essentially, compared to (3.1) we have that $\mathbf{Q}\mathbf{X}$,
$\mathbf{D}^{{\sf{{}_{obs}}}}_{j}\mathbf{X}$, and $\mathbf{Z}_{2}\mathbf{X}$
replace $\mathbf{Q}$ and $\mathbf{D}^{{\sf{{}_{obs}}}}_{j}$, and
$\mathbf{Z}_{2}$ respectively, which is similar to the standard simultaneous
method in (2.5). By frequently changing $\mathbf{X}$ in (4.1), we can
essentially solve (3.1) at reduced cost.

###  4.1 Alternating minimization

To use the simultaneous sources approach, we choose a new dimensionality
reduction matrix $\mathbf{X}$ in step 1 of Algorithm 2, and apply the ALM
iteration with updates over $\mathbf{m}$, $\mathbf{Z}_{1}$ and
$\mathbf{Z}_{2}$ based on (4.1). We apply the following steps which are
similar to the ones described before.

To solve for $\mathbf{m}$, the GN method is applied in the same way as before,
only with the low-rank source-extension
$(\mathbf{Q}+\mathbf{Z}_{1}\mathbf{Z}_{2})\mathbf{X}$, which is easily
computed in memory thanks to the sparsity of $\mathbf{Q}$ and the low-rank
structure of $\mathbf{Z}_{1}\mathbf{Z}_{2}$. This part would be
computationally similar also if we use standard FWI with extended sources in
(2.5).

To account for $\mathbf{Z}_{2}$, we compute its counterpart
$\hat{\mathbf{Z}}_{2}=\mathbf{Z}_{2}\mathbf{X}$ by direct minimization. That
is, $\hat{\mathbf{Z}}_{2}$ replaces $\mathbf{Z}_{2}$ in (3.1), and just like
$\mathbf{Z}_{2}$ is state-less (directly computed and is not updated
iteratively), $\hat{\mathbf{Z}}_{2}$ is also state-less. Moreover, it is
computed using the same formula as (3.7), only now with the reduced residual
$\hat{\mathbf{R}}_{j}=\mathbf{R}_{j}\mathbf{X}$ instead of $\mathbf{R}_{j}$.
Hence, we do not keep track of $\mathbf{Z}_{2}$ in the inversion, and compute
its reduced version $\hat{\mathbf{Z}}_{2}$ directly every iteration, given
$\mathbf{X}$.

The solution for $\mathbf{Z}_{1}$ remains the same as the minimization of
(3.8), only with $\hat{\mathbf{Z}}_{2}$ and $\hat{\mathbf{R}}_{j}$ given above
instead of $\mathbf{Z}_{2}$ and $\mathbf{R}_{j}$, respectively. We choose the
dimension of the simultaneous sources and extended sources to be similar
($p\approx n_{es}$) so that the minimization for $\mathbf{Z}_{1}$, even though
only approximated, will not over-fit the data that is reduced by the
multiplication in $\mathbf{X}$.

###  4.2 Computational costs

Because the column-dimension of $\mathbf{Z}_{1}$ is $n_{es}$, then the cost of
the update for $\mathbf{Z}_{1}$ is similar to (3.16). As in the case of
simultaneous sources with standard FWI, the real saving is in the GN
iterations. That is, we have $p$ sources instead of $n_{s}$, so compared to
(3.14) we have

(4.2) $\mbox{cost(GN-
SS)}=\sum_{frequencies}\mbox{cost(LinSetup)}+\\#iter(\mbox{GN.CG})\cdot 2\cdot
p\cdot\mbox{cost(LinSolve)},$

where $\\#iter(\mbox{GN.CG})$ is the number of CG iterations in GN. Memory-
wise, we have

(4.3) $\mbox{mem(GN-
SS)}=\sum_{frequencies}\mbox{mem(LinSetup)}+p\cdot\mathcal{O}(N).$

In short, if $p\approx n_{es}$ the solution cost for $\mathbf{Z}_{1}$ is
proportional to the GN iterations for $\mathbf{m}$ using simultaneous sources.
This is what we wanted to achieve here.

##  5 Numerical Results

In this section we demonstrate our low-rank extended sources approach with and
without the simultaneous sources technique and compare it to standard FWI for
velocity model reconstruction. For the purpose of demonstration, we do not
augment the inversions with other modalities or techniques known in the
literature. In principle, such techniques can be used in addition to our
approach in more complicated real-life experiments. We conduct our experiments
on two 2D models, one is the SEG/EAGE salt model [5], the other is the
Marmousi model [11]. For each model we include three experiments - the first
is a reconstruction using standard FWI, which is obtained by minimizing the
reduced formulation (2.5), using a few sweeps of Algorithm 1. The second
experiment is FWI with the low-rank extended sources formulation according to
(3.1), using Algorithm 2. The last experiment involves a reconstruction using
both the extended and simultaneous sources techniques, according to (4.1),
again using Algorithm 2, but this time with a new random matrix $\mathbf{X}$
at every ALM iteration as explained in Section 4. For all the settings and
models we display the resulting reconstructed model, and the misfit at each
iteration of the GN (or ALM) methods during frequency continuation sweeps.
That misfit is computed for all the sources and frequencies, and without the
source extensions, regardless of the frequency window or method used.

In addition, we demonstrate the feasibility of our method in 3D, by performing
an experiment on part of the SEG/EAGE Overthrust model [5]. In particular, we
wish to demonstrate that a 3D FWI experiment with extended sources can be
obtained using a rather standard workstation in terms of memory and
computations.

All the code in these experiments was written in the Julia language [9], as
part of the open source jInv framework [34]. Some critical parts of the code,
like the LU solver (forward + backward substitution) and matrix-vector
products for the forward modelling is written in C++ and is parallelized using
shared memory OpenMP. The experiments were computed on a workstation with
Intel Xeon Gold 5117 2GHz X 2 (14 cores per socket) with 256 GB RAM, running
on Centos 7 Linux distribution. Our code is available online at
https://github.com/JuliaInv/jInvSeismic.jl.

###  5.1 Smoothing regularization terms used for the reconstructions

The objective functions in the FWI formulations in Equations (2.5), (3.1), and
(4.1) contain a regularization term $R(\mathbf{m})$. Based on the work [41] we
apply two regularization functions - one is a high order regularization called
spline smoothing, given by

(5.1) $R_{1}(\mathbf{m})=\|\Delta_{h}(\mathbf{m}-\mathbf{m}_{ref})\|^{2}_{2}.$

The goal of this regularization is to create a smooth model from high-
frequency data. The reference model $\mathbf{m}_{ref}$ is set to be the
initial guess for the inversion, and while using this regularization we keep
$\mathbf{m}_{ref}$ fixed. We use (5.1) to obtain a good smooth model, so that
the rest of the process will result in a plausible reconstruction.

The second regularization function is a standard diffusion regularization

(5.2) $R_{2}(\mathbf{m})=\|\nabla_{h}(\mathbf{m}-\mathbf{m}_{ref})\|^{2}_{2},$

where $\nabla_{h}$ represents a discretized gradient on a nodal grid. When
using this regularization we change $\mathbf{m}_{ref}$ at each GN/ALM
iteration to be the resulting model. Updating $\mathbf{m}_{ref}$ encourages
the model to change more at each iteration, resulting in faster convergence.
In this problem, (5.2) results in a sharp reconstruction, and is suitable to
use after an initial smooth model is constructed. Total Variation regularizer
is another approach to promote sharper piece-wise constant/smooth
reconstruction.

###  5.2 Inversion algorithm and general settings

In this section we describe the general setup of our inversion algorithms
which is common for all experiments, and more specific details will be
provided for each experiment separately. For the inversions we use the
frequency continuation strategies in Algorithms 1 or 2—we use 3-4 sweeps of
the corresponding version depending on each setting and model, with a window
size of 4 frequencies for all the sweeps. At first, we apply frequency
continuation sweeps to create an initial smooth guess, using the smoothing
high order regularization (5.1). The remaining sweeps are responsible for
sharpening the model towards the true one, hence we use the standard
regularization (5.2). The continuation sweeps with the smoothing
regularization start from the first frequency, and end at the fourth frequency
($i_{start}=1,i_{end}=4$). The next sweeps using (5.2) start from the fourth
frequency, up to the last one. The number of GN or ALM iterations differ for
each inversion and are noted for each experiment separately. In each GN
iteration we apply 5 projected and preconditioned CG iterations to
approximately solve the inner Newton problem. As preconditioner for these
iterations, we use the inverse of the Hessian of the smoothing regularization
terms. This way, the low number of CG iterations also play a role in
regularization, the inversion is not so sensitive to the choice of $\alpha$
(see [15] for more details on this technique). Once a direction is found, we
apply a standard Armijo linesearch. To solve for $\mathbf{Z}_{1}$ in either
(3.1) or (4.1), we also apply 5 (quadratic) CG iterations, at each update for
$\mathbf{Z}_{1}$ in an ALM iteration. Throughout all the relevant iterations,
the sparsity level of $\mathbf{Z}_{1}$ (percentage of non-zero values) was in
the range of $2\%-20\%$, contributing for the low memory footprint.

(a) True velocity model.

(b) The reference initial model.

Figure 1: The 2D Marmousi velocity model and initial guess in [km/s].

###  5.3 Marmousi model results

Our first set of experiments is conducted using the Marmousi model presented
in Fig. 1a. For this experiment we use a grid size of $550\times 200$
representing an area of size $9.192km\times 2.904km$. We place 136 sources and
549 receivers, which are uniformly spread on the top of the grid. The data is
generated by solving (2.1) for frequencies $\omega_{j}=2\pi f_{j}$, where

$\\{f_{j}\\}=\\{3,3.5,4.0,4.5,5.0,5.5,6.5,7.5,8.5\\}\;Hz,$

with the addition of 1% Gaussian i.i.d noise. The initial model for all the
inversions is given in Fig. 1b, and is initially used as $\mathbf{m}_{ref}$ in
the smoothing regularization.

For all the Marmousi experiments we run three sweeps of frequency continuation
- one sweep over the first 4 frequencies using the smoothing regularization
(5.1), and two additional sweeps over all the frequencies using the
regularization (5.2). All sweeps are obtained with 10 GN or ALM iterations for
each outer frequency continuation iteration. This results in total of 140 GN
iterations. Then, we applied up to 100 additional GN iterations involving the
four highest frequencies only, without source extensions. The additional
iterations are needed to sharpen the reconstructed model, especially for the
extended sources versions.

(a) Standard FWI.

(b) FWI with extended sources.

(c) FWI with extended and simultaneous sources.

Figure 2: Marmousi model reconstruction using the different FWI formulations.

####  5.3.1 Standard FWI

The reconstructed model for the experiment is shown at Fig. 2a, and the misfit
history plot is shown at Fig. 3a. This result converged to a local minima
where the misfit value (2.3) equals 6,957, and did not reconstruct the model
properly. We note that in all the convergence plots, the misfit values are
computed after the inversion is over for all the sources and frequencies
(independently of the frequency continuation schedule), and without the source
extensions.

####  5.3.2 FWI with extended sources

For the next experiment we added the extended sources (section 3) technique to
the first frequency continuation sweep with ALM (Alg. 2). The extended sources
has been used only in the first sweep since it helps with obtaining a good
initial smooth guess, and once we have such a guess the standard FWI achieves
good results. We set the rank of the source extensions to be $n_{es}=16$ (that
is the number of columns in $\mathbf{Z}_{1}$). In this experiment we chose the
initial $\beta$ parameters to be $\beta_{1}=0.1,\beta_{2}=10$. Fig. 2b shows
the reconstructed model and Fig. 3b shows the misfit history of the inversion.
The experiment results in good reconstruction of the model, which is closer to
the real model that the one obtained with standard FWI. The final misfit value
is 1,961, which is more than half the value of the standard FWI.

####  5.3.3 FWI with extended and simultaneous sources

To complete the experiments for the Marmousi model, we applied the joint
extended and simultaneous sources approach (Section 4), using the same
settings as in Sec. 5.3.2. We chose the random matrix $\mathbf{X}$ in (2.5) to
be of size $n_{s}\times p$, where $p=16$, so that in all the GN/ALM
iterations, the number of sources is 16. We chose this value through trial and
error, aiming at keeping $p$ small. The reconstructed model shown in Fig. 2c,
and is very similar to the result in Fig. 2b. This similarity is also evident
in the misfit history in Fig. 3c. This shows the addition of simultaneous
sources did not damage the effectiveness of the extended sources alone in
converging to a better minimum.

(a) Standard FWI.

(b) FWI with extended sources.

(c) FWI with extended and simultaneous sources.

Figure 3: Misfit values history, using all frequencies original sources only,
for the Marmousi model reconstruction.

####  5.3.4 Computational costs

To compare the computational costs of the experiment we count the amount of
forward simulations (2.1) that we solve in the first frequency continuation
sweep of each experiment. Since the forward simulations are the most
computationally expensive part of the algorithm, this gives a good comparison
between the different algorithms.

Standard FWI | FWI + SS | FWI + ES | FWI + ES + SS
---|---|---|---
164,560 | 19,360 | 230,296 | 48,416
Table 1: Comparison of computational costs for the Marmousi experiments.

The results in Table 1 demonstrate that even with the low-rank structure, the
addition of the extended sources (denoted as FWI+ES) increases the
computational costs quite significantly compared to standard FWI, and
especially compared to FWI with simultaneous sources (denoted FWI+SS). That is
partially because of the residual computation ($\mathbf{R}_{j}$ in Eq. (3.6)),
which is obtained for the full set of sources $\mathbf{Q}$. However, with the
addition of the simultaneous sources (denoted as FWI+ES+SS) the cost is
reduced significantly. The most computationally effective algorithm is FWI
with simultaneous sources only (denoted FWI+SS), which we do not demonstrate
here because its result is similar to standard FWI. Here we demonstrate that
our FWI+ES+SS version only roughly doubles the cost of FWI+SS, which, just
like FWI, does not involve an expanded search space and therefore is less
robust.

###  5.4 SEG/EAGE salt model results

(a) True SEG/EAGE salt model.

(b) The initial reference model.

Figure 4: The 2D SEG/EAGE salt velocity model and initial guess in [km/s].

Our second batch of experiments involves the SEG/EAGE salt model, presented in
Fig. 4a. The model is described by a $600\times 300$ grid, representing an
area of size $13.5km\times 4.2km$. We placed 119 sources and 599 receivers
uniformly spread at the top of the domain grid. The data are generated by
first solving (2.1) for frequencies $\omega_{j}=2\pi f_{j}$, where

$\\{f_{j}\\}=\\{3,3.3,3.6,3.9,4.2,4.5,5,5.5,6.5\\}\;Hz.$

We then add 1% Gaussian noise to the data. The initial model for all the
inversions is given in Fig. 4b, and is initially used as $\mathbf{m}_{ref}$ in
the smoothing regularization.

For the SEG/EAGE salt model we apply four frequency continuation sweeps - two
using the first 4 frequencies with the regularization (5.1), and two using all
frequencies with the regularization (5.2) to obtain a sharp model. For the
first two sweeps we used 20 GN iterations per outer iteration (with 7 CG
iterations in each inner iteration), and for the last two we used 15 GN
iterations (with 5 CG iterations in each, as in the rest of the
configurations). To finalize the inversion, we applied up to 100 additional GN
iterations involving the four highest frequencies, without source extensions
where relevant. These are used to sharpen the reconstructed model, especially
for the extended sources versions.

(a) Standard FWI.

(b) FWI with extended sources.

(c) FWI with extended and simultaneous sources.

Figure 5: SEG/EAGE salt model reconstruction using the different FWI
formulations.

####  5.4.1 Standard FWI

The reconstructed model for the experiment is shown in Fig. 5a, and the misfit
history plot is shown in Fig. 6a. The FWI seems to only reconstruct the upper
part of the salt block, and is missing the lower part. This is a typical
behavior of FWI that we wish to overcome. The final misfit value here was
2278, and the iterations stagnated.

####  5.4.2 FWI with extended sources

In this experiment we applied our inversion strategy with the extended sources
approach, which had been applied in the first continuation sweep only. In this
experiment we chose the initial $\beta$ parameters to be
$\beta_{1}=0.01,\beta_{2}=1$. The reconstructed model is shown at Fig. 5b, and
Fig. 6b shows the misfit history of the inversion. The addition of the
extended sources resulted in recovering the whole salt block, while keeping
the extended sources at low rank. The misfit value at the last iteration was
453 - much lower than with standard FWI.

(a) Standard FWI.

(b) FWI with extended sources.

(c) FWI with extended and simultaneous sources.

Figure 6: Misfit values history, using all frequencies and original sources,
for the SEG/EAGE salt model reconstruction.

####  5.4.3 FWI with extended sources and simultaneous sources

In the last experiment for the SEG/EAGE salt model, we use both extended and
simultaneous sources, with the same parameters as in the previous section. We
chose the $\mathbf{X}$ in (2.5) to be of size $n_{s}\times p$, with $p=16$. As
before, we see that the reconstructed model, shown in Fig. 5c, is very similar
to the result in Fig. 5b. Furthermore, the misfit plots for the current run
(6c) and previous run (6b) are similar as well, with final misfits of 591
versus 453, respectively. Like in the Marmousi experiment, the addition of the
simultaneous sources hardly affected the final reconstruction, and improved
the computational efficiency.

####  5.4.4 Computational costs

To compare the computational costs involved in the SEG/EAGE salt model
experiments for the different algorithms, we again count the amount of forward
simulations (2.1) we applied in the first frequency continuation sweep. The
results in Table 2 show the same trend as in the Marmousi case. The extended
sources increases the computational costs, while the addition of the
simultaneous sources reduces the cost drastically.

Standard FWI | FWI + SS | FWI + ES | FWI + ES + SS
---|---|---|---
285,600 | 38,400 | 418,102 | 108,640
Table 2: Computational costs comparison for the SEG/EAGE salt model
experiments.

###  5.5 SEG/EAGE Overthrust 3D model results

Our final experiment is applied to the central part of the 3D SEG/EAGE
Overthrust model, presented in Fig. 7a. The model is discretized on a grid of
size $172\times 172\times 108$, representing an area of size $7.5km\times
7.5km\times 4.65km$ (the center of the original model). We place 289 sources
and 1849 receivers which are spread uniformly in a 2D array on the top of the
grid. The data is generated by solving (2.1) for frequencies $\omega_{j}=2\pi
f_{j}$, where:

$\\{f_{j}\\}=\\{2.5,3.0,3.5,4.0,5.0\\}\;Hz,$

and adding white Gaussian noise of std $1\%$ of the magnitude of the data.
Starting from the initial reference model in Fig. 7b, we applied 3 cycles of
frequency continuation using our method with extended and simultaneous
sources. In this experiment we did not run the other configurations, as they
required too extensive time and careful management of memory swaps (between
RAM and the disk) for our resources. For the first cycle we used the smoothing
regularization with the lower 4 frequencies to obtain a smooth starting model
for the next cycles. The next two cycles were applied using the standard
regularization, starting from the 4-th frequency. Those cycles were applied
without extended sources (FWI + simultaneous sources alone), since the model
obtained after the first cycle was sufficient as a smooth guess. The resulting
reconstruction presented in Fig. 7c. We notice that the reconstructed model,
while not sharp enough, did manage to catch the important structures of the
true model. The reason for the rather smooth reconstruction is that the
highest frequency we used is only 5Hz to keep the model size reasonable—higher
frequencies would require larger grids, and a significantly more expensive
inversion.

(a) True velocity model.

(b) The reference initial model.

(c) Overthrust model reconstruction using FWI with extended and simultaneous
sources.

Figure 7: The 3D Overthrust velocity model and initial guess in [km/s].

##  6 Conclusion

In this work we aimed to improve on recent approaches for solving PDE-
constrained optimization problems—approaches that expand the search space in
order to relieve the non-linearity of the objective. In particular, we
considered the recent extended sources approach for FWI, and suggested a new
reduced version of this problem, where we couple the source extensions as a
low-rank matrix. We also showed that it is possible to accelerate the
minimization of our (source-extended) objective function by using simultaneous
sources, which reduces both memory and calculation costs. Unlike the previous
full-rank approach, ours does not require the additional computations or
memory for all the sources, which overrides the advantages of simultaneous
sources, and is prohibitively expensive in large scales. Therefore, our
approach is more applicable in real life 3D scenarios.

Our results showed that it is possible to combine the extended sources and the
simultaneous sources as we propose. On the one hand we use the source-
extension to achieve better reconstructions than the standard reduced FWI, and
on the other hand we are able to enjoy manageable computations and low-memory
footprint. We demonstrated our approach on two 2D models and one 3D model. The
latter, in particular, demonstrates the advantage of our approach—we were able
to apply the source-extended approach to a 3D problem using a rather standard
workstation.

Our method has two main limitations: one is the need to a-priory choose the
dimensions of both the low-rank source extension matrix and the dimension of
the simultaneous sources technique. Furthermore, we found that the dimension
of the latter has to be at least of the same size as the dimension of the
trace estimation, to prevent over-fitting and fluctuations when solving for
the source extensions.

## References

  * [1] H. S. Aghamiry, A. Gholami, and S. Operto, Implementing bound constraints and total-variation regularization in extended full-waveform inversion with the alternating direction method of multiplier: application to large contrast media, Geophysical Journal International, 218 (2019), pp. 855–872.
  * [2] , Improving full-waveform inversion by wavefield reconstruction with the alternating direction method of multipliers, Geophysics, 84 (2019), pp. R139–R162.
  * [3] , Robust wavefield inversion via phase retrieval, Geophysical Journal International, 221 (2020), pp. 1327–1340.
  * [4] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, and J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications, 23 (2001), pp. 15–41.
  * [5] F. Aminzadeh, B. Jean, and T. Kunz, 3-D salt and overthrust models, Society of Exploration Geophysicists, 1997.
  * [6] A. Y. Aravkin, D. Drusvyatskiy, and T. van Leeuwen, Efficient quadratic penalization through the partial minimization technique, IEEE Transactions on Automatic Control, 63 (2017), pp. 2131–2138.
  * [7] H. Avron and S. Toledo, Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix, Journal of the ACM (JACM), 58 (2011), p. 8.
  * [8] S. Bernard, V. Monteiller, D. Komatitsch, and P. Lasaygues, Ultrasonic computed tomography based on full-waveform inversion for bone quantitative imaging, Physics in Medicine & Biology, 62 (2017), p. 7011.
  * [9] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review, 59 (2017), pp. 65–98.
  * [10] B. Biondi and A. Almomin, Simultaneous inversion of full data bandwidth by tomographic full-waveform inversion, Geophysics, 79 (2014), pp. WA129–WA140.
  * [11] A. Brougois, M. Bourget, P. Lailly, M. Poulet, P. Ricarte, and R. Versteeg, Marmousi, model and data, in EAEG workshop-practical aspects of seismic data inversion, European Association of Geoscientists & Engineers, 1990, pp. cp–108.
  * [12] I. Epanomeritakis, V. Akcelik, O. Ghattas, and J. Bielak, A Newton-CG method for large-scale three-dimensional elastic full-waveform seismic inversion, Inverse Problems, (2008).
  * [13] Y. A. Erlangga, C. W. Oosterlee, and C. Vuik, A novel multigrid based preconditioner for heterogeneous Helmholtz problems, SIAM J. Sci. Comput., 27 (2006), pp. 1471–1492.
  * [14] L. Guasch, O. C. Agudo, M.-X. Tang, P. Nachev, and M. Warner, Full-waveform inversion imaging of the human brain, NPJ digital medicine, 3 (2020), pp. 1–12.
  * [15] E. Haber, Computational Methods in Geophysical Electromagnetics, vol. 1, SIAM, 2014.
  * [16] E. Haber, U. Ascher, and D. Oldenburg, On optimization techniques for solving nonlinear inverse problems, Inverse problems, 16 (2000), pp. 1263–1280.
  * [17] E. Haber and U. M. Ascher, Preconditioned all-at-once methods for large, sparse parameter estimation problems, Inverse Problems, 17 (2001), p. 1847.
  * [18] E. Haber, M. Chung, and F. Herrmann, An effective method for parameter estimation with PDE constraints with multiple right-hand sides, SIAM Journal on Optimization, 22 (2012), pp. 739–757.
  * [19] E. Haber, E. Treister, and E. Holtham, Obtaining low frequencies for full waveform inversion by using augmented physics, ASEG Extended Abstracts, (2016), pp. 1–5.
  * [20] G. Huang, R. Nammour, and W. W. Symes, Volume source based extended waveform inversion, Geophysics, 83 (2018), pp. 1–139.
  * [21] J. R. Krebs, J. E. Anderson, D. Hinkley, R. Neelamani, S. Lee, A. Baumstein, and M.-D. Lacasse, Fast full-wavefield seismic inversion using encoded sources, Geophysics, 74 (2009), pp. WCC177–WCC188.
  * [22] Y. E. Li and L. Demanet, Full-waveform inversion with extrapolated low-frequency data, Geophysics, 81 (2016), pp. R339–R348.
  * [23] M. Liu, R. Kumar, E. Haber, and A. Aravkin, Simultaneous-shot inversion for PDE-constrained optimization problems with missing data, Inverse Problems, 35 (2018), p. 025003.
  * [24] Z. Liu and J. Zhang, Joint traveltime, waveform, and waveform envelope inversion for near-surface imaging, Geophysics, 82 (2017), pp. R235–R244.
  * [25] L. Métivier, R. Brossier, S. Operto, and J. Virieux, Full waveform inversion and the truncated newton method, SIAM Review, 59 (2017), pp. 153–195.
  * [26] J. Poulson, B. Engquist, S. Li, and L. Ying, A parallel sweeping preconditioner for heterogeneous 3D Helmholtz equations, SIAM J. Sci. Comput., 35 (2013), pp. C194–C212.
  * [27] R. Pratt, Seismic waveform inversion in the frequency domain, part 1: Theory, and verification in a physical scale model, Geophysics, 64 (1999), pp. 888–901.
  * [28] R. G. Pratt, C. Shin, and G. Hick, Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion, Geophys. J. Int., 133 (1998), pp. 341–362.
  * [29] Y. Rao and Y. Wang, Seismic waveform tomography with shot-encoding using a restarted l-bfgs algorithm, Scientific reports, 7 (2017), pp. 1–9.
  * [30] F. Roosta-Khorasani and U. Ascher, Improved bounds on sample size for implicit matrix trace estimators, Foundations of Computational Mathematics, 15 (2015), pp. 1187–1212.
  * [31] F. Roosta-Khorasani, K. Van Den Doel, and U. Ascher, Data completion and stochastic algorithms for PDE inversion problems with many measurements, Electron. Trans. Numer. Anal, 42 (2014), pp. 177–196.
  * [32] F. Roosta-Khorasani, K. van den Doel, and U. Ascher, Stochastic algorithms for inverse problems involving PDEs and many measurements, SIAM J. Sci. Comput., 36 (2014), pp. S3–S22.
  * [33] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena, 60 (1992), pp. 259–268.
  * [34] L. Ruthotto, E. Treister, and E. Haber, jInv – a flexible Julia package for PDE parameter estimation, SIAM J. Sci. Comput., 39 (2017), pp. S702––S722.
  * [35] F. Santosa and W. W. Symes, Linear inversion of band-limited reflection seismograms, SIAM Journal on Scientific and Statistical Computing, 7 (1986), pp. 1307–1330.
  * [36] O. Schenk and K. Gärtner, Solving unsymmetric sparse systems of linear equations with pardiso, Future Generation Computer Systems, 20 (2004), pp. 475–487.
  * [37] C. Shin and Y. Ho Cha, Waveform inversion in the Laplace—Fourier domain, Geophysical Journal International, 177 (2009), pp. 1067–1079.
  * [38] E. Somersalo and J. Kaipio, Statistical and computational inverse problems, Applied Mathematical Sciences, 160 (2004).
  * [39] E. Soubies, T.-A. Pham, and M. Unser, Efficient inversion of multiple-scattering model for optical diffraction tomography, Optics express, 25 (2017), pp. 21786–21800.
  * [40] A. Tarantola, Inverse problem theory, Elsevier, Amsterdam, 1987.
  * [41] E. Treister and E. Haber, Full waveform inversion guided by travel time tomography, SIAM J. Sci. Comput., 39 (2017), pp. S587––S609.
  * [42] E. Treister and I. Yavneh, A multilevel iterated-shrinkage approach to $l_{1}$ penalized least-squares minimization, IEEE Transactions on Signal Processing, 60 (2012), pp. 6319–6329.
  * [43] T. van Leeuwen, A. Y. Aravkin, and F. J. Herrmann, Seismic waveform inversion by stochastic optimization, International Journal of Geophysics, 2011 (2011).
  * [44] T. van Leeuwen and F. J. Herrmann, Mitigating local minima in full-waveform inversion by expanding the search space, Geophysical Journal International, 195 (2013), pp. 661–667.
  * [45] , 3D frequency-domain seismic inversion with controlled sloppiness, SIAM J. Sci. Comput., 36 (2014), pp. S192–S217.
  * [46] T. van Leeuwen and F. J. Herrmann, A penalty method for PDE-constrained optimization in inverse problems, Inverse Problems, 32 (2016), p. 015007.
  * [47] C. R. Vogel, Computational methods for inverse problems, vol. 23, SIAM, Philadelphia, 2002.
  * [48] M. Zibulevsky and M. Elad, L1-l2 optimization in signal and image processing, IEEE Signal Processing Magazine, 27 (2010), pp. 76–88.