# Why POCS works, and how to make it better

Matteo Ravasi
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
&Nick Luiken
KAUST
Thuwal, Kingdom of Saudi Arabia
<EMAIL_ADDRESS>
###### Abstract

Projection Over Convex Sets (POCS) is one of the most widely used algorithms
in geophysical data processing to interpolate seismic data. Whilst usually
described as a modification of the Gerchberg-Saxton algorithm, a formal
understanding of the underlying objective function and its implication for the
associated optimization process is lacking to date in the literature. We show
that the POCS algorithm can be interpreted as the application of the Half-
Quadratic Splitting (HQS) method to the $L_{0}$ norm of the projected sought
after data (where the projection can be any orthonormal transformation),
constrained on the available traces. Similarly, the popular, apparently
heuristic strategy of using a decaying threshold in POCS is revealed to be the
result of the continuation strategy that HQS must employ to converge to a
solution of the minimizer. Finally, we prove that for off-the-grid receivers,
the POCS algorithms must operate in an inner-outer fashion (i.e., an inverse
problem is solved at each outer iteration). In light of our new theoretical
understanding, we suggest the use of a modern solver such as the Chambolle-
Pock Primal-Dual algorithm and show that this can lead to a new POCS-like
method with superior interpolation capabilities at nearly the same
computational cost of the industry-standard POCS method.

##  1 Introduction

Motivated by economical and practical constraints, geophysical data
acquisition campaigns typically consist of a limited portion of the sought
after regular and finely sampled grids of sources and receivers. As such,
reconstructing missing traces represents one of the key early steps in the
geophysical processing pipeline.

In 2006, [1] proposed a method for reconstructing irregularly sampled data
under the name of Projection Over Convex Sets (POCS). Owing to its intuitive
nature and simplicity of implementation, POCS has become the de-facto first
choice in the seismic industry. Since its introduction to the geophysical
community, several works have focused on improving its effectiveness by
considering different transformations that render seismic data to be sparse
[2, 3, 4] and relaxation strategies for the thresholding step [5, 6]. Recent
interest has also emerged in combining the POCS iterations with neural network
based transformations [7]. Nevertheless, a formal interpretation of POCS in
the context of convex optimization is currently lacking. Whilst several
authors suggest that POCS is a modified version of the Gerchberg–Saxton
iterative algorithm [8], they fail to clearly identify the underlying
objective function that this method is trying to optimize.

By leveraging recent advances in the field of convex optimization and with the
help of proximal operators, we prove that the POCS iterations can be formally
derived by applying the Half-Quadratic Splitting (HQS) algorithm of [9] to the
$L_{0}$ norm of the projected sought after data, constrained on the available
traces. We also theoretically justify why the POCS algorithm benefits from
applying a relaxation strategy on the thresholding parameter. After laying the
mathematical foundations, we introduce two alternative solvers that can be
used to solve the objective functions that underlie the POCS method. Such
solvers are shown to be more flexible than HQS in that they can more naturally
handle the inclusion of frequency-wavenumber masking operators and off-the-
grid receivers. Numerical examples on a complex 3D synthetic dataset highlight
how such solvers can lead to improved and sometimes faster convergence than
the currently used POCS algorithm.

## Theory

### Preliminaries

Let $\mathbf{X}\in\mathbb{R}^{N_{r,y}\times N_{r,x}\times N_{t}}$ be a tensor
containing seismic data regularly sampled along a two-dimensional spatial grid
of size $N_{r,y}\times N_{r,x}$, and $\mathbf{Y}\in\mathbb{R}^{N\times N_{t}}$
be the physically recorded seismic data at $N$ arbitrarily selected receiver
locations within the grid. Their vectorized versions,
$\mathbf{x}=vec(\mathbf{X})$ and $\mathbf{y}=vec(\mathbf{Y})$, are linked via

$\textbf{y}=\textbf{R}\textbf{x},$ (1)

where $\mathbf{R}\in\mathbb{R}^{NN_{t}\times N_{r,y}N_{r,x}N_{t}}$ is a
restriction operator that selects the $N$ available traces from the regularly
sampled data. An alternative definition of the modelling operator is also
provided, where
$\tilde{\mathbf{R}}=\mathbf{R}^{H}\mathbf{R}\in\mathbb{R}^{N_{r,y}N_{r,x}N_{t}\times
N_{r,y}N_{r,x}N_{t}}$ is a masking operator (i.e., diagonal matrix containing
1s at the available spatial locations and 0s at the unavailable spatial
locations). Similarly, we define
$\tilde{\mathbf{Y}}\in\mathbb{R}^{N_{r,y}\times N_{r,x}\times N_{t}}$ as the
recorded dataset where zero traces are interleaved to the physically recorded
traces at locations of missing receivers. Its vectorized version
$\tilde{\mathbf{y}}=vec(\tilde{\mathbf{Y}})$ is therefore linked to the vector
x via $\tilde{\textbf{y}}=\tilde{\textbf{R}}\textbf{x}$. Similar to the
modelling operator, the following relation holds
$\tilde{\textbf{y}}=\textbf{R}^{H}\textbf{y}$.

Finally, the proximal operator of a possibly non-smooth function $\tau
f(\textbf{x})$ is defined as:

$prox_{\tau
f}(\mathbf{u})=\underset{\mathbf{x}}{\mathrm{min}}\;f(\mathbf{x})+\frac{1}{2\tau}||\mathbf{x}-\mathbf{u}||_{2}^{2},$
(2)

where $\tau>0$. From equation 2 we can observe that evaluating the proximal
operator involves solving an optimization problem, which compromises between
minimizing the function itself and finding a solution close to $\mathbf{u}$.
The introduction of the $L_{2}$ norm of the difference between the vector
$\mathbf{u}$ and the solution $\mathbf{x}$ renders the functional strongly
convex, so that is has a unique minimizer for every vector $\mathbf{u}$. As we
will see in the following, proximal operators are a key component in convex
optimization and closed-form solutions exist for many commonly used functions
$f$.

### The POCS algorithm

The aim of any seismic interpolation algorithm is to reconstruct the regularly
sampled data $\mathbf{X}$ from the available traces $\mathbf{Y}$. However,
solving equation 1 is a severely ill-posed inverse problem and requires prior
knowledge about the expected solution. One way to introduce such prior
knowledge is via linear sparsifying transforms,
$\mathbf{S}\in\mathbb{R}^{N_{r,y}N_{r,x}N_{t}\times N_{s}}$, that provide a
compact representation (i.e., few non-zero coefficients) of the regularly
sampled seismic data in the transformed domain. Here $N_{s}$ is used to
indicate the number of coefficients in the transformed domain.

The original POCS algorithm of [1] leverages a specific type of sparsifying
transform, namely the multi-dimensional Fourier transform, and can be
compactly expressed by the following iterations:

$\textbf{x}^{k}=\tilde{\textbf{y}}+(\textbf{I}-\tilde{\textbf{R}})\textbf{S}\mathbb{T}_{\alpha}(\textbf{S}^{H}\textbf{x}^{k-1}),$
(3)

where $\mathbb{T}_{\alpha}(x)=xH(|x|-\alpha)$ is a hard-thresholding operation
controlled via the parameter $\alpha$, and $H$ is the Heaviside step function.
Given a vector x, this operation acts in an element-wise fashion such that
elements with an absolute value smaller than $\alpha$ are truncated to zero
whilst all other elements remain untouched. The POCS algorithm owes its
popularity to its simplicity of implementation and effectiveness; in fact, the
main iteration of the algorithm in (3) can be easily interpreted as follows:
transform the current solution to a sparse domain of choice, apply a
threshold, transform back to the original domain, extract only traces at
unavailable receivers and interleave them with the true traces at available
receivers.

Note that, the POCS algorithm is agnostic to the sparsifying transform S
provided that it satisfies the following conditions:
$\textbf{S}^{H}\textbf{S}=\textbf{I}_{N_{s}}$ and
$\textbf{S}\textbf{S}^{H}=\textbf{I}_{N_{r,y}N_{r,x}N_{t}}$ (i.e.,
orthogonormality) as shown below. Moreover, whilst our theory is developed for
3D seismic data, its extension to 5D data is trivial as explained in [10].
Finally, in practical scenarios physical receivers may be placed at locations
that do not lie exactly on top of the nominal grid: the restriction operator
must therefore be substituted by a spatial interpolator. The implications of
such a change for the POCS algorithm will be later discussed.

## Deriving POCS from HQS

In this section, a mathematically sound derivation of the POCS iterations is
provided drawing upon the theory of convex optimization. To begin with, let’s
define the following objective function:

$J=f(\textbf{x})+g(\textbf{x}),$ (4)

where $f$ and $g$ are any convex functions with a known proximal operator. The
HQS algorithm minimizes the functional in equatio 4 by introducing a splitting
of the form:

$J=f(\textbf{x})+g(\textbf{z})\;\text{s.t.}\;\textbf{z}=\textbf{x},$ (5)

and relaxing the constraint using the penalty method:

$J_{HQS}=f(\textbf{x})+g(\textbf{z})+\frac{\rho}{2}||\textbf{z}-\textbf{x}||_{2}^{2}.$
(6)

The associated unconstrained minimization is solved in an alternating fashion:

$\begin{split}\textbf{x}^{k}&=\underset{\textbf{x}}{\mathrm{argmin}}\>J_{HQS}(\textbf{x},\textbf{z}^{k-1})=\underset{\textbf{x}}{\mathrm{argmin}}\;f(\textbf{x})+\frac{\rho}{2}||\textbf{z}^{k-1}-\textbf{x}||_{2}^{2}=prox_{\frac{1}{\rho}f}(\textbf{z}^{k-1})\\\
\textbf{z}^{k}&=\underset{\textbf{z}}{\mathrm{argmin}}\>J_{HQS}(\textbf{x}^{k},\textbf{z})=\underset{\textbf{z}}{\mathrm{argmin}}\;g(\textbf{z})+\frac{\rho}{2}||\textbf{z}-\textbf{x}^{k}||_{2}^{2}=prox_{\frac{1}{\rho}g}(\textbf{x}^{k}).\end{split}$
(7)

It becomes obvious from equation 7 that HQS can operate with any combination
of $f$ and $g$ functions provided their proximal operator can be easily
computed.

Let’s now define the following constrained objective function and prove that
its minimization via HQS is equivalent to the POCS algorithm:

$J_{POCS}=||\textbf{S}^{H}\textbf{x}||_{0}\;\text{s.t.}\;\textbf{y}=\textbf{R}\textbf{x},$
(8)

where $||.||_{0}$ identifies the $L_{0}$ norm. Note that this is equivalent to
writing an unconstrained objective function of the form:

$J_{POCS}=||\textbf{S}^{H}\textbf{x}||_{0}+i_{\textbf{y}=\textbf{R}\textbf{x}},$
(9)

where $i_{\textbf{y}=\textbf{R}\textbf{x}}$, is the indicator function of the
Affine set
$\\{\textbf{x}\in\mathbb{R}^{N_{r,y}N_{r,x}N_{t}}|\textbf{y}=\textbf{R}\textbf{x}\\}$.
The optimization problem associated with the functional in equation 8 can be
interpreted as finding a sparse representation of the regularly sampled
seismic data x (by minimizing the $L_{0}$ norm in the transformed domain)
whilst exactly matching the recorded data. Let
$f(\textbf{x})=||\textbf{S}^{H}\textbf{x}||_{0}$ and
$g(\textbf{x})=i_{\textbf{y}=\textbf{R}\textbf{x}}$, the POCS functional in
equation 8) can be minimized via the HQS algorithm by simply repeatedly
applying their proximal operators. Starting from $f$ and by leveraging the
fact that S is an orthonormal operator, the proximal operator has a closed-
form solution [11],

$prox_{\frac{1}{\rho}f}(\textbf{x})=\textbf{S}\mathbb{T}_{\frac{1}{\rho}}(\textbf{S}^{H}\textbf{x}).$
(10)

Similarly, since $g$ is an indicator function over an affine set, its proximal
operator can be written as [12]

$\begin{split}prox_{\frac{1}{\rho}g}(\textbf{x})&=\textbf{x}-\textbf{R}^{H}(\textbf{R}\textbf{R}^{H})^{-1}(\textbf{R}\textbf{x}-\textbf{y})\\\
&=\textbf{x}-\textbf{R}^{H}(\textbf{R}\textbf{x}-\textbf{y})\\\
&=\tilde{\textbf{y}}+(\textbf{I}-\tilde{\textbf{R}})\textbf{x},\end{split}$
(11)

where we have used the following property of a restriction operator,
$\textbf{R}\textbf{R}^{H}=\textbf{I}_{NN_{t}}$. This proximal operator
performs exactly the operation referred in the POCS literature to as the
available trace re-inserting strategy. By combining equations 10 and 11
together in the iterations of HQS, and taking $\alpha=1/\rho$, we obtain the
iterations of the original POCS algorithm (equation 3). At this point, three
key observations can be made:

  * •

following the definition in [13], POCS is an algorithm which seeks a solution
that jointly satisfies a number of properties expressed in terms of sets (and
it can be only applied to convex sets). On the other hand, the $L_{0}$ norm
used in equation 8 to obtain the iterations of Abma’s POCS is not a convex
set. We have therefore discovered that, despite its name, this algorithm is
strictly-speaking not a POCS-type algorithm;

  * •

in order to ensure convergence to a minimizer of the original objective
function, HQS must use a continuation strategy, i.e.,
$\rho^{k}\rightarrow\infty$. Looking at equation 10, this leads to decreasing
the threshold of the hard-thresholding operation during iterations. This
finding explains why an effective strategy to improve the convergence of the
POCS algorithm is to exponentially decrease the thresholding applied during
the x-update in the transformed domain [14];

  * •

several transforms commonly used in geophysics are not orthonormal. A
prominent example of non-orthonormal transform is the Radon transform. In this
case, equation 10 does not hold and a closed-form solution for the proximal
operator does not exist. Whilst [15] use POCS-like iterations in conjunction
with a parabolic Radon transform to fill missing near offsets, their algorithm
does not formally solve the objective function in equation 8. We will see in
the next section how by simply swapping HQS with another solver, the objective
function in equation 8 can still be minimized also when a non-orthonormal
transform is chosen.

## Modifications to the POCS algorithm

Equipped with a solid mathematical understanding of the POCS algorithm, our
next step is to define some alternative ways to minimize the functional in
equation 8.

Two popular methods in the field of convex optimization for the minimization
of such composite objective functions are the Alternating Direction Method of
Multipliers (ADMM - [16]) and the Chambolle-Pock Primal-Dual algorithm (PD -
[17]). ADMM takes a similar route to that of HQS in that it performs a
splitting of the original objective function. Differently from HQS, ADMM
introduces an augmented Lagrangian of the form:

$J_{ADMM}=f(\textbf{x})+g(\textbf{z})+\textbf{y}^{H}(\textbf{z}-\textbf{x})+\frac{\rho}{2}||\textbf{z}-\textbf{x}||_{2}^{2},$
(12)

where y is the Lagrange multiplier. By solving this augmented objective
function in an alternating fashion for our definition of $f$ and $g$
functions, we obtain:

$\begin{split}\textbf{x}^{k}&=prox_{\frac{1}{\rho}f}(\textbf{z}^{k-1}-\textbf{y}^{k-1})=\textbf{S}\mathbb{T}_{\frac{1}{\rho}}(\textbf{S}^{H}(\textbf{z}^{k-1}-\textbf{y}^{k-1}))\\\
\textbf{z}^{k}&=prox_{\frac{1}{\rho}g}(\textbf{x}^{k}+\textbf{y}^{k-1})=\tilde{\textbf{y}}+(\textbf{I}-\tilde{\textbf{R}})(\textbf{x}^{k}+\textbf{y}^{k-1})\\\
\mathbf{y}^{k+1}&=\mathbf{y}^{k}+\mathbf{x}^{k+1}-\mathbf{z}^{k+1}.\end{split}$
(13)

The presence of a Lagrangian term equation 12 ensures that at convergence the
constraint $\textbf{z}=\textbf{x}$ is verified and the original problem is
solved. Theoretically, the scalar $\rho$ can therefore be kept fixed
throughout iterations, relieving the user from defining a cooling strategy as
required by HQS.

On the other hand, the Primal-Dual algorithm is designed to solve problems in
the form of $f(\textbf{Kx})+g(\textbf{x})$ where K is any linear operator. It
proceeds finding the saddle point of the functional:

$\max_{\textbf{y}}\min_{\textbf{x}}J_{PD}=\max_{\textbf{y}}\min_{\textbf{x}}\textbf{y}^{H}(\textbf{Kx})+g(\textbf{x})+f^{*}(\textbf{y}),$
(14)

where $f^{*}$ is the convex conjugate of $f$. This problem can be solved by
the following iterations:

$\begin{split}\mathbf{y}^{k+1}&=prox_{\mu
f^{*}}(\mathbf{y}^{k}+\mu\mathbf{K}\bar{\mathbf{x}}^{k})\\\
\mathbf{x}^{k+1}&=prox_{\tau
g}(\mathbf{x}^{k}-\tau\mathbf{K}^{H}\mathbf{y}^{k+1})\\\
\bar{\mathbf{x}}^{k+1}&=2\mathbf{x}^{k+1}-\mathbf{x}^{k}.\end{split}$ (15)

where $\tau$ and $\mu$ represent the step-lengths of the sub-gradients of $f$
and $g$. To ensure convergence $\tau\mu L^{2}<1$, where
$L^{2}=||\textbf{K}||_{2}^{2}=\lambda_{max}(\mathbf{K}^{H}\mathbf{K})$ is the
spectral radius of the linear operator K.

If we choose $\textbf{K}=\textbf{I}$,
$f(\textbf{x})=||\textbf{S}^{H}\textbf{x}||_{0}$, and
$g(\textbf{x})=i_{\textbf{y}=\textbf{R}\textbf{x}}$, the iterations in
equation 15 can be employed to optimize the functional in equation 8. To
evaluate step 1 in equation 15, the proximal operator of $f^{*}$ is required:
this can be simply obtained from the proximal operator of the original
function $f$ via the Moreau identity [12], $prox_{\mu
f^{*}}(\mathbf{x})=\mathbf{x}-\frac{1}{\mu}prox_{\frac{f}{\mu}}(\frac{\mathbf{x}}{\mu})$.
All other expression for the proximal operators are the same as those provided
in equations 10 and 11, with the only difference that the inputs in the PD
iterations vary from those of HQS and ADMM. However, a different choice that
can lead to increased reconstruction capabilities is made here: more
specifically, by defining $\textbf{K}=\textbf{S}^{H}$,
$f(\textbf{x})=||\textbf{x}||_{0}$ and
$g(\textbf{x})=i_{\textbf{y}=\textbf{R}\textbf{x}}$, the functional in
equation 8 is again obtained and solved via the Primal-Dual algorithm.
Similarly, the ADMM algorithm can be adapted to such a scenario: the resulting
algorithm takes the name of Linearized-ADMM (L-ADMM) as described in detail in
[12]. Nevertheless, as previously reported in both the context of natural
image reconstruction problems [18] and seismic inversion [19], our numerical
examples will show the superior capabilities of PD against L-ADMM.

Since the transform S is moved out of the function $f$, the Primal-Dual
algorithm can be equally applied to equation 8 with orthogonal and non-
orthogonal sparsifying transforms. This represents a benefit compared to the
HQS formulation. Another benefit of not requiring K to be orthogonal is
represented by the fact that a masking operator $\textbf{M}_{fk}$ operating in
the frequency-wavenumber can be included to enforce the spectrum of the
reconstructed data to lie inside the expected signal cone (i.e., erase energy
from the aliasing regions). We can define
$\textbf{K}=\textbf{M}_{fk}\textbf{S}^{H}$ such that the problem to solve
becomes:

$J_{Masked-
POCS}=||\textbf{M}_{fk}\textbf{S}^{H}\textbf{x}||_{0}\;\text{s.t.}\;\textbf{y}=\textbf{R}\textbf{x}.$
(16)

We note that the idea of including a masking operation within the original
POCS iterations in order to handle regularly sampled data was originally
proposed by [20]. However, since the resulting function $f$ does not have a
closed-form proximal operator, their iterative scheme is not guaranteed to
solve the objective function in equation 16.

## POCS with off-the-grid receivers

In this final section we discuss the implications of off-the-grid receivers
for the POCS algorithm. As briefly mentioned in the Introduction section, the
restriction operator R in equation 1 is not suitable for scenarios where
physical receivers do not lie exactly on on the nominal recording grid. A
bilinear (or sinc) interpolation operator $\textbf{B}\in R^{NN_{t}\times
N_{r,y}N_{r,x}N_{t}}$ is therefore required to handle such a scenario.
Contrary to the restriction operator, this operator does not satisfy the
condition $\textbf{B}\textbf{B}^{H}=\textbf{I}_{NN_{t}}$ required by an Affine
set to have a closed-form solution for its proximal operator (equation 11). As
a consequence, for each proximal evaluation, the problem
$(\textbf{B}\textbf{B}^{H})\textbf{z}=\textbf{B}\textbf{x}-\textbf{y}$ must be
solved. This can be accomplished using a gradient-based iterative solver
(e.g., CGLS, LSQR). Finally, the evaluation of the proximal operator reduces
to $prox_{\frac{1}{\rho}g}(\textbf{x})=\textbf{x}-\textbf{B}^{H}\textbf{z}$.
Although an approximate solution is usually enough and therefore a few
iterations suffice, this slightly increases the cost of the overall
reconstruction process. We note that this differs from what has been
previously suggested in the literature by [21]: their extended POCS (EPOCS)
method ignores the $\textbf{B}\textbf{B}^{H}$ term, which can however deviate
from an identity matrix as later shown in the Numerical Results section.

## Numerical Results

The original POCS algorithm and the proposed modifications are now applied to
a synthetic 3D seismic dataset modeled using a modified version of the
SEG/EAGE Salt and Overthrust Model to which a water column of 300m has been
added on top. The dataset is created using an acoustic finite-difference
modelling code with free-surface activated and a source signature with flat
spectrum between 3 and 45 Hz. The source is placed at a depth of 10m in the
middle of the model and receivers are placed along the seafloor in a regular
grid of size $176\times 80$ with spacing $dx_{R}=20m$ and $dy_{R}=20m$. The
total recording time is 2sec. Figure 1 displays the regularly sampled shot
gather and a decimated version obtained by randomly selecting 40% of the
receiver lines along the X-axis. The reconstruction produced after 80
iterations of the original POCS algorithm with the exponential thresholding
decay proposed in [5] is shown in Figure 2a. The reconstruction obtained from
the Primal-Dual algorithm with constant thresholding is displayed in Figure
2b. In both cases, interpolation is performed in small spatio-temporal
overlapping patches of size $32\times 32\times 32$ with an overlap of $8\times
8\times 6$. Moreover, the seismic shot gather is time-shifted by
$t(y_{R},x_{R})=[(y_{S}-y_{R})^{2}+(x_{S}-x_{R})^{2}+(z_{S}-z_{R})]^{p}/v_{sea}+t_{ZO,off}$,
with $p=0.43$, $v_{sea}=1500m/s$ and $t_{ZO,off}=50ms$, in order to reduce
aliasing prior to reconstruction and shifted back afterwards. Figure 3a
displays the signal-to-noise (SNR) over iterations for the two methods
displayed in Figure 2 as well as for the L-ADMM algorithm (blue line).
Moreover, the final SNRs are reported in the legends of Figure 3. Similarly,
SNR curves are displayed for three other subsampling strategies, namely
dithered subsampling [22] along the X-axis, irregular sampling along both the
X- and Y-axis, and dithered subsampling along both axes. In this figure, we
can observe two different behaviors: when subsampling is performed along the
X-axis, PD both converges faster than HQS and the overall solution is also
improved. On the other hand, when subsampling is performed in a fully random
fashion along the entire grid of receivers, the quality of the seismic data
reconstructed using the new PD algorithm is more in line with that of the
original POCS algorithm, however faster convergence is retained. Note that
whilst creating fully random subsampled data is known to be obtained from a
compressive sensing point of view [23], applying subsampling along a single
spatial direction is much easier to achieve in practical applications: the
clear improvement of our PD-POCS algorithm over the original POCS algorithm is
therefore of significant relevance for the seismic and geophysical industry.

As a second example, we consider the problem of off-the-grid receivers. We use
the same synthetic dataset and create a geometry that contains 60% of the
original receivers (Figure 4a); however, in this case, their position does not
fall exactly on the nominal grid (see Figure 4b). The reconstructed wavefield
using the PD method and 2 iterations of LSQR for the evaluation of the
proximal operator of $g$ is displayed in Figure 4c. The corresponding SNR as a
function of iterations is compared in Figure 4e with that of the EPOCS
algorithm and the PD method with inaccurate evaluation of the proximal of $g$,
i.e. $(\textbf{B}\textbf{B}^{H})^{-1}\approx\mathbf{I}$. We observe that
whilst EPOCS seems to be less affected than PD by such an approximation, the
proposed method outperforms it both in terms of convergence speed and overall
quality of reconstruction. This result highlights the importance of correctly
handling the off-diagonal elements in $(\textbf{B}\textbf{B}^{H})^{-1}$(Figure
4d).

Figure 1: a) Full data, b) Decimated data obtained by selecting 40% of
receiver lines along the X-axis. c-d) Frequency-wavenumber spectra of the full
and decimated data, respectively. Figure 2: Reconstruction with a) POCS
algorithm with exponential decaying threshold, b) PD algorithm with constant
threshold. c-d) Corresponding reconstruction errors. Figure 3: SNRs as
function of iterations for four different subsampling schemes. Figure 4: Off-
the-grid interpolation example. a) Subsampled data, b) Geometry (black dots:
regular grid, green dots: available receivers), c) PD reconstruction, d)
$(\mathbf{BB}^{H})^{-1}$ matrix, and e) SNRs a function of iterations.

## Conclusion

We have presented new insights into the mathematical foundations of the POCS
algorithm based on the theory of proximal operators. The aim of this work is
two-fold: for the math savvy readers, our derivation provides a theoretically
sound explanation for why POCS works in the context on seismic interpolation.
The POCS iterations come essentially from a HQS-based minimization of the
$L_{0}$ norm of the reconstructed seismic data in a transformed domain of
choice, constrained on the actual spatial-temporal measurements. The
connection with the HQS method also explains why better convergence is
observed when employing a varying thresholding as a function of iterations:
this is nothing more than employing the continuation strategy, which is
strictly required by HQS to successfully converge. For the practitioners, we
have suggested a number of changes to the underlying objective function and
optimizer that can be easily implemented in existing codebases; our numerical
example shows how such changes can increase the speed of convergence as well
as the overall quality of the reconstructed wavefield and better handle the
case of off-the-grid receivers.

## Acknowledgements

The authors thank KAUST for supporting this research. All numerical examples
have been produced using the PyLops library [24] and are available at
https://github.com/DIG-Kaust/PyPOCS.

## References

  * [1] Ray Abma and Nurul Kabir.  3d interpolation of irregular data with a POCS algorithm.  Geophysics, 71(6):E91–E97, 11 2006.
  * [2] Pengliang Yang, Jinghuai Gao, and Wenchao Chen.  Curvelet-based POCS interpolation of nonuniformly sampled seismic records.  Journal of Applied Geophysics, 79:90–99, 4 2012.
  * [3] Benfeng Wang, Ru-Shan Wu, Yu Geng, and Xiaohong Chen.  Dreamlet-based interpolation using POCS method.  Journal of Applied Geophysics, 109:256–265, 10 2014.
  * [4] Liu Jicheng, Chou Yongxin, and Zhu Jianjiang.  Interpolating seismic data via the POCS method based on shearlet transform.  Journal of Geophysics and Engineering, 15(3):852–876, mar 7 2018\.
  * [5] Jian-Jun Gao, Xiao-Hong Chen, Jing-Ye Li, Guo-Chang Liu, and Jian Ma.  Irregular seismic data reconstruction based on exponential threshold model of POCS method.  Applied Geophysics, 7(3):229–238, 9 2010.
  * [6] Jianjun Gao, Aaron Stanton, Mostafa Naghizadeh, Mauricio D. Sacchi, and Xiaohong Chen.  Convergence improvement and noise attenuation considerations for beyond alias projection onto convex sets reconstruction.  Geophysical Prospecting, 61:138–151, oct 15 2012.
  * [7] Min Jun Park, Joseph Jennings, Bob Clapp, and Biondo Biondi.  Seismic data interpolation using a POCS-guided deep image prior.  In SEG Technical Program Expanded Abstracts 2020. Society of Exploration Geophysicists, sep 30 2020.
  * [8] R. W. Gerchberg.  A practical algorithm for the determination of phase from image and diffraction plane pictures.  Optik, 35:237–246, 1972.
  * [9] D. Geman and Chengda Yang.  Nonlinear image recovery with half-quadratic regularization.  IEEE Transactions on Image Processing, 4(7):932–946, 7 1995.
  * [10] Ray Abma.  Issues in multidimensional interpolation.  In SEG Technical Program Expanded Abstracts 2009. Society of Exploration Geophysicists, 1 2009.
  * [11] Daniel O’Connor and Lieven Vandenberghe.  Primal-Dual Decomposition by Operator Splitting and Applications to Image Deblurring.  SIAM Journal on Imaging Sciences, 7(3):1724–1754, 1 2014.
  * [12] Neal Parikh.  Proximal Algorithms.  Foundations and Trends in Optimization, 1(3):127–239, 2014.
  * [13] William Menke.  Applications of the POCS inversion method to interpolating topography and other geophysical fields.  Geophysical Research Letters, 18(3):435–438, 3 1991.
  * [14] Zi-Jian Ge, Jing-Ye Li, Shu-Lin Pan, and Xiao-Hong Chen.  A fast-convergence POCS seismic denoising and reconstruction method.  Applied Geophysics, 12(2):169–178, 6 2015.
  * [15] M.M. Nurul Kabir and D.J. Verschuur.  Restoration of missing offsets by parabolic Radon transform1.  Geophysical Prospecting, 43(3):347–368, 4 1995.
  * [16] Stephen Boyd.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.  Foundations and Trends in Machine Learning, 3(1):1–122, 2010.
  * [17] Antonin Chambolle and Thomas Pock.  A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging.  Journal of Mathematical Imaging and Vision, 40(1):120–145, dec 21 2010.
  * [18] Felix Heide, Markus Steinberger, Yun-Ta Tsai, Mushfiqur Rouf, Dawid Pająk, Dikpal Reddy, Orazio Gallo, Jing Liu, Wolfgang Heidrich, Karen Egiazarian, Jan Kautz, and Kari Pulli.  Flexisp.  ACM Transactions on Graphics, 33(6):1–13, nov 19 2014.
  * [19] Matteo Ravasi and Claire Birnie.  A joint inversion-segmentation approach to assisted seismic interpretation.  Geophysical Journal International, 228(2):893–912, sep 23 2021\.
  * [20] Hua Zhang, Hengqi Zhang, Junhu Zhang, Yaju Hao, and Benfeng Wang.  An anti-aliasing POCS interpolation method for regularly undersampled seismic data using curvelet transform.  Journal of Applied Geophysics, 172:103894, 1 2020.
  * [21] Tao Jiang, Bin Gong, Feng Qiao, Yifeng Jiang, Anbo Chen, David Hren, and Zhaobo Meng.  Compressive seismic reconstruction with extended POCS for arbitrary irregular acquisition.  In SEG Technical Program Expanded Abstracts 2017. Society of Exploration Geophysicists, aug 17 2017.
  * [22] Gilles Hennenfent and Felix J. Herrmann.  Simply denoise: Wavefield reconstruction via jittered undersampling.  Geophysics, 73(3):V19–V28, 5 2008.
  * [23] E.J. Candes and M.B. Wakin.  An Introduction To Compressive Sampling.  IEEE Signal Processing Magazine, 25(2):21–30, 3 2008.
  * [24] Matteo Ravasi and Ivan Vasconcelos.  Pylops—A linear-operator Python library for scalable algebra and optimization.  SoftwareX, 11:100361, 1 2020.