# Randomized Source Sketching for Full Waveform Inversion

Stéphane Operto, Kamal Aghazade, Hossein Aghamiry, Ali Gholami

► To cite this version:

Stéphane Operto, Kamal Aghazade, Hossein Aghamiry, Ali Gholami. Randomized Source Sketching for Full Waveform Inversion. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60, pp.1-11. 10.1109/TGRS.2021.3131039. [hal-03478761]

HAL Id: hal-03478761

https://hal.science/hal-03478761v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Randomized Source Sketching for Full Waveform Inversion

Kamal Aghazade, Hossein S. Aghamiry<sup>©</sup>, Member, IEEE, Ali Gholami<sup>†</sup>, and Stéphane Operto

**Abstract**—Partial differential equation (PDE)-constrained optimization problems such as seismic full waveform inversion (FWI) frequently arise in the geoscience and related fields. For such problems, many observations are usually gathered by multiple sources, which form the right-hand sides of the PDE constraint. Solving the inverse problem with such massive datasets is computationally demanding, in particular, when dealing with a large number of model parameters. This article proposes a novel randomized source sketching method for the efficient resolution of multisource PDE-constrained optimization problems. We first formulate the different source-encoding strategies used in seismic imaging into a unified framework based on a randomized sketching. To this end, the source dimension of the problem is projected in a smaller domain by a suitably defined projection matrix that gathers the physical sources in super-sources through a weighted summation. This reduction in the number of physical sources decreases significantly the number of PDE solves while suitable sparsity-promoting regularization can efficiently mitigate the footprint of the crosstalk noise to maintain the convergence speed of the algorithm. We implement the randomized sketching method in an extended search-space formulation of frequency-domain FWI, which relies on the alternating-direction method of multipliers (ADMMs). Numerical examples carried out with a series of well-documented 2-D benchmarks demonstrate that the randomized sketching algorithm reduces the cost of large-scale problems by at least one order of magnitude compared to the original deterministic algorithm.

**Index Terms**—Alternating direction methods of multipliers (ADMMs), partial differential equation (PDE), randomized source sketching, seismic full waveform inversion (FWI), sparsity-promoting regularization.

## I. INTRODUCTION

PARTIAL differential equation (PDE)-constrained optimization problems frequently arise in the geoscience and related fields whenever the relation between the parameters of interest and the observations (forward problem) is described by a PDE such as the wave equation in seismology [1], Poisson-type equations in dc resistivity [2], and Maxwell's equations in ground-penetrating radar (GPR) [3].

This article deals with full waveform inversion (FWI), which is governed by the wave equation. FWI is an attractive inversion scheme with resolution capabilities at the wavelength

scale to build different subsurface parameters from sparse measurements of the propagated wavefields [1], [4]–[6]. The inverse problem is usually formulated as the minimization of a misfit distance between the calculated and observed data subject to a set of equality constraints [7].

Traditional FWI methods (also known as reduced approaches) strictly enforce the wave-equation constraint to be satisfied during iterations by eliminating the wavefield, leading to a monovariate unconstrained optimization problem for model parameters only. It is usually solved iteratively using gradient-based numerical methods in the time or frequency domain [6]. This study is focused on the frequency-domain formulation of FWI, where seismic modeling requires the solution of large and sparse indefinite linear systems with multiple right-hand sides. These local optimizations methods, however, are prone to cycle skipping, which drives the inversion toward spurious local minima when the initial model does not allow one to match the recorded data with an error smaller than half a period [8].

To manage the cycle-skipping issue, a new generation of algorithms known as the extended waveform inversion append extra degrees of freedom to the problem, which bring flexibility for fitting the data with inaccurate velocity models and hence reduce the cycle skipping pathology [7], [9]–[11]. In such extended-space approaches, wavefields or sources are processed as optimization variables in addition to the model parameters, which implies that the PDE is solved approximately (i.e., with a relaxation) unlike the case of the reduced approach. Among them, the wavefield reconstruction inversion (WRI) [10] and iteratively refined WRI (IR-WRI) [11] showed promising results for accurate parameter estimation from rough starting models. The WRI method is based on a penalty formulation of the constrained problem, which requires a fine-tuning of the penalty parameter and has a slow convergence rate. IR-WRI overcomes these drawbacks by using the augmented Lagrangian (AL) function, which has several advantages over the penalty counterpart [12]. This study is thus focused on increasing the computational efficiency of IR-WRI based on randomized sketching methods, usually referred to as simultaneous sources or source encoding strategies.

Generally, waveform inversion is computationally intensive, which makes 3-D field-scale case studies challenging. The main computational burden results from the solution of the PDE with multiple right-hand sides. In the frequency domain, this PDE is typically solved with direct or iterative methods [13], [14]. In the latter case which may be more suitable for very large problems (more than 100 million degrees of freedom), the computational complexity theoretically scales to

Manuscript received August 6, 2021; revised September 27, 2021; accepted November 23, 2021. Date of publication November 30, 2021; date of current version February 25, 2022. This work was supported in part by the Waveform Inversion of Node Data (WIND) Consortium (https://www.geoazur.fr/WIND), sponsored by Chevron, Shell, and Total. (Corresponding author: Ali Gholami.)

Kamal Aghazade and Ali Gholami are with the Institute of Geophysics, University of Tehran, Tehran 1417466191, Iran (e-mail: aghazade.kamal@ut.ac.ir; agholami@ut.ac.ir).

Hossein S. Aghamiry and Stéphane Operto are with Geoazur, CNRS-IRD-OCA, University Côte d'Azur, 06100 Valbonne, France (e-mail: aghamiry@geoazur.unice.fr; operto@geoazur.unice.fr).

Digital Object Identifier 10.1109/TGRS.2021.3131039
the number of right-hand sides (source terms). Even if block
iterative methods and recycling can speed up the processing
of multiple right-hand sides [15], [16], it remains beneficial to
mitigate the computational burden of multisource simulations
for dense acquisitions. Therefore, along with the developments
of efficient solvers, some cost-effective strategies based on
dimensionality reduction methods have been designed to han-
dle computational challenges of waveform inversion both in
terms of memory requirement and speed up purposes. Data
decimation (see [5], [17]), data compression, and determin-
istic/random stacking methods [18] are preliminary works
toward this goal. Decimation methods may be unfavorable in
areas with inadequate illumination [19]. Deterministic stacking
of sources, as well as data compression methods, suffer from
crosstalk noise artifacts due to the unbalanced summation
of the individual sources [20]. As a deterministic encoding
strategy, the plane-wave encoding reduces the dimensionality
of data by inverting the dominant components of the data
after transformation into the plane-wave domain (linear Radon
transform) [21]–[25]. However, the linear Radon transform is
not able to optimally focus the seismic events in the trans-
form domain. This makes the performance of the plane-wave
encoding highly dependent on the number of ray parameters
within an appropriate range [24], [26].

Random projection methods consider a linear combination
of the sources by applying some random codes and are
known to be more effective compared with other deterministic
methods [19], [27]–[30].

Reducing the number of physical sources to a few
super-sources decreases the computational cost dramati-
cally [19], [28], [30], [31], but still two main challenges
remain, namely the crosstalk noise induced by the source
blending and the slow convergence rate of the algorithm
when significant reduction of the source dimensionality is
performed. In what regards the first issue, one may investigate
different random codes [30], [32], [33], regenerate random
codes at each iteration [33], use random source positions [34],
or employ a regularization method [35]–[37]. Regarding the
convergence issue, one may employ second-order optimization
algorithms to accelerate the convergence rate [19], [38], [39].
Castellanos et al. [19] show that for non-convex problems,
random encoding may explore areas in the model space that
is not accessible to the deterministic counterpart. Ben-Hadj-
Ali et al. [30] show that implementation of multi-frequency
group inversion and partial assembling of the sources make
the inversion more robust against both crosstalk and phys-
ical noises. Finally, [40] concluded that inverting partially
overlapped frequency groups yields high-quality results of the
source encoding method.

The driving motivation of this article is to increase the
computational efficiency of the recent IR-WRI algorithm
with randomized sketching strategies [41], [42]. Randomized
IR-WRI first requires to define the associated AL function
in an appropriate transform domain of lower dimension. For
this, the source dimension of the problem is projected by
a sketching matrix in the transform domain, whereby the
physical sources are replaced by super-sources, each can be
a combination of the physical sources. For a satisfactory

recovery performance of the lower-dimensional problem, the
Gram matrix's expectation associated with the sketching
matrix approaches identity. In this case, the projected AL
function converges to the original function. Two categories
of sketching matrices are proposed. The first category includes
unstructured random matrices with independent and identically
distributed (i.i.d.) entries, including those from the Gaussian
and Bernoulli distributions. The second category includes
structured random matrices that are built by randomizing
orthogonal matrices. As special cases, we use matrices asso-
ciated with the Fourier transform, discrete cosine transform
(DCT), discrete wavelet transform (DWT), noiselet trans-
form [43], and Walsh-Hadamard transform [44].

There are three key contributions of this article as follows.

1) We develop a general algorithm for randomization of extended waveform inversion in the frame of the AL method.

2) This randomized scheme admits the existing encoding methods as special cases.

3) We demonstrate with numerical examples that the ran-
domized algorithm provides improved performance over
existing deterministic algorithms.

II. PRELIMINARIES

The basic definitions and symbols used in this article are
collected in this section and are as follows. Real, complex,
integers, and natural numbers are denoted by $\mathbb{R}$, $\mathbb{C}$, $\mathbb{Z}$, and $\mathbb{N}$,
respectively. The imaginary unit is denoted by $i = (-1)^{1/2}$.
The number of discrete model parameters, the number of
receivers, the number of sources, and the number of sketched
sources are denoted by $n$, $m$, $p$, and $q$, respectively. The
angular frequency is denoted by $\omega = 2\pi f$. The identity
matrix and diagonal matrices are denoted by $\mathbf{I}$ and $\text{Diag}(\bullet)$,
respectively, and $\nabla^2$ is the Laplace operator. For a matrix
$\mathbf{X} = [\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_n]$, $\mathbf{x}_j$ denotes the $j$th column of $\mathbf{X}$ and $\mathbf{x}(i)$
represents the $i$th entry of the column vector $\mathbf{x}$. The conjugate
transpose of a matrix/vector is denoted by superscript $T$.
The Frobenius norm of a matrix $\mathbf{X}$ is defined as $\|\mathbf{X}\|_F = (\sum_{i,j} |X(i,j)|^2)^{1/2}$ and for two matrices $\mathbf{X}$ and $\mathbf{Y}$ of the
same dimension the inner product is defined as $\langle \mathbf{X}, \mathbf{Y} \rangle = \text{tr}(\mathbf{X}^T \mathbf{Y}) = \sum_{i,j} \bar{X}(i,j) Y(i,j)$, where the bar shows complex
conjugate and $(i,j)$ denote row-and-column index of matrix.
$(\bullet)^k$ represents the value of $(\bullet)$ at iteration $k$.

We also define mode error, ME = $100 \times \|\mathbf{m}^* - \mathbf{m}\|_2 / \|\mathbf{m}^*\|_2$ as a convergence measure of the methods during iterations, wherein $\mathbf{m}^*$ and $\mathbf{m}$ are true and estimated models, respectively.

III. METHOD

In this section, we develop a new algorithm of IR-WRI [11]
based upon a randomized sketching scheme. We formulate
the AL method in the transform domain of lower dimension
and review a quite general formulation that admits existing
source-encoding strategies as special cases.

A. Extended FWI With AL Method

Multisource FWI problem with a general regularization
term $\mathcal{R}(\mathbf{m}) : \mathbb{R}^{n \times 1} \rightarrow \mathbb{R}$ can be expressed as the following
equality-constrained bi-convex optimization problem [7], [45]:

$$ 
\begin{array}{l}
\min_{\mathbf{U},\mathbf{m} \in \mathcal{B}} \mathcal{R}(\mathbf{m}) \\
\text{s.t. } \mathbf{A}(\mathbf{m})\mathbf{U} = \mathbf{B} \\
\qquad \mathbf{P}\mathbf{U} = \mathbf{D}
\end{array}
\qquad (1) $$

where the columns of $\mathbf{U} \in \mathbb{C}^{n \times p}$ gather $p$ wavefields represented by vectors of dimension $n$. The columns of $\mathbf{D} \in \mathbb{C}^{m \times p}$ contain the associated $p$ observation vectors of dimension $m$. The matrix $\mathbf{P} \in \mathbb{R}^{m \times n}$ is the observation operator that samples the wavefields at the receiver locations. Through this article, we assume a fixed-spread acquisition where $\mathbf{P}$ is independent of the sources. The matrix $\mathbf{A}(\mathbf{m}) \in \mathbb{C}^{n \times n}$ is the discretized Helmholtz operator with appropriate boundary conditions [46] and the columns of $\mathbf{B} \in \mathbb{C}^{n \times p}$ gather the source terms. Finally, the vector $\mathbf{m} \in \mathbb{R}^{n \times 1}$ gathers the model parameters (squared slowness), and $\mathcal{B}$ is a convex set defined by known lower and upper bounds of the model. In this article, we develop our algorithm using the frequency-domain formulation of FWI, but the resulting algorithm can also be applied in the time domain.

The traditional FWI solves (1) by eliminating the wavefields, that is, $\mathbf{U} = \mathbf{A}(\mathbf{m})^{-1}\mathbf{B}$, leading to the following unconstrained monovariate optimization problem:

$$ \min_{\mathbf{m} \in \mathcal{B}} \mathcal{R}(\mathbf{m}) + \frac{\mu}{2} \|\mathbf{P}\mathbf{A}(\mathbf{m})^{-1}\mathbf{B} - \mathbf{D}\|_{\mathrm{F}}^{2} \quad (2) $$

where $\mu > 0$ is the penalty term. Gradient-based optimization methods are usually used to minimize this objective function where the gradient of the misfit function is efficiently computed with the adjoint-state method [47]. However, there is a general agreement on the sensitivity of FWI to the accuracy of the starting model when gradient-based optimization methods are used to solve it. The AL method was proposed to mitigate this sensitivity.

The AL function of the constrained problem (1) is given by [11], [48]

$$ 
\begin{aligned}
\mathcal{L}(\mathbf{U}, \mathbf{m}, \mathbf{\Lambda}, \mathbf{\Gamma}) = & \mathcal{R}(\mathbf{m}) + \frac{\alpha}{2} \|\mathbf{A}(\mathbf{m})\mathbf{U} - \mathbf{B}\|_{\mathrm{F}}^{2} \\
& + \frac{\beta}{2} \|\mathbf{P}\mathbf{U} - \mathbf{D}\|_{\mathrm{F}}^{2} - \langle \mathbf{\Lambda}, \mathbf{A}(\mathbf{m})\mathbf{U} - \mathbf{B} \rangle \\
& - \langle \mathbf{\Gamma}, \mathbf{P}\mathbf{U} - \mathbf{D} \rangle
\end{aligned}
\qquad (3) $$

where $\alpha, \beta > 0$ are penalty parameters. The matrices $\mathbf{\Lambda} \in \mathbb{C}^{n \times p}$ and $\mathbf{\Gamma} \in \mathbb{C}^{m \times p}$ gather the Lagrange multipliers associated with the PDE constraints and the observation constraints, respectively, each column of these matrices being tied to a particular source. The AL function (3) recasts FWI as the following min-max problem:

$$ \min_{\mathbf{U},\mathbf{m}} \max_{\mathbf{\Lambda},\mathbf{\Gamma}} \mathcal{L}(\mathbf{U}, \mathbf{m}, \mathbf{\Lambda}, \mathbf{\Gamma}). \qquad (4) $$

IR-WRI solves (4) iteratively with the alternating direction method of multipliers (ADMM) [49], where each primal and dual variables are updated in alternating mode as

$$ \mathbf{U}_{k+1} = \arg \min_{\mathbf{U}} \mathcal{L}(\mathbf{U}, \mathbf{m}_k, \mathbf{\Lambda}_k, \mathbf{\Gamma}_k) \qquad (5a) $$

$$ \mathbf{m}_{k+1} = \arg \min_{\mathbf{m}} \mathcal{L}(\mathbf{U}_{k+1}, \mathbf{m}, \mathbf{\Lambda}_k, \mathbf{\Gamma}_k) \qquad (5b) $$

$$ \mathbf{\Lambda}_{k+1} = \mathbf{\Lambda}_k - \alpha(\mathbf{A}(\mathbf{m}_{k+1})\mathbf{U}_{k+1} - \mathbf{B}) \qquad (5c) $$

$$ \mathbf{\Gamma}_{k+1} = \mathbf{\Gamma}_k - \beta(\mathbf{P}\mathbf{U}_{k+1} - \mathbf{D}) \qquad (5d) $$

beginning with $\mathbf{\Lambda}_0 = \mathbf{0}$ and $\mathbf{\Gamma}_0 = \mathbf{0}$, and $k$ stands for the iteration number. The dual updating in (5c) and (5d) can be interpreted as gradient ascent steps with respect to $\mathbf{\Lambda}$ and $\mathbf{\Gamma}$ (with step lengths $\alpha$ and $\beta$) applied to the dual function or to the AL function while the primal variables are kept fixed [12]. We now provide the solution of the wavefield and model subproblems.

1) **Wavefield Subproblem:** Solving the wavefield subproblem (5a) by differentiating the AL function with respect to $\mathbf{U}$ and equating the result to zero gives the closed-form expression of $\mathbf{U}_{k+1}$ as the solution of the following linear system with multiple right-hand sides:

$$ (\alpha \mathbf{A}_k^\mathrm{T} \mathbf{A}_k + \beta \mathbf{P}^\mathrm{T} \mathbf{P}) \mathbf{U}_{k+1} = \alpha \mathbf{A}_k^\mathrm{T} \left( \mathbf{B} - \frac{1}{\alpha} \mathbf{\Lambda}_k \right) + \beta \mathbf{P}^\mathrm{T} \left( \mathbf{D} - \frac{1}{\beta} \mathbf{\Gamma}_k \right). \qquad (6) $$

The computational burden of solving this system is proportional to the number of right-hand sides when iterative solvers, which are the only choice for large-scale problems, are used. In Section IV, we use randomized sketching methods to reduce the dimensionality of the problem and decrease the computational burden of the inverse problem.

2) **Model Subproblem:** The objective function for the model parameters in (5b) reduces to

$$ \min_{\mathbf{m} \in \mathcal{B}} \mathcal{R}(\mathbf{m}) + \frac{\alpha}{2} \left\| \mathbf{A}(\mathbf{m})\mathbf{U}_{k+1} - \mathbf{B} + \frac{1}{\alpha} \mathbf{\Lambda}_k \right\|_{\mathrm{F}}^2. \qquad (7) $$

Compared to the classical FWI, (2), in which the nonquadratic misfit function is defined in the data space, the misfit function is defined in the source space in (7) and is quadratic with respect to the model parameters. Considering the special structure of the Helmholtz operator

$$ \mathbf{A}(\mathbf{m}) = \nabla^2 + \omega^2 \mathrm{Diag}(\mathbf{m}) \qquad (8) $$

we can rewrite (7) in the following equivalent form:

$$ \min_{\mathbf{m} \in \mathcal{B}} \mathcal{R}(\mathbf{m}) + \frac{\alpha \omega^2}{2} \sum_{j=1}^{p} \| \mathrm{Diag}(\mathbf{u}_j^{k+1}) \mathbf{m} - \mathbf{y}_j^k \|_{2}^{2} \qquad (9) $$

where $\mathbf{y}_j^k = (1/\omega^2)(\mathbf{b}_j - (1/\alpha)\mathbf{\lambda}_j^k - \nabla^2 \mathbf{u}_j^{k+1})$, and $\mathbf{b}_j, \mathbf{\lambda}_j^k$ and $\mathbf{u}_j^{k+1}$ are the jth column of matrices $\mathbf{B}, \mathbf{\Lambda}_k$, and $\mathbf{U}_{k+1}$, respectively. The optimization problem (9) can be solved efficiently with different regularization function, which may force smoothness, sparseness, and blockiness. Alternatively, a combination of them can be used to compute the regularized solution [50]. Applications of such regularizers in IR-WRI are reviewed in [45] and [51]. In this study, this regularization term is essential to mitigate the footprint of data over-fitting resulting from crosstalk noise.

## IV. IR-WRI WITH SOURCE SKETCHING

For two matrices $\mathbf{X}, \mathbf{Y} \in \mathbb{C}^{m \times p}$ and sketching matrix $\mathbf{S} \in \mathbb{C}^{p \times q}$ such that $\mathbf{S}\mathbf{S}^\mathrm{T} = \mathbf{I}$ (the identity matrix) we have

$$ \langle \mathbf{X}, \mathbf{Y} \rangle = \langle \mathbf{X}, \mathbf{Y}\mathbf{I} \rangle \qquad (10a) $$

$$ = \langle \mathbf{X}, \mathbf{Y}\mathbf{S}\mathbf{S}^\mathrm{T} \rangle \qquad (10b) $$

$$ = \langle \mathbf{X}\mathbf{S}, \mathbf{Y}\mathbf{S} \rangle \qquad (10c) $$

$$ = \langle \mathbf{\hat{X}}, \mathbf{\hat{Y}} \rangle \qquad (10d) $$
where the hat shows row-transformed matrix obtained by right multiplication by the sketching matrix $\mathbf{S}$. Furthermore, the AL function in (3) can be written as

$$ \mathcal{L}(\mathbf{U}, \mathbf{m}, \mathbf{\Lambda}, \mathbf{\Gamma}) = \mathcal{R}(\mathbf{m}) + \frac{\alpha}{2} \langle \delta\mathbf{B}, \delta\mathbf{B} \rangle + \frac{\beta}{2} \langle \delta\mathbf{D}, \delta\mathbf{D} \rangle - \langle \mathbf{\Lambda}, \delta\mathbf{B} \rangle - \langle \mathbf{\Gamma}, \delta\mathbf{D} \rangle \quad (11) $$

where $\delta\mathbf{D} = \mathbf{PU} - \mathbf{D}$ and $\delta\mathbf{B} = \mathbf{A}(\mathbf{m})\mathbf{U} - \mathbf{B}$. Using (10c), the AL function can be defined in the transform domain as

$$ \hat{\mathcal{L}}(\hat{\mathbf{U}}, \mathbf{m}, \hat{\mathbf{\Lambda}}, \hat{\mathbf{\Gamma}}) = \mathcal{R}(\mathbf{m}) + \frac{\alpha}{2} \|\mathbf{A}(\mathbf{m})\hat{\mathbf{U}} - \hat{\mathbf{B}}\|_{\mathrm{F}}^{2} \\ + \frac{\beta}{2} \|\mathbf{P}\hat{\mathbf{U}} - \hat{\mathbf{D}}\|_{\mathrm{F}}^{2} - \langle \hat{\mathbf{\Lambda}}, \mathbf{A}(\mathbf{m})\hat{\mathbf{U}} - \hat{\mathbf{B}} \rangle \\ - \langle \hat{\mathbf{\Gamma}}, \mathbf{P}\hat{\mathbf{U}} - \hat{\mathbf{D}} \rangle \quad (12) $$

where $\hat{\mathbf{B}} = \mathbf{BS}$, $\hat{\mathbf{D}} = \mathbf{DS}$, $\hat{\mathbf{U}} = \mathbf{US}$, $\hat{\mathbf{\Lambda}} = \mathbf{\Lambda}\mathbf{S}$, and $\hat{\mathbf{\Gamma}} = \mathbf{\Gamma}\mathbf{S}$, and the IR-WRI algorithm in (5) can be readily applied to the AL function in the transform domain. For $\mathbf{B} \in \mathbb{C}^{n \times p}$ and $\mathbf{S} \in \mathbb{R}^{p \times q}$, $\hat{\mathbf{B}} \in \mathbb{C}^{n \times q}$ thus for $q \ll p$ we expect significant computational saving if we perform optimization over $\mathcal{L}$ rather than $\mathcal{L}$ while at the same time being able to maintain the good quality of the estimates. This indeed results because, if $q \ll p$, the sketched PDE (6) is of a much smaller dimension than the original one and hence easier to solve. However, for $q < p$, the equality $\mathbf{SS}^{\mathrm{T}} = \mathbf{I}$ is no longer be satisfied. Thus, we deal with cases where $\mathbf{SS}^{\mathrm{T}} \approx \mathbf{I}$, which requires us to work with an approximate AL function $\hat{\mathcal{L}}$. For a fixed matrix $\mathbf{S}$, one can build the projected source/data matrices $\hat{\mathbf{B}}/\hat{\mathbf{D}}$ once and for all before the inversion and then apply IR-WRI to the projected source and data. In this case, however, the imprint of crosstalk noises will propagate into the model update unless $q$ is sufficiently large such that $\mathbf{SS}^{\mathrm{T}} \approx \mathbf{I}$. To address the crosstalk noise effects while maintaining the computational efficiency of the algorithm at the desired level, we apply a randomization scheme. That is, we regenerate the projection matrix $\mathbf{S}$ at each IR-WRI iteration. In this case, the physical constraint errors, i.e., $\mathbf{PU} - \mathbf{D}$ and $\mathbf{AU} - \mathbf{B}$, may be removed more efficiently in iterations by keeping the dual variables in the original space. To do this, we rewrite the transformed AL function (12) in an equivalent form by using the equality (10b) for the Lagrangian terms

$$ \hat{\mathcal{L}}(\hat{\mathbf{U}}, \mathbf{m}, \mathbf{\Lambda}, \mathbf{\Gamma}) = \mathcal{R}(\mathbf{m}) + \frac{\alpha}{2} \|\mathbf{A}(\mathbf{m})\hat{\mathbf{U}} - \hat{\mathbf{B}}\|_{\mathrm{F}}^{2} \\ + \frac{\beta}{2} \|\mathbf{P}\hat{\mathbf{U}} - \hat{\mathbf{D}}\|_{\mathrm{F}}^{2} - \langle \mathbf{\Lambda}, (\mathbf{A}(\mathbf{m})\hat{\mathbf{U}} - \hat{\mathbf{B}})\mathbf{S}^{\mathrm{T}} \rangle \\ - \langle \mathbf{\Gamma}, (\mathbf{P}\hat{\mathbf{U}} - \hat{\mathbf{D}})\mathbf{S}^{\mathrm{T}} \rangle. \quad (13) $$

Although (12) and (13) are formally equivalent by virtue of equality (10b), they don't lead to the same convergence path in iterations when $\mathbf{S}$ is regenerated at each iteration. Unlike (12), (13) calculates the constraint errors in the full space rather than in the sketching domain at each iteration $k$ by “deblending” the sketched constraint errors with the $\mathbf{S}^{\mathrm{T}}$, i.e., an approximate inverse of the current $\mathbf{S}$.

The randomized IR-WRI is summarized in Algorithm 1. Generally, the sketching matrix $\mathbf{S}$ and the scalar parameter $q$ can be seen as additional parameters of the randomized IR-WRI compared to the traditional IR-WRI. In Sections IV-A and IV-B, we provide possible choices of

Algorithm 1 IR-WRI With Randomized Source Sketching

Require: $\mathbf{D}, \mathbf{B}$

1: set $\mathbf{\Lambda}_0 = \mathbf{0}$, $\mathbf{\Gamma}_0 = \mathbf{0}$

2: repeat

3: Sample an independent copy $\mathbf{S}$

4: Set $\hat{\mathbf{B}} = \mathbf{BS}$, $\hat{\mathbf{D}} = \mathbf{DS}$, $\hat{\mathbf{\Lambda}}_k = \mathbf{\Lambda}_k\mathbf{S}$, $\hat{\mathbf{\Gamma}}_k = \mathbf{\Gamma}_k\mathbf{S}$

5: $\hat{\mathbf{U}}_{k+1} = \mathbf{H}^{-1}[\alpha\mathbf{A}_k^T(\hat{\mathbf{B}} - \frac{1}{\alpha}\hat{\mathbf{\Lambda}}_k) + \beta\mathbf{P}^T(\hat{\mathbf{D}} - \frac{1}{\beta}\hat{\mathbf{\Gamma}}_k)]$

6: $\mathbf{m}_{k+1} = \operatorname*{argmin}_{\mathbf{m}} \mathcal{R}(\mathbf{m}) + \frac{\alpha}{2}\|\mathbf{A}(\mathbf{m})\hat{\mathbf{U}}_{k+1} - \hat{\mathbf{B}} + \frac{1}{\alpha}\hat{\mathbf{\Lambda}}_k\|_{\mathrm{F}}^2$

7: $\mathbf{\Lambda}_{k+1} = \mathbf{\Lambda}_k - \alpha(\mathbf{A}(\mathbf{m}_{k+1})\hat{\mathbf{U}}_{k+1} - \hat{\mathbf{B}})\mathbf{S}^{\mathrm{T}}$

8: $\mathbf{\Gamma}_{k+1} = \mathbf{\Gamma}_k - \beta(\mathbf{P}\hat{\mathbf{U}}_{k+1} - \hat{\mathbf{D}})\mathbf{S}^{\mathrm{T}}$

9: $k = k + 1$.

10: until stopping conditions are satisfied.

the sketching matrix $\mathbf{S}$, which are built either randomly or from randomized orthonormal matrices.

A. Random Sketching Matrices

We would like the random matrix $\mathbf{S}$ to satisfy

$$ \mathbb{E}(\mathbf{S}\mathbf{S}^{\mathrm{T}}) = \mathbf{I} \quad (14) $$

where $\mathbb{E}$ denotes expectation over random matrix $\mathbf{S}$ with a pre-defined distribution. A truncated expectation has been widely used for stochastic estimation of the trace of a matrix [52]. The stochastic trace estimation has also been applied to randomized/stochastic FWI problems [53], [54]. In the following, we review different forms of $\mathbf{S}$, which are referred to as random projection matrices because each element of them is randomly generated.

1) Gaussian Sketch: The $p \times q$ matrix $\mathbf{S}$ reads

$$ \mathbf{S}(i, j) = \frac{1}{\sqrt{q}} G_{ij} \quad (15) $$

where $G_{ij}$ is identically and independently sampled from a standard normal distribution $N(0, 1)$.

2) Bernoulli/Rademacher Sketch: The $p \times q$ Bernoulli sketching matrix $\mathbf{S}$ is

$$ \mathbf{S}(i, j) = \frac{1}{\sqrt{q}} B_{ij} \quad (16) $$

where $B_{ij}$ is identically and independently sampled from a symmetric Bernoulli distribution $\mathrm{Pr}(B_{ij} = \pm 1) = (1/2)$.

3) Random Phase Sketch: In this case, the $p \times q$ sketching matrix $\mathbf{S}$ is defined as

$$ \mathbf{S}(i, j) = \frac{1}{\sqrt{q}} e^{i\varphi_{ij}} \quad (17) $$

where $\varphi_{ij}$ are random numbers in the interval $[0, 2\pi]$.

4) Count Sketch: In this case, the $p \times q$ matrix has a single randomly chosen non-zero entry in each row [55]. Each non-zero entry is a Rademacher random variable $\pm 1$ with probability $1/2$. Intuitively, the product $\mathbf{BS}$ for count sketch matrix $\mathbf{S}$ can be computed in three steps [56]. First, hash each source with a discrete value uniformly sampled from $\{1, 2, \dots, q\}$. Second, multiply each source by a Rademacher random variable $\pm 1$ with probability $1/2$. Third, sum up the sources with the same hash value. This is similar to a random sampling of the sources as in Section IV-B1. However, the advantage of count sketching is that all of the sources are used
<figure><img src="image_2.png" /><figcaption>Fig. 1. Samples of sketching matrices for $p = 50$ and $q = 20$. The color ranges between -1 (black) and +1 (white).</figcaption></figure>

at each iteration, but unlike the Gaussian, Bernoulli, and phase random projections, each blended source is a combination of a subset of the physical sources.

## B. Randomized Orthogonal Matrices

In order to define a randomized orthogonal matrix, we begin by an orthogonal matrix $\mathbf{H}$ of size $p \times p$ then sketching matrix $\mathbf{S}$ is defined as

$$ \mathbf{S} = \Pi \mathbf{H} \mathbf{R} \quad (18) $$

where each column of the $p \times q$ matrix $\mathbf{R}$ is chosen randomly from the set of all $p$ canonical basis vectors (the columns of the identity matrix). The matrix $\Pi$ can be a permutation matrix that reorders the rows of the matrix $\mathbf{H}$ randomly, or it can be a diagonal matrix with diagonal entries equal to +1 or -1 with probability 1/2 [41]. In what follows, we review some orthogonal matrices that are commonly used in applications.

1) **Identity Matrix:** The simplest choice for $\mathbf{H}$ is the identity matrix, $\mathbf{H} = \mathbf{I}$. In this case, in each iteration of the randomized IR-WRI, we use a number of $q$ sources that are selected uniformly at random from the $p$ physical sources.

2) **Discrete Fourier Transform:** The matrix $\mathbf{H}$ reads

$$ \mathbf{H}(i, j) = \frac{1}{\sqrt{q}} e^{\frac{(2\pi(i-1)(j-1))}{p}}, \quad i, j = 1, \dots, p. \quad (19) $$

3) **Discrete Cosine Transform:** The matrix $\mathbf{H}$ reads

$$ \mathbf{H}(i, j) = \begin{cases} \frac{1}{\sqrt{p}} \cos\left(\frac{\pi(2i-1)(j-1)}{2p}\right), & \text{if } j=1 \\ \sqrt{\frac{2}{p}} \cos\left(\frac{\pi(2i-1)(j-1)}{2p}\right), & \text{if } j \neq 1 \end{cases} \quad (20) $$

$i, j = 1, \dots, p$. The DCT is a widely used tool for signal compression.

4) **Discrete Wavelet Transform:** It is another effective transform for signal compression and multiscale analysis [57]. The basis functions for DWT associated with the quadrature mirror filters $h(k)$ and $g(k)$ are constructed at the scale $j \in \mathbb{Z}$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Gram matrices associated with the samples of sketching matrices shown in Fig. 1 ($p = 50$ and $q = 20$). The color scale ranges between -1 (black) and +1 (white).</figcaption></figure>

and location $k \in \mathbb{Z}$ as $\varphi_{j_0,k}(x) = (2^{j_0})^{1/2}\varphi(2^{j_0}x - k)$ and $\psi_{j,k}(x) = (2^j)^{1/2}\psi(2^j x - k)$ where the scaling function $\varphi(x)$ and the wavelet function $\psi(x)$ are obtained by solving the following two-scale difference equations:

$$ \left\{ \begin{aligned} \varphi(x) &= \sqrt{2} \sum_{k} h(k) \varphi(2x-k) \\ \psi(x) &= \sqrt{2} \sum_{k} g(k) \varphi(2x-k). \end{aligned} \right. \quad (21) $$

The DWT matrix requires $p$ to be a power of 2. For $p$ that is not a power of 2, one may take a random subsample of rows and columns of the matrix that is a power of 2.

5) **Noiselet Transform:** It is originally presented in [43] as a measuring matrix for compressed sensing of signals that have a sparse representation in the Haar wavelet domain because noiselets are maximally incoherent to the Haar basis and have a fast implementation algorithm. The noiselet matrix of order $p$ is defined recursively by

$$ \mathbf{H}_p = \frac{1-\ell}{2} \Pi \begin{bmatrix} \ell \mathbf{H}_{p/2} & \mathbf{H}_{p/2} \\ \mathbf{H}_{p/2} & \ell \mathbf{H}_{p/2} \end{bmatrix} \quad (22) $$

beginning with $\mathbf{H}_1 = \mathbf{1}$, where $\Pi$ is a permutation matrix such that for $\mathbf{y} = \Pi \mathbf{x}$ we have $\mathbf{y}(2i - 1) = \mathbf{x}(i)$ and $\mathbf{y}(2i) = \mathbf{x}(i + (n/2))$ for $i = 1, 2, \dots, (p/2)$. For example, the matrix $\Pi$ for $p = 8$ can be explicitly expressed as

$$ \Pi = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}. \quad (23) $$

The noiselet thus requires $p$ to be a power of 2. For $p$ that is not a power of 2, one may take a random subsample of rows and columns of the matrix that is a power of 2.

6) **Hadamard Transform:** The Walsh-Hadamard matrix of order $p$ is defined recursively by [44] as

$$ \mathbf{H}_p = \frac{1}{\sqrt{2}} \begin{bmatrix} \mathbf{H}_{p/2} & \mathbf{H}_{p/2} \\ \mathbf{H}_{p/2} & -\mathbf{H}_{p/2} \end{bmatrix} \quad (24) $$

beginning with $\mathbf{H}_1 = \mathbf{1}$. This matrix also requires $p$ to be a power of 2. For $p$ that is not a power of 2, one may take a random subsample of rows and columns of the matrix that is a power of 2.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Error between the identity matrix and truncated expectation of sketching matrices for $p = 50$ and $q = 20$ versus level of truncation $N$ (number of realizations).</figcaption></figure>

power of 2. Similar to the Bernoulli and count sketch matrices, the entries of the Walsh–Hadamard matrix are either –1 or +1. However, the rows and columns of the Walsh–Hadamard matrix are mutually orthogonal.

Figs. 1 and 2 compare samples of sketching matrices $\mathbf{S}$ and associated Gram matrices $\mathbf{SS}^T$ for the aforementioned methods. Since we require the expectation of the Gram matrix to approach the identity matrix, the associated error curve for different sketching matrices is plotted in Fig. 3. Except for the identity and DWT, other matrices are well behaved with a little difference in their convergence rate.

## V. NUMERICAL EXAMPLES

We assess the performance of the IR-WRI method with randomized source sketching for three 2-D mono-parameter benchmark models: Society of Exploration Geophysicists (SEG)/European Association of Geoscientists and Engineers (EAGE) Overthrust [58], 2004 BP model [59], and 2007 BP model created by Hemang Shah from BP Exploration Operation Company in 2007. Frequency-domain finite-difference modeling with a 9-point stencil and perfectly matched layer (PML) [46] is used for discretizing the wave equation in this article. Also, a fixed penalty parameter is selected for all the numerical tests following the instruction given in [11]. Additionally, we utilized bound constraints introduced in [45] and [60] with the minimum and maximum of the true model as lower and upper bounds. Furthermore, in order to make the performance of different (deterministic and stochastic) methods comparable, we normalize the sketching matrix $\mathbf{S}$ at each iteration such that its maximum singular value is unity. We also utilize the traditional plane-wave encoding. It has been applied to both FWI and WRI in a deterministic setting [21], [23]–[25]. It is based on transforming the source dimension by the linear Radon transform and then inverting the coefficients associated with a predefined range of angles. The plane wave encoding can be properly formulated with a fixed sketching matrix (see Appendix). However, we found that, for different models, satisfactory results can be obtained only when the coefficients span over a wide range of angles, which translates to a large number of plane-wave sources. In order to increase the computational efficiency while maintaining the reconstruction accuracy, we use a randomized plane-wave

<figure><img src="image_5.png" /><figcaption>Fig. 4. (Top) SEG/EAGE overthrust velocity model. (Middle) 2004 BP velocity model. (Bottom) 2007 BP velocity model. (Left column) True model and (right column) initial models used for the inversion.</figcaption></figure>

encoding, in which only $q \ll p$ random plane-wave sources are used at each iteration.

### A. SEG/EAGE Overthrust Example

We assess the performance of IR-WRI with randomized source sketching against a 2-D section of the 3-D SEG/EAGE overthrust model [see Fig. 4 (top left)]. The model includes 801 × 187 grid points with a grid interval of 25 m in both spatial directions. The surface acquisition is conducted with $p = 134$ sources and $m = 801$ receivers with an interval of 150 and 25 m, respectively. A smooth version of the true velocity model is used as the initial model [Fig. 4 (top right)], and the inversion is performed over two frequency paths 3–6.5 and 3–13 Hz with a frequency interval of 0.5 Hz in which the final estimated model of the first path is used as an initial model for the second path. Moreover, the maximum number of iterations per frequency is set to 10, except for the first two frequencies for which we perform 30 iterations. We perform the inversion with both deterministic and randomized IR-WRI. We remind that all $p$ sources are inverted at each iteration in the deterministic case, while in the randomized IR-WRI, only $q$ sketched sources are inverted at each iteration. The randomized algorithm reduces to the deterministic algorithm if the sketching matrix is identity, $\mathbf{S} = \mathbf{I}$. We consider the deterministic inversion result as a reference model to assess the results obtained by the randomized inversion. Figs. 5 and 6 show the inversion result for $q = 2$ and $q = 10$, respectively, where the deterministic result (with $p = 134$) is shown for comparison. We see that even with $q = 2$, the randomized inversion performs satisfactorily with most of the sketching matrices. Among them, DCT performed best, and the identity and DWT performed worst (of all). Fig. 7 shows the evolution of the ME (averaged over ten realizations) of different methods for a range of $q$ values between 2 and 20. We observe that for large $q$, all error curves approach the deterministic error (shown by the horizontal dashed line). Furthermore, sketching with DCT, DFT, noiselet, count, Gaussian, and Bernoulli performed almost similarly, and their quality is nearly equal to that of deterministic inversion for $q > 12$.
<figure><img src="image_6.png" /><figcaption>Fig. 5. Overthrust model. Inversion results from different methods with $p = 134$ and $q = 2$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Overthrust model. Inversion results from different methods with $p = 134$ and $q = 10$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. ME curve versus $q$ for different methods associated with the overthrust model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Same as Fig. 5 but (BM3D) regularization applied.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Same as Fig. 6 but with noisy data.</figcaption></figure>

Comparing the results obtained for $q = 2$ and $q = 10$ (Figs. 5 and 6) shows that the former results suffer from crosstalk noises which mainly takes the form of random noise in the final model. In order to suppress such noise, we applied the adaptive block-matching and 3-D filtering (BM3D) regularization [51], [61] for the case $q = 2$. The regularized inversion results are shown in Fig. 8. We observe that the crosstalk noise has been significantly attenuated, showing that we can use a very low number of $q$ while maintaining the quality of the reconstruction with appropriate regularization.

We continued by examining the effect of random noise in the data on the randomized sketching algorithm. For this, the data are contaminated by random noises of signal-to-

noise ratio (SNR) 10 dB and inverted with different methods while we set $q = 10$. The inversion results are shown in Fig. 9. We see that the randomized inversion stably generated results with the same quality as those of the deterministic algorithm. In the next section, we test the algorithm with the challenging 2004 BP salt model.

## B. 2004 BP Example

Here, we assess the performance of IR-WRI with randomized source sketching against the challenging 2004 BP salt model [see Fig. 4 (middle left)]. The model includes 100 × 369 grid points with a grid interval of 30 m in both spatial directions. The surface acquisition is conducted with equally spaced $p = 62$ sources and $m = 185$ receivers. A smooth version of the true velocity model is used as the initial model [see Fig. 4, (middle right)], and the inversion is performed over two frequency paths 3–8 and 3–12 Hz with a frequency interval of 0.5 Hz in which the final estimated model of the first path is used as an initial model for the second path. Moreover, the maximum number of iterations per frequency is set to 10, except for the first two frequencies for which we perform 25 and 15 iterations, respectively. We perform the inversion with both deterministic and randomized IR-WRI. Figs. 10 and 11 show the randomized inversion results for $q = 2$ and $q = 10$, respectively, while the deterministic result (with $p = 62$) is shown for comparison. Again, we see that the randomized inversion with $q = 2$ was able to recover the main structures of the model for all of the sketching matrices while all of results obtained with $q = 10$ have the same quality as those obtained with the deterministic inversion. Fig. 12 shows the evolution of the ME (averaged over ten realizations)
<figure><img src="image_11.png" /><figcaption>Fig. 10. Inversion results for the 2004 BP model for different methods with $p = 62$ and $q = 2$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Inversion results for the 2004 BP model for different methods with $p = 62$ and $q = 10$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. ME curve versus $q$ for different methods associated with the 2004 BP model.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Same as Fig. 10 but with (BM3D) regularization applied.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Inversion results for the 2007 BP model for different methods with $p = 263$ and $q = 2$.</figcaption></figure>

of different methods for a range of $q$ values between 2 and 20. Again, the quality of all methods (except the identity) is nearly equal to that of the deterministic inversion for $q > 12$. We also apply the regularization for the case $q = 2$ (see Fig. 13). Comparing this figure with Fig. 10 clearly shows the positive impact of regularization in reducing the crosstalk noise. In the next section, we test the algorithm with the 2007 BP salt model.

### C. 2007 BP Example

Finally, we assess the performance of IR-WRI with randomized source sketching against the large 2007 BP model (Fig. 4, bottom left). The model includes 151 × 1050 grid points with

a grid interval of 75 m in both spatial directions. The surface acquisition is conducted with equally spaced $p = 263$ sources and $m = 1050$ receivers. A smooth version of the true velocity model is used as the initial model [see Fig. 4 (bottom right)], and the inversion is performed over three frequency paths 1–2, 1–3.5, and 1–4 Hz with a frequency interval of 0.5 Hz in which the final estimated model of each path is used as an initial model for the next path. Moreover, the maximum number of iterations per frequency is set to 10, except for the first frequency for which we perform 20 iterations. We perform the inversion by both deterministic and randomized IR-WRI. Figs. 14 and 15 show the inversion results for $q = 2$ and $q = 10$, respectively, while the deterministic inversion result (with $p = 263$) is shown for comparison. This model is more challenging to reconstruct. Generally, the quality of the reconstructions for $q = 2$ is poorer compared to the previous benchmarks. The result obtained with the randomized identity sketching matrix is meaningless. The other methods recover the top of the central salt body, but the left salt is poorly recovered. For $q = 10$ (see Fig. 15), however, all methods performed satisfactorily. We can see from the ME (averaged over ten realizations) curves in Fig. 16 that except identity, Hadamard, and DWT, the error of other methods approaches that of the deterministic for $q > 10$. We also see that the randomized noiselet performs even better than the original deterministic inversion for $q > 8$. Randomized plane wave sketching shows similar behavior for $q > 16$.

In order to examine the effect of regularization for improving the quality of reconstruction for the 2007 BP model, we implemented all the methods for $q = 2$ with BM3D
<figure><img src="image_16.png" /><figcaption>Fig. 15. Inversion results for the 2007 BP model for different methods with $p = 263$ and $q = 10$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. ME curve versus $q$ for different methods associated with the 2007 BP model.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Same as Fig. 14 but with (BM3D) regularization.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Same as Fig. 15 but with (BM3D) regularization.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. Evolution of ME curve versus IR-WRI iteration for different methods in recovering the inversion results in Fig. 18.</figcaption></figure>

<table><caption>TABLE I
TOTAL NUMBER OF PDES SOLVED AND SPEED UP (IN PERCENT) FOR THE DETERMINISTIC IR-WRI AND ITS RANDOMIZED SOURCE SKETCHING VARIANT</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Randomized Sketching</th></tr><tr><th>q=2</th><th>q=10</th></tr></thead><tbody><tr><td rowspan="2">Ovethrust</td><td>PDE Solved</td><td>41540</td><td>620</td><td>3100</td></tr><tr><td>Speed up (%)</td><td>-</td><td>98.50</td><td>92.53</td></tr><tr><td rowspan="2">2004 BP</td><td>PDE Solved</td><td>21080</td><td>680</td><td>3400</td></tr><tr><td>Speed up (%)</td><td>-</td><td>96.77</td><td>83.87</td></tr><tr><td rowspan="2">2007 BP</td><td>PDE Solved</td><td>42080</td><td>320</td><td>1600</td></tr><tr><td>Speed up (%)</td><td>-</td><td>99.23</td><td>96.20</td></tr></tbody></table>

regularization (see Fig. 17). For all of the sketching methods, regularization greatly improves the quality (compare Figs. 14 and 17).

For a more comprehensive assessment of the regularization, we also apply it with the randomized algorithms for $q = 10$. The inversion results after 160 iterations are shown in Fig. 18. Besides, the associated ME curves versus iteration are shown in Fig. 19. Interestingly, most of the randomized algorithms with regularization outperformed the original deterministic algorithm with regularization. This shows that even with less computational burden, the randomized algorithm can explore areas in model space that are not accessible to the deterministic counterpart. In Fig. 20, we demonstrate the best inversion results obtained by the randomized source sketching method

(with $q = 10$) along with true models and corresponding deterministic results. This representation shows the highest similarity between the deterministic inversion results and its randomized source sketching version.

## D. Computational Efficiency

The computational gain of randomized source sketching can be analyzed by the total number of PDEs solves. The total number of PDE solves for each method is equal to the number of inversion iterations multiplied by the number of (right-hand sides) source terms in the data assimilated system. For the deterministic and randomized algorithms, this value equals $N \times p$ and $N \times q$, respectively, wherein $N$ denotes the total number of iterations performed. Assuming an equal
<figure><img src="image_21.png" /><figcaption>Fig. 20. Comparison between deterministic inversion and randomized source sketching method. (Left column) true models, (middle column) inversion results obtained by the deterministic, and (right column) randomized source sketching method with $q = 10$.</figcaption></figure>

number of iterations, the speedup achieved by the sketching is given by

$$S = \left(1 - \frac{q}{p}\right) \times 100\%. \qquad (25)$$

The speedup gained for each example presented above is reported in Table I and is one to two order of magnitudes.

## VI. CONCLUSION

We improve the computational efficiency of an extended-space formulation of FWI based on wavefield reconstruction (IR-WRI) using randomized source sketching or encoding. The source dimension of the problem is projected to a lower-dimensional space by using a random sketching matrix, which is regenerated at each iteration. This allowed reducing the number of working sources significantly compared to the original problem and hence the number of PDE solutions. We showed that the proposed algorithm is quite general and includes the existing source encoding methods as special cases. Accordingly, we proposed a randomized variant of the traditional deterministic plane-wave encoding. We assess the performance of the proposed method with three 2-D benchmark models. We conclude that the proposed randomized algorithm reduces the computational cost of the wavefield inversion by more than one order of magnitude for noise-free and noisy data. Higher speedup may be obtained for 3-D acquisitions. We showed that sparsity-promoting regularization can reduce the crosstalk noise which arise when the source dimension is reduced significantly. Furthermore, regularization makes the randomized algorithm able to explore areas of the model space that are not accessible to the deterministic counterpart.

# APPENDIX

## PLANE WAVE ENCODING

We review the plane wave theory (linear radon transform) and construction of the associated projection matrix. It has

been used for encoding seismic sources in FWI and WRI problems in a deterministic setting [21], [23]–[25]. For a time domain seismic data $d(t, x)$ with spatial variable $x$ (e.g., the source coordinate and in this case $d(t, x)$ will be a common receiver gather) and time $t$, the linear Radon transform of $d$ can be defined by the following integral [62]:

$$\hat{d}(\tau, p) = \int d(t = \tau + \bar{p}x, x)dx \qquad (26)$$

where $\tau$ is the intercept time and $\bar{p}$ is a ray parameter. The integration is computed over the range of $x$. Due to the time-invariant nature of this integral equation, it can be computed for each frequency slice $\omega$ separately. After a proper discretization of the equation, we get for a single-frequency slice of the Radon coefficient

$$\hat{d}(\omega, \bar{p}) = \sum_{x} d(\omega, x) e^{i\omega \bar{p} x}. \qquad (27)$$

Accordingly, the AL (12) in the plane-wave domain is obtained by defining the sketch matrix $S$ as

$$S(i, j) = \begin{cases} e^{i\omega\bar{p}(j)(x_i(i) - x_{\min})}, & \text{if } \bar{p}(j) \ge 0 \\ e^{i\omega\bar{p}(j)(x_i(i) - x_{\max})}, & \text{if } \bar{p}(j) < 0 \end{cases} \qquad (28)$$

[24], where $\bar{p}$ denotes the ray parameter vector and $x_{\min}/x_{\max}$ denote the minimum/maximum offset. Using the encoding matrix (28), the $p$ original sources are combined linearly to be replaced by $q$ plane-wave sources, where $q$ is the number of ray parameters. The success of the plane-wave encoding relies on the range and number of ray parameters which can be determined by the following condition [23]:

$$q \ge x_{\max} f_{\max} (\bar{p}_{\max} - \bar{p}_{\min}) \qquad (29)$$

where $q$ is the minimum required number of ray parameters, $x_{\max}$ and $f_{\max}$ denote the maximum offset and frequency, respectively, and $\bar{p}_{\max}/\bar{p}_{\min}$ denote the maximum/minimum ray parameter.
ACKNOWLEDGMENT

The authors are grateful to the OPAL infrastructure from Observatoire de la Côte d'Azur (CRIMSON) for providing resources and support. This work was granted access to the HPC resources of IDRIS under the allocation A0050410596 made by GENCI.

REFERENCES

[1] R. G. Pratt, C. Shin, and G. J. Hick, “Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion,” *Geophys. J. Int.*, vol. 133, no. 2, pp. 341–362, May 1998.

[2] T. Ha, S. Pyun, and C. Shin, “Efficient electric resistivity inversion using adjoint state of mixed finite-element method for Poisson’s equation,” *J. Comput. Phys.*, vol. 214, no. 1, pp. 171–186, May 2006.

[3] F. Lavoué, R. Brossier, L. Métivier, S. Garambois, and J. Virieux, “Two-dimensional permittivity and conductivity imaging by full waveform inversion of multiooffset GPR data: A frequency-domain quasi-Newton approach,” *Geophys. J. Int.*, vol. 197, pp. 248–268, Apr. 2014.

[4] A. Tarantola, “Inversion of seismic reflection data in the acoustic approximation,” *Geophysics*, vol. 49, no. 8, pp. 1259–1266, 1984.

[5] L. Sirgue and R. G. Pratt, “Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies,” *Geophysics*, vol. 69, no. 1, pp. 231–248, Jan. 2004.

[6] J. Virieux and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” *Geophysics*, vol. 74, no. 6, pp. WCC1–WCC26, Nov. 2009.

[7] E. Haber, U. M. Ascher, and D. Oldenburg, “On optimization techniques for solving nonlinear inverse problems,” *Inverse Prob.*, vol. 16, no. 5, p. 1263, 2000.

[8] J. Virieux, A. Asnaashari, R. Brossier, L. Métivier, A. Ribodetti, and W. Zhou, “An introduction to full waveform inversion,” in *Encyclopaedia Exploration Geophysics*. Tulsa, OK, USA: Society of Exploration Geophysicists, 2017, p. R1.

[9] P. M. van den Berg and R. E. Kleinman, “A contrast source inversion method,” *Inverse Prob.*, vol. 13, no. 6, p. 1607, Dec. 1997.

[10] T. van Leeuwen and F. J. Herrmann, “Mitigating local minima in full-waveform inversion by expanding the search space,” *Geophys. J. Int.*, vol. 195, no. 1, pp. 661–667, 2013.

[11] H. S. Aghamiry, A. Gholami, and S. Operto, “Improving full-waveform inversion by wavefield reconstruction with the alternating direction method of multipliers,” *Geophysics*, vol. 84, no. 1, pp. R125–R148, Jan. 2019.

[12] D. P. Bertsekas, “On penalty and multiplier methods for constrained minimization,” *SIAM J. Control Optim.*, vol. 14, no. 2, pp. 216–235, Feb. 1976.

[13] I. S. Duff, A. M. Erisman, and J. K. Reid, *Direct Methods for Sparse Matrices*, 2nd ed. Oxford, U.K.: Clarendon Press, 2017.

[14] Y. Saad, *Iterative Methods for Sparse Linear Systems*. Philadelphia, PA, USA: SIAM, 2003.

[15] M. L. Parks, E. de Sturler, G. Mackey, D. D. Johnson, and S. Maiti, “Recycling Krylov subspaces for sequences of linear systems,” *SIAM J. Sci. Comput.*, vol. 28, no. 5, pp. 1651–1674, Oct. 2006.

[16] P. Jolivet and P.-H. Tournier, “Block iterative methods and recycling for improved scalability of linear solvers,” in *Proc. Int. Conf. High Perform. Comput., Netw., Storage Anal. (SC)*, Nov. 2016, pp. 190–203.

[17] X. Li, A. Y. Aravkin, T. van Leeuwen, and F. J. Herrmann, “Fast randomized full-waveform inversion with compressive sensing,” *Geophysics*, vol. 77, no. 3, pp. A13–A17, May 2012.

[18] T. M. Habashy, A. Abubakar, G. Pan, and A. Belani, “Source-receiver compression scheme for full-waveform seismic inversion,” *Geophysics*, vol. 76, no. 4, pp. R95–R108, Jul. 2011.

[19] C. Castellanos, L. Métivier, S. Operto, R. Brossier, and J. Virieux, “Fast full waveform inversion with source encoding and second-order optimization methods,” *Geophys. J. Int.*, vol. 200, no. 2, pp. 720–744, Feb. 2015.

[20] Q. Zhang, W. Mao, H. Zhou, H. Zhang, and Y. Chen, “Hybrid-domain simultaneous-source full waveform inversion without crosstalk noise,” *Geophys. J. Int.*, vol. 215, no. 3, pp. 1659–1681, 2018.

[21] D. Vigh and E. W. Starr, “3D prestack plane-wave, full waveform inversion,” *Geophysics*, vol. 73, pp. VE135–VE144, 2008.

[22] N. D. Whitmore, “An imaging hierarchy for common-angle plane wave seismograms,” Ph.D. dissertation, Graduate School, Univ. Tulsa, Tulsa, OK, USA, 1995.

[23] Y. Zhang, J. Sun, C. Notfors, S. H. Gray, L. Chernis, and J. Young, “Delayed-shot 3D depth migration,” *Geophysics*, vol. 70, no. 5, pp. E21–E28, Sep. 2005.

[24] Y. Tao and M. K. Sen, “Frequency-domain full waveform inversion with plane-wave data,” *Geophysics*, vol. 78, no. 1, pp. R13–R23, Jan. 2013.

[25] S. Kim, W. Chung, Y. S. Kim, and C. Shin, “Seismic wavefield reconstruction inversion using a plane-wave encoding strategy,” *J. Geophys. Eng.*, vol. 17, no. 6, pp. 1037–1048, Dec. 2020.

[26] F. Liu, D. W. Hanson, N. D. Whitmore, R. S. Day, and R. H. Stolt, “Toward a unified analysis for source plane-wave migration,” *Geo-physics*, vol. 71, no. 4, pp. S129–S139, Jul. 2006.

[27] L. A. Romero, D. C. Ghiglia, C. C. Ober, and S. A. Morton, “Phase encoding of shot records in prestack migration,” *Geophysics*, vol. 65, no. 2, pp. 426–436, Mar. 2000.

[28] J. R. Krebs et al., “Fast full-wavefield seismic inversion using encoded sources,” *Geophysics*, vol. 74, no. 6, pp. WCC177–WCC188, Nov. 2009.

[29] A. J. Berkhout and G. Blacquière, “Blended acquisition with dispersed source arrays, the next step in seismic acquisition?” in *Proc. SEG Tech. Program Expanded Abstr.*, Jan. 2011, pp. 16–19.

[30] H. Ben-Hadj-Ali, S. Operto, and J. Virieux, “An efficient frequency-domain full waveform inversion method using simultaneous encoded sources,” *Geophysics*, vol. 76, no. 4, pp. R109–R124, Jul. 2011.

[31] A. Schiemenz and H. Igel, “Accelerated 3-D full-waveform inversion using simultaneously encoded sources in the time domain: Application to Valhall ocean-bottom cable data,” *Geophys. J. Int.*, vol. 195, no. 3, pp. 1970–1988, Dec. 2013.

[32] G. T. Schuster, X. Wang, Y. Huang, W. Dai, and C. Boonyasiriwat, “Theory of multisource crosstalk reduction by phase-encoded statics,” *Geophys. J. Int.*, vol. 184, no. 3, pp. 1289–1303, Mar. 2011.

[33] E. Haber, M. Chung, and F. Herrmann, “An effective method for parameter estimation with PDE constraints with multiple right-hand sides,” *SIAM J. Optim.*, vol. 22, no. 3, pp. 739–757, Jul. 2012.

[34] C. Boonyasiriwat and G. T. Schuster, “3D multisource full-waveform inversion using dynamic random phase encoding,” in *Proc. SEG Tech. Program Expanded Abstr.*, Jan. 2010, pp. 1044–1049.

[35] Z. Xue, Y. Chen, S. Fomel, and J. Sun, “Seismic imaging of incomplete data and simultaneous-source data using least-squares reverse time migration with shaping regularization,” *Geophysics*, vol. 81, no. 1, pp. S11–S20, Jan. 2016.

[36] Y. Chen, H. Chen, K. Xiang, and X. Chen, “Preserving the discontinuities in least-squares reverse time migration of simultaneous-source data,” *Geophysics*, vol. 82, no. 3, pp. S185–S196, May 2017.

[37] A. Guitton and E. Díaz, “Attenuating crosstalk noise with simultaneous source full waveform inversion,” *Geophys. Prospecting*, vol. 60, no. 4, pp. 759–768, 2012.

[38] A. Y. Anagaw and M. D. Sacchi, “Full waveform inversion with simultaneous sources using the full Newton method,” in *Proc. SEG Tech. Program Expanded Abstr.*, Sep. 2012, pp. 1–5.

[39] Y. Rao and Y. Wang, “Seismic waveform tomography with shot-encoding using a restarted L-BFGS algorithm,” *Sci. Rep.*, vol. 7, no. 1, pp. 1–9, Dec. 2017.

[40] A. Y. Anagaw and M. D. Sacchi, “Comparison of multifrequency selection strategies for simultaneous-source full-waveform inversion,” *Geophysics*, vol. 79, no. 5, pp. R165–R181, Sep. 2014.

[41] M. Pilanci and M. J. Wainwright, “Randomized sketches of convex programs with sharp guarantees,” *IEEE Trans. Inf. Theory*, vol. 61, no. 9, pp. 5096–5115, Sep. 2015.

[42] R. M. Gower and P. Richtárik, “Randomized quasi-Newton updates are linearly convergent matrix inversion algorithms,” *SIAM J. Matrix Anal. Appl.*, vol. 38, no. 4, pp. 1380–1409, Jan. 2017.

[43] R. Coifman, F. Geshwind, and Y. Meyer, “Noiselets,” *Appl. Comput. Harmon. Anal.*, vol. 10, no. 1, pp. 27–44, 2001.

[44] R. K. Yarlagadda and J. E. Hershey, *Hadamard Matrix Analysis and Synthesis: With Applications to Communications and Signal/Image Processing*, vol. 383. New York, NY, USA: Springer, 2012.

[45] H. S. Aghamiry, A. Gholami, and S. Operto, “Compound regularization of full-waveform inversion for imaging piecewise media,” *IEEE Trans. Geosci. Remote Sens.*, vol. 58, no. 2, pp. 1192–1204, Feb. 2020.

[46] Z. Chen, T. Wu, and H. Yang, “An optimal 25-point finite difference scheme for the Helmholtz equation with PML,” *J. Comput. Appl. Math.*, vol. 236, no. 6, pp. 1240–1258, Oct. 2011.

[47] R.-E. Plessix, “A review of the adjoint-state method for computing the gradient of a functional with geophysical applications,” *Geophys. J. Int.*, vol. 167, no. 2, pp. 495–503, Nov. 2006.
[48] M. J. Powell, "A method for nonlinear constraints in minimization problems," Optim., pp. 283-298, 1969.

[49] S. Boyd, N. Parikh, and E. Chu, *Distributed Optimization and Statistical Learning Via the Alternating Direction Method of Multipliers*. Norwell, MA, USA: Now Publishers, 2011.

[50] A. Gholami and S. M. Hosseini, "A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals," *Signal Process.*, vol. 93, no. 7, pp. 1945-1960, 2013.

[51] H. S. Aghamiry, A. Gholami, and S. Operto, "Full waveform inversion by proximal Newton method using adaptive regularization," *Geophys. J. Int.*, vol. 224, no. 1, pp. 169-180, Nov. 2020, doi: 10.1093/gji/ggaa434.

[52] M. F. Hutchinson, "A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines," *Commun. Statist. Simul. Comput.*, vol. 18, no. 3, pp. 1059-1076, Jan. 1989.

[53] T. van Leeuwen, A. Y. Aravkin, and F. J. Herrmann, "Seismic waveform inversion by stochastic optimization," *Int. J. Geophy.*, vol. 2011, pp. 1-18, Mar. 2011.

[54] M. Louboutin and F. J. Herrmann, "Ultra-low memory seismic inversion with randomized trace estimation," 2021, arXiv:2104.00794.

[55] M. Charikar, K. Chen, and M. Farach-Colton, "Finding frequent items in data streams," *Theor. Comput. Sci.*, vol. 312, no. 1, pp. 3-15, 2004.

[56] S. Wang, "A practical guide to randomized matrix computations with MATLAB implementations," 2015, arXiv:1505.07570.

[57] S. Mallat, *A Wavelet Tour of Signal Processing*. Amsterdam, The Netherlands: Elsevier, 1999.

[58] F. Aminzadeh, J. Brac, and T. Kunz, "3-D salt and overthrust models," in *SEG/EAGE 3-D Modeling Series*, no. 1. USA: Soc. Explor. Geophysicists, 1997.

[59] F. J. Billette and S. Brandsberg-Dahl, "The 2004 BP velocity benchmark," in *Proc. 67th EAGE Conf. Exhibit.*, 2005, p. B035.

[60] 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," *Geophys. J. Int.*, vol. 218, no. 2, pp. 855-872, 2019.

[61] A. Danielyan, V. Katkovnik, and K. Egiazarian, "BM3D frames and variational image deblurring," *IEEE Trans. Image Process.*, vol. 21, no. 4, pp. 1715-1728, Apr. 2012.

[62] G. Beylkin, "Discrete Radon transform," *IEEE Trans. Acoust., Speech, Signal Process.*, vol. ASSP-35, no. 2, pp. 162-172, Feb. 1987.

<figure><img src="image_22.png" /></figure>

**Kamal Aghazade** received the B.Sc. degree in physics from the University of Tabriz, Tabriz, Iran, in 2013, and the M.Sc. degree in geophysics from the University of Tehran, Tehran, Iran, in 2015, where he is currently pursuing the Ph.D. degree.

His research interests include linear and nonlinear optimization and their applications in seismic inversion.

Mr. Aghazade is a member of the Society of Exploration Geophysicists (SEG) and the European Association of Geoscientists and Engineers (EAGE).

<figure><img src="image_23.png" /></figure>

**Hossein S. Aghamiry** (Member, IEEE) received the Ph.D. degree in seismology from the Institute of Geophysics, University of Tehran, Tehran, Iran, and the Ph.D. degree in earth sciences from the University Côte d'Azur, Valbonne, France, in the framework of joint Ph.D. student, in 2019.

He is currently a Post-Doctoral Fellow with Geoazur, Valbonne. His research interests include linear and nonlinear optimization and their applications in seismic processing, imaging, and inversion.

Dr. Aghamiry is a member of the Society of Exploration Geophysicists (SEG), the European Association of Geoscientists and Engineers (EAGE), and the Society for Industrial and Applied Mathematics (SIAM).

<figure><img src="image_24.png" /></figure>

**Ali Gholami** received the Ph.D. degree in geophysics from the University of Tehran, Tehran, Iran, in 2010.

He is currently a Professor of geophysics with the Institute of Geophysics, University of Tehran. His research interests include inverse problems, statistical signal processing, compressed sensing, regularization methods, and application of these methods to seismic processing, imaging, and inversion.

<figure><img src="image_25.png" /></figure>

**Stéphane Operto** received the degree in engineering from ENSAM, Paris, France, in 1990, and the Ph.D. degree in marine geophysics from the University Pierre et Marie Curie (Paris 6), Paris, in 1995.

He was a Post-Doctoral Fellow at the Institute for Geophysics, The University of Texas at Austin, Austin, TX, USA, in 1995, and at the Paris School of Mines, Institute for Geophysics, Paris, in 1996. In 1998, he joined CNRS as a Research Scientist. His professional interests are wave-propagation modeling and full-waveform inversion (FWI).

Dr. Operto is a member of the Society of Exploration Geophysicists (SEG), the European Association of Geoscientists and Engineers (EAGE), the Society for Industrial and Applied Mathematics (SIAM), and the American Geophysical Union (AGU).