SuppReferences

# CLAMP: Majorized Plug-and-Play for Coherent 3D LIDAR Imaging

Tony G. Allen, David J. Rabb,
Gregery T. Buzzard, Charles A. Bouman Tony G. Allen is with the Department of
Mathematics, Purdue University, West Lafayette, IN 47907, USA (email:
allen450@purdue.edu).David J. Rabb is with the Sensors Directorate, Air Force
Research Laboratory, WPAFB, OH 45433, USA (email:
david.rabb@us.af.mil).Gregery T. Buzzard is with the Department of
Mathematics, Purdue University, West Lafayette, IN 47907, USA (email:
<EMAIL_ADDRESS>and was partially supported by NSF CCF-1763896.Charles A.
Bouman is with the School of Electrical and Computer Engineering, Purdue
University, West Lafayette, IN 47907, USA (email<EMAIL_ADDRESS>and was
partially supported by the Showalter Trust.(Corresponding author: Tony G.
Allen)

###### Abstract

Coherent LIDAR uses a chirped laser pulse for 3D imaging of distant targets.
However, existing coherent LIDAR image reconstruction methods do not account
for the system’s aperture, resulting in sub-optimal resolution. Moreover,
these methods use majorization-minimization for computational efficiency, but
do so without a theoretical treatment of convergence.

In this paper, we present Coherent LIDAR Aperture Modeled Plug-and-Play
(CLAMP) for multi-look coherent LIDAR image reconstruction. CLAMP uses multi-
agent consensus equilibrium (a form of PnP) to combine a neural network
denoiser with an accurate physics-based forward model. CLAMP introduces an
FFT-based method to account for the effects of the aperture and uses
majorization of the forward model for computational efficiency. We also
formalize the use of majorization-minimization in consensus optimization
problems and prove convergence to the exact consensus equilibrium solution.
Finally, we apply CLAMP to synthetic and measured data to demonstrate its
effectiveness in producing high-resolution, speckle-free, 3D imagery.

###### Index Terms:

LIDAR, coherent imaging, majorization-minimization, plug-and-play, consensus
equilibrium, deep neural networks.

##  I Introduction

Coherent LIDAR is a method for three-dimensional (3D) imaging in which a
target is flood-illuminated by a frequency-modulated laser source. In multi-
look LIDAR, multiple snapshots of the reflected wave are recorded on a 2D
detector array by a spatial-heterodyne interferometric system [1]. This
process yields a collection of 2D measurements over the duration of the
chirped waveform. A 3D complex-valued volume can be recovered by a
demodulation and filtering process [2, 3, 4]. This volume can be further
processed to produce a 3D image of the surfaces in the scene.

There are several challenges in creating high-quality 3D images from multi-
look LIDAR data. Recent efforts have focused primarily on correcting phase
aberrations caused by non-concurrent measurements [5], object motion [6], or
atmospheric turbulence [7, 8, 9, 10]. However, beyond addressing phase
aberrations, reconstruction algorithms also need to reduce speckle, noise, and
aperture-induced blur.

The most straightforward method to reconstruct multi-look coherent LIDAR
images, often called speckle averaging, involves averaging the squared
magnitude of back-projected, demodulated data from multiple looks. Assuming
the speckle decorrelates between looks, this averaging reduces speckle in the
final reconstruction. However, this method does not incorporate regularization
or priors.

The papers of Pellizzari et al. [11, 12, 13, 14, 15] introduced model-based
iterative reconstruction (MBIR) for 2D coherent LIDAR. This approach uses a
Bayesian framework to estimate the speckle-free, real-valued reflectivity of a
target rather than the speckle-corrupted, complex-valued reflectance. The
real-valued reflectivity roughly corresponds to the image seen using broad-
spectrum illumination; since the reflectivity is relatively smooth, it can be
more strongly regularized to reduce the effective dimensionality of the
estimation problem. While estimating the reflectivity introduces a seemingly
intractable, nonlinear relationship between measured data and the estimated
image, Pellizzari et al. overcomes this problem by using the Expectation-
Maximization (EM) algorithm [16, 17, 18].

An important benefit of this EM-based MBIR approach is that advanced image
priors, such as neural networks, can be incorporated through the use of Plug-
and-Play (PnP) algorithms [19, 20, 21], thereby further improving
reconstruction quality. However, PnP methods critically depend on a high
quality image denoiser for use as the prior model.

In contrast to 2D MBIR problems, which can use 2D image denoisers as prior
models, MBIR reconstruction of 3D LIDAR data requires the formulation of a 3D
prior model that accurately represents the point-cloud of object and material
surfaces. To simplify the creation of such a denoiser, some previous research
assumes there is only a single surface depth at each 2D pixel [22, 23], while
other research uses Markov chain Monte Carlo methods [24] or point cloud
denoisers [25].

Alternatively, Ziabari [26] and Majee [27] introduced the idea of using 2.5D
denoisers for use as PnP priors in 3D X-ray CT reconstruction. This allows
much simpler neural networks to be trained on a small sets of slices within a
3D volume to denoise the center slice. These 2.5D denoisers are then applied
in multiple orientations and combined using the Multi-Agent Consensus
Equilibrium (MACE) framework [28] to provide PnP regularization for a 3D
volume. However, this approach has not been applied to 3D LIDAR image
reconstruction.

Another challenge in LIDAR reconstruction is to efficiently correct the blur
caused by the finite extent of the system’s aperture. The aforementioned MBIR
methods [11, 12, 13, 14, 15] make a simplifying assumption that the forward
model operator is an orthogonal matrix. This greatly reduces computational
complexity, but effectively disregards the system aperture. On a similar
problem in polarimetric synthetic-aperture radar, [29] uses the same EM
algorithm approach, but addresses the aperture problem by using a greedy
graph-coloring-based matrix probing method to do an approximate matrix
inversion, which can be computationally expensive.

Additionally, from a theoretical point of view, a variety of LIDAR
reconstruction algorithms combine the EM algorithm with PnP or ADMM [30];
however, convergence results for this type of combined algorithm are limited.
In [31], the authors prove an ergodic convergence rate when using general
majorization-minimization [32, 33, 34] principles within ADMM, but they do not
prove non-ergodic convergence of the iterates.

In this paper, we present Coherent LIDAR Aperture Modeled Plug-and-Play
(CLAMP), a majorized PnP algorithm for multi-look coherent LIDAR imaging that
accounts for the imaging system’s aperture. This work builds on our previous
publication [35].

Figure 1: Diagram of the CLAMP algorithm. CLAMP uses surrogate-based, multi-
look coherent LIDAR forward agents to fit to data and account for the system
aperture. These are incorporated with multi-slice fusion 2.5D neural network
denoisers in a consensus equilibrium framework to produce a 3D reconstruction
with reduced blur, noise, and speckle.

Figure 1 graphically illustrates the essential concept behind CLAMP in which
the information from each of $L$ “looks” of a LIDAR sensor are integrated
along with a neural network prior model of the 3D scene to form a 3D
reconstruction that has reduced noise, speckle and increased resolution
relative to traditional reconstruction methods. This integration is done using
the MACE algorithm, an extension of PnP, that provides a computationally
efficient method for integrating neural network prior models with physics-
based forward models. Specifically, we make the following novel contributions:

  * •

Accurate LIDAR aperture modeling using a computationally efficient FFT-based
method;

  * •

Incorporation of a 2.5D convolutional neural network prior in 3D multi-look
coherent LIDAR reconstruction;

  * •

Formulation of a unified framework for majorization-minimization within MACE
with theoretical guarantees of convergence to an exact MACE solution;

  * •

Results on simulated and measured data that show CLAMP yields 3D images with
reduced speckle and improved resolution.

A key component of CLAMP is the use of majorization-minimization, or surrogate
function optimization, in the MACE framework. We use convex surrogate
functions in place of the proximal agents in MACE and show that this method
converges to the exact MACE solution. This surrogate-based MACE applies to
many imaging problems and modalities and encompasses several common
approximation methods and existing algorithms as special cases.

##  II Background

In this section, we provide an overview of the multi-look coherent LIDAR
imaging problem, develop the models used in our approach, and introduce the
MACE framework.

###  II-A Coherent LIDAR Imaging Measurement System

Figure 2: The measurement process and conventional data processing pipeline
for a 3D coherent LIDAR imaging system using image plane recording geometry.
(a) Multiple chirped laser pulses, or looks, illuminate a target. Each look is
assumed to generate an independent speckle pattern, $g_{\ell}$. Reflected
light propagates to an optical system with aperture, $a$, where a lens focuses
an image onto a focal plane array. The image is interfered with a delayed
local reference beam to form a hologram, which is sampled in time to create
$N_{t}$ frames for each look. (b) The 2D DFT of each hologram frame is taken
to separate the frequency components of the target and reference beam. The
pupil images are then windowed, centered, and zero-padded to form a 3D stack
of complex fields, $y_{\ell}$ for the ${\ell}^{th}$ look. (c) A noisy 3D image
of the target can be reconstructed from $y_{\ell}$ by the inverse 3D DFT and
averaging over $\ell=1,\ldots,L$.

Figure 2(a) illustrates a typical coherent LIDAR imaging system with an image-
plane recording geometry. In this system, a series of $L$ chirped laser
pulses, or looks, illuminates a target with reflectivity, $r$. We assume that
the target is stationary and that the speckle realizations and measurement
noise samples are statistically independent from look to look. For each look,
$\ell=1,\ldots,L$, the corresponding complex reflectance, $g_{\ell}$
propagates to the aperture, $a$, where a thin lens forms an image of the
target on a focal plane array. The resulting image is mixed with a delayed,
off-axis, local reference beam, which results in an interference pattern, or
hologram, on the focal plane array.

During a single look, a series of $N_{t}$ hologram images or frames are
captured by the sensor in time. Intuitively, each of these 2D holograms
corresponds to a slightly different laser wavelength due to the chirp slope.
This stack of 2D holograms will then be used to recover depth information from
the scene [2, 4, 3].

Figure 2(b) and (c) illustrate the conventional processing steps to
reconstruct a 3D image of the target from this stack of 2D holograms and shows
how the 3D image can be decomposed for visualization purposes into a depth
image and a reflectivity image at the corresponding depths. First, the 2D DFT
is taken for each frame in a look. Since the frames are real valued, each DFT
results in a complex image that contains two circular disks corresponding to
the complex conjugate circular apertures in the pupil plane. One of the two
circular apertures is then embedded in a bounding square array (shown with a
red dotted line) and windowed to give the pupil image. This image is zero-
padded to increase the effective sampling rate in the space and time domains
and yields a larger 3D stack of pupil images, $y_{\ell}$. Specifically, we
zero-pad the pupil image by a factor of $q$ in each of the three dimensions;
so the final array, $y_{\ell}$, is $q^{3}$ times the size of the un-padded
images. The last step in conventional reconstruction is to take the inverse 3D
DFT of each $y_{\ell}$ and average over $\ell$, which yields a candidate
reconstruction in cross-range and depth.

The final step in conventional reconstruction limits reconstruction quality in
two key ways: lost spatial and temporal resolution due to the missing
information in the frequency domain, and noise due to speckle.

In contrast, our algorithm replaces the conventional final 3D DFT and
averaging step with an iterative model-based strategy which employs a physics-
based forward model in conjunction with neural network denoisers to ensure the
reconstructed image is high-resolution and free of speckle and noise.

###  II-B Forward Model

In order to derive a forward model for our coherent LIDAR system, we first
relate the complex reflectivity from the target, $g_{\ell}\in\mathbb{C}^{n}$,
to the pre-processed measurements, $y_{\ell}\in\mathbb{C}^{n}$, from the
previous section. Using standard assumptions [36], we have that

$y_{\ell}=Ag_{\ell}+w_{\ell},$ (1)

where $A$ is a linear operator that models the imaging system, and
$w_{\ell}\sim\text{CN}(0,\sigma_{w}^{2}I)$ is circularly-symmetric complex
Gaussian white noise, where $z\sim\text{CN}(\mu,\Gamma)$ denotes the
$n$-dimensional complex density

$p(z)=\frac{1}{\pi^{n}\mbox{det}(\Gamma)}\exp(-(z-\mu)^{H}\Gamma^{-1}(z-\mu)).$

Assuming that the target subtends a small angle and the chirp exhibits a
narrow fractional bandwidth, the matrix $A$ can be approximated as

$A=\mathcal{D}(a)F,$ (2)

where $F$ is the orthonormal 3D DFT, $\mathcal{D}(\cdot)$ represents a
diagonal matrix with entries given by its argument, and the vector $a$ is a
binary vector that encodes the aperture of the imaging system after zero-
padding.

Our ultimate goal will be to recover $r\in\mathbb{R}^{n}$, the real-valued
speckle-free reflectivity at each 3D voxel. While the data $y_{\ell}$ directly
relates to $g_{\ell}$ by (1), images reconstructed from estimates of
$g_{\ell}$ are degraded by speckle. In order to address this, we need a
forward model that relates $r$ to $g_{\ell}$.

We will assume that each look produces a sample of $g_{\ell}$ that is
conditionally independent given $r$. Using a standard model of fully developed
speckle [37, 38], we can model the conditional distribution of $g_{\ell}$
given $r$ as being a circularly-symmetric complex Gaussian random variable
with conditional distribution

$g_{\ell}|r\sim\text{CN}(0,\mathcal{D}(r))\ ,$ (3)

where $\mathcal{D}(r)$ is a diagonal covariance with entries $r$. We can then
directly model the data $y_{\ell}$ in terms of $r$ by composing the forward
model (1) with the speckle model of (3) to yield

$y_{\ell}|r\sim\text{CN}(0,A\mathcal{D}(r)A^{H}+\sigma_{w}^{2}I)\ ,$ (4)

where $y_{\ell}$ for $\ell=1,\dots,L$ are assumed conditionally independent
given $r$.

###  II-C MACE Formulation

Our goal is to reconstruct $r$ from a collection of $L$ coherent LIDAR looks,
$\\{y_{\ell}\\}_{\ell=1}^{L}$. In order to do this, we will use the MACE
framework in which multiple agents, such as data-fitting operators or
denoisers, are balanced to form a single, consistent reconstruction [28]. We
give an overview of MACE here and provide details of the agents used in CLAMP
in Section III.

For our problem with $L$ looks of data, we will have $L$ data-fitting agents,
$F_{\ell}(w_{\ell})$ for $\ell=1,\ldots,L$. Specifically, $F_{\ell}$ is a
proximal map given by

$F_{\ell}(w_{{\ell}})=\underset{r}{\operatorname{argmin}}\left\\{f_{\ell}(r)+\frac{1}{2\sigma^{2}}\left\lVert
r-w_{\ell}\right\rVert^{2}\right\\}.$ (5)

where $f_{\ell}(r)=-\log p\left(y_{\ell}|r\right)$ is the negative log
likelihood associated with the ${\ell}^{th}$ look. These operators have an
intuitive interpretation of taking a full reconstruction $w_{\ell}$ as input
and returning an image close to $w_{\ell}$ that better fits the data from look
$\ell$.

We will also use three denoising agents, $H_{1}$, $H_{2}$, and $H_{3}$, which
are 2.5D denoisers [26] applied respectively along $xy$, $xz$, and $yz$ slices
of the 3D reconstruction. Each denoising agent takes a full 3D reconstruction
as input and generates a denoised 3D reconstruction as output; more details
are given in Section III-B.

To formulate the MACE solution, we concatenate these $L+3$ agents into a
single operator $\mathbf{F}$ given by

$\mathbf{F}(\mathbf{w})\\!=\\!\left[F_{1}(w_{1}),\ldots,F_{L}(w_{L}),H_{1}(w_{L+1}),\cdots,H_{3}(w_{L+3})\right],$
(6)

where $\mathbf{w}=\left[w_{1},\ldots,w_{L+3}\right]$ is the stacked input. We
also define the operator

$\mathbf{G}(\mathbf{w})=\left[\overline{\mathbf{w}},\ldots,\overline{\mathbf{w}}\right]$
(7)

where
$\overline{\mathbf{w}}=\frac{1}{2}\left[\frac{1}{L}\sum_{\ell=1}^{L}w_{\ell}+\frac{1}{3}\sum_{i=L+1}^{L+3}w_{i}\right]$
is a weighted average of the components, with equal weight given to the data
agents and the denoising agents.

The MACE reconstruction is then determined by the solution to the equilibrium
equation

$\mathbf{F}(\mathbf{w}^{*})=\mathbf{G}(\mathbf{w}^{*})\ .$ (8)

The final reconstruction, $r^{*}$, can be obtained from any single component
of $\mathbf{w}^{*}$, however, we typically use
$r^{*}=\overline{\mathbf{w}}^{*}$ for numerical robustness.

The MACE equation of (8) is commonly solved as a fixed point of the operator
$\mathbf{T}=\left(2\mathbf{G}-\mathbf{I}\right)\left(2\mathbf{F}-\mathbf{I}\right)$
by the iterations

$\mathbf{w}\leftarrow(1-\rho)\mathbf{w}+\rho\mathbf{T}\mathbf{w},$ (9)

where $\rho\in(0,1)$. This approach is summarized in Algorithm 1.

Algorithm 1 MACE

1:Input: Initialize $\mathbf{w}$

2:while not converged do

3: $\mathbf{r}\leftarrow\mathbf{F}(\mathbf{w})$

4: $\mathbf{x}\leftarrow 2\mathbf{r}-\mathbf{w}$

5:
$\mathbf{w}\leftarrow\mathbf{w}+2\rho\left(\mathbf{G}(\mathbf{x})-\mathbf{r}\right)$

6:end while

7:Output: $r^{*}=\overline{\textbf{w}}$

##  III CLAMP

In this section, we present the Coherent LIDAR Aperture Modeled Plug-and-Play
(CLAMP) algorithm for reconstructing 3D images from multi-look coherent LIDAR
measurements. CLAMP is built on the MACE framework using $L+3$ agents. The
first $L$ agents are forward agents that enforce fidelity to the measurements,
$y_{\ell}$, for $\ell=1,\ldots,L$, and the last three agents are 2.5D CNN
denoisers that regularize the image in all three dimensions. The consensus
equilibrium of these agents yields an image that is both consistent with all
measurements and regularized in all three dimensions. The details of these
agents are given in the following sections, and the full algorithm is given in
Algorithm 2.

###  III-A EM Surrogate Forward Agents

Unfortunately, direct application of the MACE iterations in (9) is
computationally infeasible. This is because the forward model in (4) quickly
becomes intractable for even moderately large images due to the nonlinear
relationship between $r$ and $y_{\ell}$.

In order to overcome this problem, we will take the approach first proposed by
Pellizzari, et. al [11, 12, 13, 14, 15] and use the EM algorithm to compute
surrogate functions to the exact negative log likelihood functions. These
surrogate functions are given by

$\hat{f}_{\ell}(r;r_{\ell}^{\prime})=\mathbb{E}_{g|y_{\ell},r_{\ell}^{\prime}}\left[-\log
p\left(y_{\ell},g\middle|r\right)\right],$ (10)

where the expectation is taken over the conditional distribution of $g$ given
$y_{\ell}$ and $r_{\ell}^{\prime}$, an approximation of $r$. Assuming (1) and
(3), and using $y_{\ell}|g,r=y_{\ell}|g$, [11] shows that this conditional
distribution of $g$ given $y_{\ell},r_{\ell}^{\prime}$ is Gaussian with the
form

$\displaystyle
p(g|y_{\ell},r_{\ell}^{\prime})=\frac{1}{Z}\exp\left\\{-(g-\mu_{\ell})^{H}\tilde{C}_{\ell}^{-1}(g-\mu_{\ell})\right\\},$
(11)

where $Z$ is a normalizing constant and $\mu_{\ell}$ and $\tilde{C}_{\ell}$
are the mean and covariance given by

$\displaystyle\mu_{\ell}$
$\displaystyle=\underset{g\in\mathbb{C}^{n}}{\operatorname{argmin}}\left\\{\frac{1}{2\sigma_{w}^{2}}\lVert
y_{\ell}-Ag\rVert^{2}+\frac{1}{2}g^{H}\mathcal{D}\left(\frac{1}{r_{\ell}^{\prime}}\right)g\right\\},$
(12) $\displaystyle\tilde{C}_{\ell}$
$\displaystyle=\left[\frac{1}{\sigma_{w}^{2}}A^{H}A+\mathcal{D}\left(\frac{1}{r_{\ell}^{\prime}}\right)\right]^{-1}.$
(13)

However, using the surrogate function in (10) is also intractable for large
images since computing the covariance, $\tilde{C}_{\ell}$, requires a large
matrix inversion, while computing the mean, $\mu_{\ell}$, of (12) directly
requires the solution of a large quadratic optimization problem. We address
each of these issues separately.

First, to improve the tractability of computing $\tilde{C}_{\ell}$, we make
the approximation that $A^{H}A=\alpha I$, where $\alpha=\lVert a\rVert_{1}/n$
has the interpretation of the amount of light passing through the aperture. In
this case, $\tilde{C}_{\ell}$ is replaced by diagonal matrix $C_{\ell}$ with
diagonal entries $c_{\ell,j}$, given by

$c_{\ell,j}=\frac{\sigma_{w}^{2}r_{\ell,j}^{\prime}}{\alpha
r_{\ell,j}+\sigma_{w}^{2}}.$ (14)

Next, to evaluate (12) efficiently and accurately, we incorporate a gradient
step in each iteration of the EM algorithm. To avoid division by zero in the
objective in (12), we add the noise floor $\sigma_{w}^{2}/\alpha$ to
$r_{\ell}^{\prime}$. This yields the regularized objective

$h(g;y_{\ell},r_{\ell}^{\prime})=\frac{1}{2\sigma_{w}^{2}}\lVert
y_{\ell}-Ag\rVert^{2}+\frac{1}{2}g^{H}\mathcal{D}\left(\frac{1}{r_{\ell}^{\prime}+\sigma_{w}^{2}/\alpha}\right)g.$
(15)

The gradient step is then given by

$\displaystyle d$ $\displaystyle\leftarrow-\nabla
h(\mu_{\ell};y_{\ell},r_{\ell}^{\prime})$ (16) $\displaystyle\mu_{\ell}$
$\displaystyle\leftarrow\mu_{\ell}+\gamma_{*}d,$ (17)

where the optimal step-size, $\gamma_{*}$, is computed using a line search,
which has a closed-form solution [21, Chapter 5.3].

Finally, using these forms for $\mu_{\ell,j}$ and $c_{\ell,j}$, the EM
surrogate of (10) can be calculated [11] to be

$\displaystyle\hat{f}_{\ell}(r;\mu_{\ell},c_{\ell})$
$\displaystyle=\sum_{j=1}^{n}\left\\{\log
r_{j}+\frac{|\mu_{\ell,j}|^{2}+c_{\ell,j}}{r_{j}}\right\\}.$ (18)

From this, the CLAMP forward agent is the proximal map

$\hat{F_{\ell}}(v;\mu_{\ell},c_{\ell})=\underset{r}{\operatorname{argmin}}\left\\{\hat{f_{\ell}}(r;\mu_{\ell},c_{\ell})+\frac{1}{2\sigma^{2}}\left\lVert
r-v\right\rVert^{2}\right\\},$ (19)

which can be computed in closed form by solving for a root of a cubic equation
for each voxel in the image. However, in practice, to further simplify
computation and speed up implementation, we employ an additional quadratic
surrogate of $\hat{f}_{\ell}$. The details of which are given in Appendix A.

###  III-B 2.5D CNN Prior Agents

For the CLAMP prior agents, we take an approach similar to Multi-Slice Fusion
[27], in which three prior agents, $H_{1},H_{2},H_{3},$ are used. These agents
are denoted together as H in Algorithm 2.

These three prior agents are 2.5D CNNs, each of which denoises 2D slices
perpendicular to one of the 3 coordinate axes. Each of these networks takes as
input a block of 17 slices of a 3D image and then outputs the denoised center
slice. One prior agent is obtained from such a network by choosing a
coordinate axis and applying the network to denoise each slice perpendicular
to this coordinate axis. MACE enforces a consensus among the three prior
agents, resulting in an image that is regularized in all three dimensions.
With this approach, local 3D information is leveraged to denoise the image
while reducing the memory requirements of the network.

We trained a 2.5D CNN to denoise 3D images from ShapeNet [39] with 10%
additive white Gaussian noise and then applied this network in 3 orientations
to form the 3 agents. The details of the network architecture and training are
given in Appendix B.

Algorithm 2 CLAMP

1:Input: $y_{\ell}$ for $\ell=1,\ldots,L$

2:Initialize: $\mathbf{w}$, $\mathbf{r}$, $\mu$

3:while not converged do

4: for $\ell=1,\dots,L$ do c

5: $\forall j$, $c_{\ell,j}\leftarrow\frac{\sigma_{w}^{2}r_{\ell,j}}{\alpha
r_{\ell,j}+\sigma_{w}^{2}}$ $\triangleright$ Update $C_{\ell}$

6: $d\leftarrow-\nabla h(\mu_{\ell};y_{\ell},r_{\ell})$ $\triangleright$
Update $\mu_{\ell}$

7: $\gamma\leftarrow\operatorname*{argmin}_{\gamma}\left\\{h(\mu_{\ell}-\gamma
d;y_{\ell},r_{\ell})\right\\}$

8: $\mu_{\ell}\leftarrow\mu_{\ell}+\gamma d$

9: $r_{\ell}\leftarrow\hat{F}_{\ell}\left(w_{\ell};\mu_{\ell},c_{\ell}\right)$
$\triangleright$ Update $r_{\ell}$ by (19)

10: end for

11: $r_{L+1},r_{L+2},r_{L+3}$ $\triangleright$ Apply prior agents

12: $\leftarrow\textbf{H}(w_{L+1},w_{L+2},w_{L+3})$

13: $\mathbf{r}\leftarrow\left[r_{1},\ldots,r_{L+3}\right]$ $\triangleright$
Do MACE iteration

14: $\mathbf{x}\leftarrow 2\mathbf{r}-\mathbf{w}$

15:
$\mathbf{w}\leftarrow\mathbf{w}+2\rho\left(\mathbf{G}(\mathbf{x})-\mathbf{r}\right)$

16:end while

17:Output: $r^{*}=\overline{\textbf{w}}$

##  IV Majorized MACE Theory

A key component of CLAMP is the use of a surrogate function through the EM
algorithm. In this section, we present a generalized theory of majorization-
minimization within the MACE framework with guaranteed convergence to an exact
MACE solution as defined in (8).

Although our theory is limited to the case where all agents are proximal maps,
it applies in practice to many widely used denoisers. Many PnP algorithms
(including CLAMP) use denoising agents that are not necessarily proximal maps.
However, typical CNN denoisers are trained to be minimum mean squared error
(MMSE) denoisers, which were shown to be proximal maps in [40]. This
equivalence was also demonstrated experimentally in [41], in which the authors
show that the convergence of a PnP algorithm with an exact MMSE denoiser
agrees remarkably well with the convergence of the same algorithm with a CNN
denoiser. This highlights the practical relevance of our theory.

In our setting, we modify the operator $F_{i}$ by approximating its objective
function, $f_{i}$, with a surrogate function, $\hat{f_{i}}$, whose proximal
map, $\hat{F_{i}}$, can be computed more efficiently. We make a few standard
assumptions on the surrogate functions, given in Definition 1, which are
adapted from those made in [33, 42, 43]. Commonly used surrogates, such as the
Jensen surrogate used in the EM algorithm or quadratic approximations of twice
differentiable functions meet these requirements.

######  Definition 1 (Surrogate).

Let $\Theta\subset\mathbb{R}^{n}$ be convex, and
$f:\mathbb{R}^{n}\to\mathbb{R}$ be a convex function. A function
$\hat{f}:\mathbb{R}^{n}\to\mathbb{R}$ is a surrogate of $f$ near
$\xi\in\Theta$ when the following conditions hold:

  * •

Majorization: we have $\hat{f}(x)\geq f(x)$ for all $x\in\mathbb{R}^{n}$.

  * •

Smoothness: the approximation error $e\triangleq\hat{f}-f$ is differentiable,
and its gradient is $L$-Lipschitz continuous. Moreover, we have that
$e(\xi)=0$ and $\nabla e(\xi)=0$,

  * •

Strong convexity: $\hat{f}$ is $p$-strongly convex with $p\geq L$.

When the point $\xi$ is relevant, we write the surrogate as $\hat{f}(x;\xi)$,
and we denote the set of such surrogates as $\mathcal{S}_{L,p}(f,\xi)$.

We do not require the objective function, or its surrogate, to be
differentiable, nor do we require the surrogate function to be continuous as a
function of $\xi$. These are common assumptions in the literature, but are not
necessary for our theory.

Majorization-minimization schemes work by alternately minimizing the surrogate
function and then updating the surrogate itself, which gives updates of the
form

$r^{(k+1)}\leftarrow\underset{x}{\operatorname{argmin}}\left\\{\hat{f}\left(r;r^{(k)}\right)\right\\}.$
(20)

The conditions in Definition 1 ensure that the original objective is
decreasing, $f\left(r^{(k+1)}\right)\leq f\left(r^{(k)}\right)$, and that the
iterates $r^{(k)}$ converge to a minimizer of $f$.

Algorithm 3 applies these principles within the MACE framework. Assuming we
have $N$ agents, on the $k$-th iteration, we compute a surrogate $\hat{f_{i}}$
of the objective function $f_{i}$ at $r_{i}^{(k)}$ for $i=1,\dots,N$. We
define the concatenated operators

$\hat{\mathbf{F}}^{(k)}\left(\mathbf{w};\mathbf{r}^{(k)}\right)=\left[\hat{F}_{1}^{(k)}(w_{1};r_{1}^{(k)}),\dots,\hat{F}_{N}^{(k)}(w_{N};r_{N}^{(k)})\right]$
(21)

and

$\hat{\mathbf{T}}^{(k)}=\left(2\mathbf{G}-I\right)\left(2\hat{\mathbf{F}}^{(k)}\left(\cdot\,;\mathbf{r}^{(k)}\right)-I\right),$
(22)

where $\hat{F}_{i}^{(k)}$ is the proximal map of $\hat{f}_{i}^{(k)}$. Akin to
the fixed point iterations in (9), we analyze the fixed point of the system

$\displaystyle\mathbf{w}^{(k+1)}$
$\displaystyle=\rho\hat{\mathbf{T}}^{(k)}\mathbf{w}^{(k)}+(1-\rho)\mathbf{w}^{(k)},$
(23) $\displaystyle\mathbf{r}^{(k+1)}$
$\displaystyle=\hat{\mathbf{F}}^{(k)}(\mathbf{w}^{(k)};\mathbf{r}^{(k)}).$

Our main result is Theorem 1, which shows that Algorithm 3 converges to an
exact MACE solution. Part (i), states that any fixed point
$(\mathbf{w}^{*},\mathbf{r}^{*})$ of (23) is an exact MACE solution as in (8).
Part (ii) states that convergence to a fixed point is guaranteed if the
surrogate is strongly convex. Note that if
$\mathbf{r}^{(k+1)}=\mathbf{r}^{(k)}$, then
$\hat{\mathbf{F}}^{(k+1)}=\hat{\mathbf{F}}^{(k)}$ and
$\hat{\mathbf{T}}^{(k+1)}=\hat{\mathbf{T}}^{(k)}$, so the functions are also
fixed at such a fixed point.

######  Theorem 1 (Majorized-MACE Solution).

Let $f_{i}:\mathbb{R}^{n}\to\mathbb{R},i=1,\dots,N$, be convex, and
$\hat{f_{i}}^{(k)}\in\mathcal{S}_{L_{i},p_{i}}\left(f_{i},r_{i}^{(k)}\right)$.
Let $F_{i}$ and $\hat{F_{i}}^{(k)}$ be the proximal maps of $f_{i}$ and
$\hat{f_{i}}^{(k)}$, respectively, and $\hat{\mathbf{T}}^{(k)}$ be the
operator defined in (22). Then:

  1. (i)

Any fixed point of the augmented system (23) is a solution to the original
MACE equation given in (8).

  2. (ii)

The iterates defined by (23) converge to a fixed point, and hence a MACE
solution, if one exists.

###### Proof.

Proof is in Appendix C. ∎

##  V Results

We evaluate the performance of CLAMP using synthetic data (Section V-D) and
experimental data (Section V-E), and we compare CLAMP reconstructions with the
standard back-projection method, or speckle averaged image,
$\frac{1}{L}\sum_{\ell=1}^{L}\lvert A^{H}y_{\ell}\rvert^{2}$. To highlight the
importance of our aperture model, we also compare to a version of CLAMP with
the aperture modeled as $a=\mathbf{1}$. In this case, we use the same
algorithm as described in Algorithm 2, but with the exact solutions of
$\mu_{\ell}$ and $C_{\ell}$.

Algorithm 3 Majorized-MACE

1:Input: Initialize $\mathbf{w}^{(0)}$, $\mathbf{r}^{(0)}$

2:$k=0$

3:while not converged do

4: for $i=1,\dots,N$ do

5: Compute surrogate
$\hat{f_{i}}^{(k)}\in\mathcal{S}_{L_{i},p_{i}}\left(f_{i},r_{i}^{(k)}\right)$

6: end for

7:
$\mathbf{r}^{(k+1)}=\hat{\mathbf{F}}^{(k)}(\mathbf{w}^{(k)};\mathbf{r}^{(k)})$

8: $\mathbf{x}=2\mathbf{r}^{(k+1)}-\mathbf{w}^{(k)}$

9:
$\mathbf{w}^{(k+1)}=\mathbf{w}^{(k)}+2\rho\left(\mathbf{G}(\mathbf{x})-\mathbf{r}^{(k+1)}\right)$

10: $k=k+1$

11:end while

###  V-A Methods

For each CLAMP reconstruction in this section, we ran the corresponding
algorithm for 250 iterations with $\rho=0.5$. Each algorithm was initialized
with $r=\frac{1}{L}\sum_{\ell=1}^{L}\lvert A^{H}y_{\ell}\rvert^{2}$ and
setting $\mathbf{w}$, $\mathbf{r}$, to be stacked copies of $r$. Each
$\mu_{\ell}$ was initialized as $\frac{1}{\alpha}A^{H}y_{i}$. The average
$\overline{\mathbf{w}}$ was computed so that the $L$ forward agents equally
shared half of the weight and the 3 prior agents equally shared the other
half. The forward model proximal parameter, $\sigma^{2}$, was empirically
chosen in a range of $[0.001,1.0]$ to produce the best image quality.

To measure convergence of CLAMP, we define the convergence error based on the
MACE equation (8) as

$\text{Convergence
Error}=\frac{\left\lVert\hat{\mathbf{F}}(\mathbf{w})-\mathbf{G}(\mathbf{w})\right\rVert}{\left\lVert\mathbf{G}(\mathbf{w})\right\rVert}.$
(24)

Additionally, to show convergence of the $\mu_{\ell}$ updates, we define the
relative residual averaged across looks as

$\mu\text{-Residual}=\frac{1}{L}\sum_{\ell=1}^{L}\frac{\left\lVert\tilde{C}_{\ell}^{-1}\mu_{\ell}-\frac{1}{\sigma_{w}^{2}}A^{H}y_{\ell}\right\rVert}{\left\lVert\frac{1}{\sigma_{w}^{2}}A^{H}y_{\ell}\right\rVert},$
(25)

where the numerator comes from the exact solution of (12). Plots of CLAMP’s
convergence behavior are shown in Section V-E.

In the synthetic data experiments, we compare the reconstructed images to a
ground truth image. Since the ground truth image is a different size than the
reconstructed images, we cannot compare them without interpolation. Instead,
we convert the images to point clouds. We convert a 3D image, $r$, to a point
cloud, $P\subseteq\mathbb{R}^{4}$, such that each point
$p=(\textbf{s}_{p},r_{p})\in P$ corresponds to a voxel in $r$ with a value
larger than a given threshold. In this work, we use the natural threshold, the
noise floor, $\sigma_{w}^{2}/\alpha$. The first component,
$\textbf{s}_{p}\in\mathbb{R}^{3}$, of $p$ represents the spatial coordinates
of the point, and $r_{p}$ represents that point’s reflectivity.

We quantify the spatial accuracy of the reconstructions by computing the
Euclidean distance between points in the reconstruction and points in the
ground truth point cloud. For each reconstruction, $\hat{r}$, we transform the
image into a point cloud $P=\left\\{(\textbf{s}_{p},\hat{r}_{p})\right\\}$.
Averaging over each point in $P$, we compute the Euclidean distance to its
closest point in the ground truth point cloud,
$Q=\left\\{(\textbf{t}_{q},r_{q})\right\\}$. Let
$\tilde{n}(p)=\operatorname*{argmin}_{q\in
Q}||\textbf{s}_{p}-\textbf{t}_{q}||$ be the point in $Q$ nearest to $p$,
which, in practice, we determine by a $k$-d tree search [44]. The average
Euclidean distance is then given by

$\text{Euclidean Distance}=\frac{1}{|P|}\sum_{p\in
P}||\textbf{s}_{p}-\textbf{t}_{\tilde{n}(p)}||.$ (26)

To avoid corrupting this metric with outliers (points in the reconstruction
with no nearby neighbor in the ground truth), we remove any point with a
distance to the ground truth more then 3 times the Rayleigh criterion for
resolution ($\qty{1.5}{}$\mathrm{c}\mathrm{m}$$).

Using this method, we can also compute the normalized root mean squared error
(NRMSE) of a reconstruction as

$\text{NRMSE}(\alpha^{*}\hat{r},r)=\sqrt{\frac{\sum_{p\in
P}(\alpha^{*}\hat{r}_{p}-r_{\tilde{n}(p)})^{2}}{\sum_{p\in
P}|r_{\tilde{n}(p)}|^{2}}},$ (27)

where $\alpha^{*}$ is computed as

$\alpha^{*}=\underset{\alpha}{\operatorname{argmin}}\left\\{\text{NRMSE}(\alpha\hat{r},r)\right\\}.$

This multiplicative factor accounts for any scaling differences in the
reconstructed images and allows for a fair comparison between the different
reconstructions.

Finally, to further quantify the improvement of resolution brought by CLAMP,
we use the Fourier Shell Correlation (FSC) [45] function. The FSC function
provides a measure of similarity between the reconstructed image and the
ground truth image over spherical shells in Fourier space. Generally, higher
resolution reconstructions will show larger correlations at higher
frequencies. To make this precise, an effective resolution can be determined
by the intersection of the FSC curve with a threshold curve [46].

Figure 3: Ground truth target used for synthetic data generation. The 3D image
is illustrated as a pair of depth (upper-left) and reflectivity (lower-right)
images. The reflectivity image is computed as the maximum reflectivity value
along the depth dimension, and the range image is the location of such
maximum.

###  V-B Synthetic Data Generation

In the synthetic data experiment, we generated data by simulating the multi-
look coherent LIDAR imaging process described in Section II-B. We repeated
this experiment ten times to ensure statistical robustness of our results.

Each data simulation included nine looks at a 3D model of a train, shown in
Figure 3 as a pair of depth and reflectivity images. The reflectivity image is
given by the maximum value along the depth dimension, while the depth is given
by the location of the maximum.

The train was assumed to be 1 meter in size in each dimension. Assuming the
train was a lambertian surface, the reflectivity of each point on the target
was computed as proportional to the cosine of the angle between the surface
normal and the viewing direction. The reflectivity was then normalized so that
the brightest point in the image had unit reflectivity. A speckle realization
was generated by multiplying the reflectivity of each voxel in the image by a
complex Gaussian random variable with unit variance. The speckle realization
was then propagated to the hologram plane by applying the forward model with
parameters listed in the ‘Simulation’ column of Table I. The aperture was
modeled as a centered, circular aperture with a diameter equal to $50\%$ of
the hologram grid length. Finally, the data is corrupted by adding complex
white Gaussian noise with variance $\sigma_{w}^{2}=10^{-3}$.

###  V-C Experimental Data Measurement

In addition to synthetic data, we evaluate CLAMP on data collected at the Air
Force Research Laboratory. The experiment consisted of two targets, a toy car
and a hexagonal nut, which were painted with matte white paint so that the
surfaces had a uniform Lambertian reflectance. In order to measure multiple,
statistically independent speckle realizations, the targets were placed on a
high-precision rotation stage and rotated slightly after each measurement.

The targets were illuminated by a linear-frequency modulated waveform with
central wavelength of 1550 (193.4THz) and chirp rate of 117.9/s. The laser was
split so that 95% of the power was transmitted to the target and the remaining
5% was used as a reference beam for holographic imaging. The reference beam
was placed off-axis in the plane of the exit pupil and pointed at the center
of a focal plane array, where it interfered with a focused image of the
target. The interference pattern was recorded on a $640\times 512$ pixel
InGaAs photodetector with 20μm pixel pitch at 8-bit resolution. The camera’s
region of interest was narrowed to fit the image of the target, thereby
allowing faster frame rates. The toy car was measured by a $128\times 128$
array of pixels at a frame rate of 17.60, and the hexagonal nut was measured
by an $80\times 80$ array at a frame rate of 28.92. The noise variance was
estimated to be $\sigma_{w}^{2}=0.0012$ for the toy car and
$\sigma_{w}^{2}=0.0025$ for the hexagonal nut. This estimation was done by
computing the mean power within the pupil region of the hologram and dividing
by the mean power outside of the pupil region. This is similar to the noise
variance used in the synthetic data experiment, which allows for a comparison
of the reconstruction quality between the two experiments.

Further details of the experimental setup, system hardware, and calibration
methods can be found in [7, Chapter 4]. We summarize important experimental
parameters in Table I.

TABLE I: Parameters used in our multi-look coherent LIDAR simulation and experimental measurements of a hexagonal nut and toy car. Item | Simulation | Hex. nut | Toy car | Units
---|---|---|---|---
Distance | 52.9 |  2.64 |  2.64 | $\mathrm{m}$
Central wavelength |  1550 |  1550 |  1550 | $\mathrm{n}\mathrm{m}$
Chirp rate | — |  117.9 |  117.9 |  $\mathrm{THz}$/$\mathrm{s}$
Chirp duration | — |  2.0 |  2.0 | $\mathrm{ms}$
Frame rate | — |  28.9 |  17.6 | $\mathrm{kHz}$
Frames, $N_{t}$ | $64$ | $46$ | $28$ | samples
Frequency step size |  0.15 |  6.7 |  4.2 | $\mathrm{GHz}$
Aperture diameter |  6.4 |  6.4 |  6.4 | $\mathrm{m}\mathrm{m}$
Focal length | — |  18.2 |  18.2 | $\mathrm{c}\mathrm{m}$
Pixel pitch | — |  20.0 |  20.0 |  $\mu$m
Hologram grid size | $(128,128)$ | $(80,80)$ | $(128,128)$ | pixels
Noise variance, $\sigma_{w}^{2}$ | $0.001$ | $0.0025$ | $0.0012$ | —

###  V-D Synthetic Data Results

Figure 4: Depth and reflectivity images of the synthetic toy train image
reconstructed from nine looks by the standard backprojection method (Speckle
Average) and CLAMP with and without aperture modelling over a variety of zero-
padding factors. In each case, the CLAMP reconstructions present less speckle
and improved contrast over the Speckle Average image, and the CLAMP with
aperture model reconstruction greatly improves image sharpness and ability to
resolve finer features, especially in the $1.5\times$ and $2\times$ Nyquist
sampling rate cases (e.g., within the red circle).  Figure 5: Fourier Shell
Correlation plots for synthetic data reconstructions sampled at twice the
Nyquist rate. The FSC curves describe the correlation in frequencies with the
ground truth; larger values indicate higher correlation. The solid lines
represent the mean values over the ten independent data simulations. The
shaded regions how the range of values across all ten simulations. The CLAMP
with aperture model reconstruction generally displays a larger FSC and
intersects the $1/2$-bit threshold curve at a higher frequency.

In Figure 4, we compare the speckle-averaged image and CLAMP reconstructions
with and without an aperture model. The images shown in this figure are the
reconstructions from one of the ten simulated experiments, chosen at random.
For each method, we reconstruct an image at three sampling rates. The sampling
rate is determined in the pre-processing step (Figure 2(c)) by zero-padding by
the factor $q$ in each dimension. When $q=1$, the reconstructed image is
sampled at the Nyquist rate, as labeled. Increasing this factor increases the
sampling rate of the image and allows for more detailed inspection of high
frequency components in the reconstructions. We show reconstructions using
$q=1,1.5,$ and $2$, i.e. sampling at $1\times$, $1.5\times$, and $2\times$ the
Nyquist rate. Note that increasing the sampling rate does not increase the
amount of data, nor does it increase the inherent resolution of the image.
However, at the higher sampling rates, the CLAMP with aperture model
reconstructions display finer, high-resolution features such as those in the
red circle, that are not present in the other reconstructions.

To quantify this improvement in resolution, Figure 5 shows the FSC curve for
the each reconstruction sampled at $2\times$ the Nyquist rate. The solid line
represents the average FSC values across all ten simulated experiments, while
the shaded region illustrates the range of these values across the
experiments. Also shown in the plot is the $1/2$-bit threshold curve,
intersection with which is a common definition of resolution. As one might
expect, the FSC values from the Speckle Average and CLAMP without aperture
model reconstructions closely resemble each other and intersect the 1/2-bit
threshold curve at approximately $0.5$. However, unequivocally, the full CLAMP
reconstruction shows higher FSC values, and intersects the threshold at a
higher frequency. This suggests that including the aperture model in CLAMP
improves the resolution of the reconstruction, which is consistent with the
images shown in Figure 4.

In Table II, we show the point cloud NRMSE and Euclidean distance error for
each reconstruction and sampling rate. The numbers presented in the table
represent the average across the ten simulated experiments. The standard
deviation across the simulated experiments was small at less than $1\%$ of the
average in each case. For each sampling rate, the CLAMP reconstructions
greatly outperform the Speckle Average reconstruction. Without over-sampling,
the metrics for the two CLAMP reconstructions are comparable. However, at
higher sampling rates, the CLAMP with aperture model reconstruction shows the
lowest NRMSE and Euclidean distance error.

TABLE II: Mean NRMSEs and Euclidean distances (meters) between point cloud reconstructions and ground truth point cloud in our synthetic data experiment. |  | Over-sampling
---
Factor
| Speckle
---
Average
| CLAMP
---
(no aperture
model)
| CLAMP
---
(with aperture
model)
Reflectivity NRMSE | None | 0.867 | 0.704 | 0.722
1.5$\times$ | 0.847 | 0.682 | 0.427
2$\times$ | 0.796 | 0.672 | 0.433
Euclidean Distance | None | 0.017 | 0.012 | 0.013
1.5$\times$ | 0.018 | 0.013 | 0.008
2$\times$ | 0.018 | 0.013 | 0.009

###  V-E Experimental Data Results

Figure 6: Left: Depth and reflectivity images of a hexagonal nut Speckle
Average and CLAMP reconstructions from nine looks of data. Before processing,
the data is zero-padded by various amounts. This results in various sampling
rates, labeled on the left. The CLAMP reconstructions with aperture modelling
show the images with the least speckle and sharpest edge boundary. Top right:
The reflectivity profile along the depth dimension on the boundary of the
thread of the nut (at the pixel circled in red). The CLAMP with aperture model
(green) shows distinct, resolvable peaks corresponding to two adjacent
threads, while the other methods do not. Bottom right: The reflectivity values
of the pixels along the red lines in Figure 6. The CLAMP with aperture model
(green) reconstruction shows deep and regular peaks and valleys when compared
to the other reconstructions.

In Figure 6, we show the depth and reflectivity maps of reconstructions of a
hexagonal nut from nine looks at various sampling rates. For each sampling
rate, the CLAMP reconstructions show large improvements in reducing speckle
noise. At higher sampling rates, the full CLAMP method produces a sharper
image. This is qualitatively verified by two 1D reflectivity profiles of the
$2\times$ Nyquist images, shown in Figure 6. One profile (top) shows the
reflectivity along the depth dimension at a pixel (circled in red) that is on
the boundary of a thread. One would expect to see two peaks — one for each
thread. In the methods that do not account for the aperture, these two peaks
are barely resolvable, while in the full CLAMP reconstruction, two completely
resolved peaks are apparent. The second reflectivity profile shows the
reflectivity in the red line across the threads. The full CLAMP reconstruction
shows deep and sharp valleys at a regular interval, as one would expect from a
precisely machined object. In comparison, the other methods show wide, shallow
or irregular peaks and valleys.

Figure 7: Depth and reflectivity images of a toy car reconstructed from nine
looks. The CLAMP reconstruction with aperture modelling has the least amount
of speckle and most discernible features.

Figure 8: Convergence error and $\mu$-residual of CLAMP with aperture modeling
for synthetic and experimental data. CLAMP converged to an error of less than
$10^{-3}$ in all cases.

Figure 7 shows reconstructions from nine looks at a toy car at twice the
Nyquist sampling rate. Relative to the other methods, the full CLAMP
reconstruction reconstructs and image with a low amount of speckle and highly
resolvable features.

In Figure 8, we show the the convergence error of CLAMP, defined in (24), and
a normalized $\mu$-residual averaged across looks, defined in (25), at each
iteration of CLAMP for each dataset with $q=2$. All instances of CLAMP
converge to an error of less than $10^{-3}$. Similarly, the $\mu$-residual
converges, indicating our iterative approach to solving (12) is converging to
the exact solution.

##  VI Conclusion

In this paper, we introduced CLAMP, a majorized PnP algorithm for multi-look
coherent LIDAR image reconstruction. CLAMP is built on the MACE framework and
consists of three main components: (1) an accurate surrogate physics-based
model of coherent LIDAR, (2) a computationally efficient method to invert the
blurring effects of the aperture and (3) a 2.5D neural network denoiser.
Together, these components enable CLAMP to produce high-resolution images with
reduced speckle.

To demonstrate its effectiveness, we applied CLAMP to synthetic and measured
coherent LIDAR data. Our results show that CLAMP produces high-quality 3D
images with reduced speckle and noise compared to speckle-averaging, while
maintaining fidelity to multi-look data. Furthermore, ablating the aperture
model in CLAMP highlights the importance of our method in improving
resolution.

Finally, we formalized the use of surrogate optimization in the MACE framework
and proved convergence to the exact consensus equilibrium solution for a
general class of surrogate functions. Our formalization provides a theoretical
foundation for the use of majorization-minimization in consensus optimization
problems, is general and applicable to many other imaging problems, and
includes many existing image reconstruction algorithms as special cases.

##  Appendix A Quadratic Surrogate

To further reduce the computational cost of the proximal operator in (19), we
follow [47] and construct a quadratic first-order surrogate,
$\widehat{Q}_{\ell}$, of $\hat{f}_{\ell}$, which is then a first-order
surrogate of the original function $f_{\ell}$. For the remainder of this
section, we drop the $\ell$ subscript for simplicity, but note that the
following surrogate is formed for each look $\ell=1,\ldots,L$.

To construct the surrogate, we use the approach in [48, Section IV], which has
two primary components: determining the interval of validity for the
surrogate, which depends on $r^{\prime}$ and $v$, and determining the
surrogate itself. We first obtain a convex surrogate, $Q$, of $\hat{f}$ by a
first-order approximation at $r^{\prime}$ of the logarithm in (18):

$Q\left(r;r^{\prime}\right)=\sum\limits_{j=1}^{n}\frac{r_{j}}{r^{\prime}_{j}}+\frac{1}{r_{j}}\left(c_{j}+\left|\mu_{j}\right|^{2}\right).$
(28)

Treating each $j$ separately, it suffices to form a surrogate the function

$\rho_{j}(r_{j})=\frac{r_{j}}{r_{j}^{\prime}}+\frac{c_{j}+\left|\mu_{j}\right|^{2}}{r_{j}}+\frac{1}{2\sigma^{2}}(r_{j}-v_{j})^{2}.$
(29)

We use (11) and (12) of [48] to fit a quadratic surrogate of $\rho_{j}$ by a
linear interpolation of its derivatives given by

$\widehat{\rho}_{j}(r_{j})=\frac{a_{j}}{2}r_{j}^{2}+b_{j}r_{j}.$ (30)

Summing over $j$ yields the surrogate $\widehat{Q}$ for $Q$.

The surrogate $\widehat{\rho}_{j}$ is constructed by a linear interpolation of
the derivative $\rho^{\prime}_{j}$ at two points: $r_{j}^{\prime}$ and
$\xi_{j}$. The point $\xi_{j}$ is chosen so that it is in the direction of the
minimum of $\rho_{j}$ from $r_{j}^{\prime}$. When
$\rho_{j}^{\prime}(r_{j}^{\prime})>0$, we will need $r_{j}<r_{j}^{\prime}$ to
minimize $\rho$, while if $\rho_{j}^{\prime}(r_{j}^{\prime})<0$ we will need
$r_{j}>r_{j}^{\prime}$. Hence we select a scaling factor $\beta>1$ and choose
$\xi_{j}=r_{j}^{\prime}/\beta$ in the first case and $\xi_{j}=\beta
r_{j}^{\prime}$ in the second. This defines an interval of validity of the
surrogate between $r_{j}^{\prime}$ and $\xi_{j}$, which we denote
$\mathcal{I}_{j}$.

The reflectivity update in (19) is then replaced by the proximal map of
$\widehat{Q}$ at $v$, which is

$r^{\prime}\leftarrow\underset{r_{j}\in\mathcal{I}_{j},\forall
j}{\operatorname{argmin}}\left\\{\widehat{Q}(r;r^{\prime})+\frac{1}{2\sigma^{2}}\
\lVert{r-v}\rVert^{2}\right\\}.$ (31)

Since the objective is separable, this simplifies to

$r^{\prime}_{j}\leftarrow\frac{v_{j}-\sigma^{2}b_{j}}{1+\sigma^{2}a_{j}}$ (32)

but clipped to lie in the interval $\mathcal{I}_{j}$. Finally, we empirically
choose $\beta=1+\frac{2}{\log(k)}$ where $k$ is the CLAMP iteration number,
which ensures that $\beta$ decays to 1 more slowly than the convergence of
$r^{\prime}$.

##  Appendix B 2.5D CNN Architecture

Our 2.5D CNN architecture is similar to DnCNN [49]. With 17 total 2D
convolutional layers, the first layer has 64 filters of size $3\times 3\times
17$, spectral normalization, and ReLU activation. The next 15 layers have 64
filters of size $3\times 3\times 64$, batch normalization, spectral
normalization, and ReLU activation. The final layer has 1 filter of size
$3\times 3\times 64$ and spectral normalization. The network uses a skip
connection between input and output and subtracts estimated noise from the
input image to produce the denoised image. Spectral normalization is used to
satisfy the Lipschitz condition in [50], which aids the convergence of PnP
algorithms.

The denoiser was trained on 800 models from the ShapeNet repository [39]
rendered on a $256\times 256\times 256$ grid. We used the Adam optimizer with
an exponentially decaying learning rate and the MSE loss function. The initial
learning rate was set to $10^{-2}$ and decayed by a factor of 0.025 every 300
training steps. The batch size was set to 24 and the network was trained for
100 epochs on a single NVIDIA A100 GPU. An additional 200 images were used for
validation, and network weights were saved when the validation loss was lowest
during training.

##  Appendix C Proof of Theorem 1

In this section, we prove the convergence of majorized MACE to a fixed point
of the original MACE equation when the surrogate functions meet the conditions
of Definition 1.

We begin by stating a few lemmas before proving the main results. First, a
crucial, perhaps intuitive, fact is that the objective function and its
surrogate have the same subgradients at the point of approximation. This is
stated in Lemma 1 and is a result of the smoothness assumption in Definition
1.

######  Lemma 1.

Let $\hat{f}\in\mathcal{S}_{L,p}\left(f,\xi\right)$ (see Definition 1). Then
$\partial f(\xi)=\partial\hat{f}(\xi)$.

###### Proof.

Proof is given in the supplemental material. We note that this result holds
even if $\hat{f}$ is not strongly convex. ∎

Another critical fact is given in Lemma 2, which is adapted from [31, Lemma
2].

######  Lemma 2.

Let $\hat{f}\in\mathcal{S}_{L,p}\left(f,\xi\right)$ (see Definition 1), and
let $u\in\partial\hat{f}(x)$ be any subgradient of $\hat{f}$ at $x$. Then for
any $y$,

$\displaystyle f(x)-f(y)-$ $\displaystyle\left\langle u,x-y\right\rangle$
$\displaystyle\leq\frac{1}{2}\left(L\left\lVert{y-\xi}\right\rVert^{2}-p\left\lVert{x-y}\right\rVert^{2}\right).$

###### Proof.

Proof is given in the supplemental material. ∎

We are now ready to prove part (i) of Theorem 1.

###### Proof of Theorem 1 part (i).

Let $(\mathbf{w}^{*},\mathbf{r}^{*})$ be a fixed point of (23). Since
$\mathbf{r}^{*}=\hat{\mathbf{F}}^{*}\mathbf{w}^{*}$, where
$\hat{\mathbf{F}}^{*}$ is the stacked proximal operator of surrogates
$\hat{f}_{i}(x;r_{i}^{*})$, we have

$-\frac{1}{\sigma^{2}}\left(r_{i}^{*}-w_{i}^{*}\right)\in\partial\hat{f}_{i}(r^{*};r^{*}),$
(33)

and by Lemma 1,

$-\frac{1}{\sigma^{2}}\left(r_{i}^{*}-w_{i}^{*}\right)\in\partial
f_{i}(r^{*}),$ (34)

for $i=1,\dots,N$.

Thus, $\mathbf{r}^{*}=\mathbf{F}\mathbf{w}^{*}$. Since
$\mathbf{w}^{*}=\left(2\mathbf{G}-I\right)\left(2\hat{\mathbf{F}}^{*}-I\right)\mathbf{w}^{*}$,
we get
$\mathbf{w}^{*}=\left(2\mathbf{G}-I\right)\left(2\mathbf{F}-I\right)\mathbf{w}^{*}$,
proving $\mathbf{w}^{*}$ is a MACE solution. ∎

In order to prove Theorem 1 part (ii), we prove convergence of an algorithm
that is equivalent to Algorithm 3. First, it is noted in [28] and elsewhere
that using Mann iterations to solve the MACE equations, as written in (9), is
equivalent up to a change of variables to a form of ADMM [30].

We present the equivalency and prove convergence for $\rho=1/2$, though with
extra bookkeeping, and using the appropriate variant of ADMM, one can
generalize this result for $\rho\in(0,1)$.

To begin, consider one iteration of Algorithm 3, written as in (23), using
$\rho=1/2$,

$\displaystyle\mathbf{w}^{(k+1)}$
$\displaystyle=\frac{1}{2}(2\mathbf{G}-\mathbf{I})(2\hat{\mathbf{F}}^{(k)}-\mathbf{I})\mathbf{w}^{(k)}+\frac{1}{2}\mathbf{w}^{(k)},$
$\displaystyle\mathbf{r}^{(k+1)}$
$\displaystyle=\hat{\mathbf{F}}^{(k)}\mathbf{w}^{(k)},$

where we recall that
$\hat{\mathbf{F}}^{(k)}\mathbf{w}^{(k)}=\hat{\mathbf{F}}^{(k)}(\mathbf{w}^{(k)};\mathbf{r}^{(k)})$.
Since $\mathbf{G}$ is linear and $\mathbf{G}\mathbf{G}=\mathbf{G}$, some
algebra shows that applying $\mathbf{G}$ to both sides of the first equation
yields $\mathbf{G}\mathbf{w}^{(k+1)}=\mathbf{G}\mathbf{r}^{(k+1)}$, and hence
$\mathbf{G}\mathbf{w}^{(k+1)}=\mathbf{G}\hat{\mathbf{F}}^{(k)}\mathbf{w}^{(k)}$
from the second equation.

Using this to re-express the first equation, recalling that
$\overline{\mathbf{r}}$ is one averaged image in $\mathbf{G}\mathbf{r}$, and
noting that $\mathbf{w}$ and $\mathbf{r}$ may be updated in either order, the
$i$-th component $(r_{i},w_{i})$ satisfies

$\displaystyle r_{i}^{(k+1)}$ $\displaystyle=\hat{F}_{i}^{(k)}w_{i}^{(k)},$
(35) $\displaystyle w_{i}^{(k+1)}$
$\displaystyle=w_{i}^{(k)}+(\overline{\mathbf{r}}^{(k+1)}-r_{i}^{(k+1)})$ (36)
$\displaystyle\quad+(\overline{\mathbf{r}}^{(k+1)}-\overline{\mathbf{r}}^{(k)}).$
(37)

Now, we introduce the variables
$u_{i}^{(k)}=\overline{\mathbf{r}}^{(k)}-w_{i}^{(k)}$ and
$z^{(k+1)}=\overline{\mathbf{r}}^{(k+1)}+\overline{\mathbf{u}}^{(k)}$. Since
$\overline{\mathbf{w}}^{(k+1)}=\overline{\mathbf{r}}^{(k+1)}$, we have
$\overline{\mathbf{u}}^{(k+1)}=0$ for $k\geq 1$, and so for $k\geq 3$, the
updates (35)–(37) can be written equivalently as

$\displaystyle r_{i}^{(k+1)}$
$\displaystyle=\hat{F}_{i}^{(k)}\left(z^{(k)}-u_{i}^{(k)}\right),$ (38)
$\displaystyle z^{(k+1)}$
$\displaystyle=\frac{1}{N}\sum\limits_{i=1}^{N}\left(r_{i}^{(k+1)}+u_{i}^{(k)}\right),$
(39) $\displaystyle u_{i}^{(k+1)}$
$\displaystyle=u_{i}^{(k)}+(r_{i}^{(k+1)}-z^{(k+1)}),$ (40)

where (40) follows from (37) using
$w_{i}^{(k)}=\overline{\mathbf{r}}^{(k)}-u_{i}^{(k)}$ and simplifying. The
updates (38)–(40) are a form of consensus ADMM [30, Chapter 7] to solve the
problem,

minimize $\displaystyle f(\mathbf{r})=\sum_{i=1}^{N}f_{i}(r_{i})$ (41) subject
to $\displaystyle r_{i}=z,$

except the updates of $r_{i}$ in (38) are given by proximal maps of surrogate
functions instead of proximal maps of $f_{i}$.

In Theorem 2, we prove that despite this difference, the iterations (38)-(40),
converge to a Karush-Kuhn-Tucker (KKT) point of (41), which is sufficient for
optimality. A KKT point can be found from the Lagrangian of (41),

$\mathcal{L}(\mathbf{r},z,\mathbf{u})=f(\mathbf{r})+\frac{1}{\sigma^{2}}\sum_{i=1}^{N}\langle
u_{i},r_{i}-z\rangle,$ (42)

where $\mathbf{u}$ is the dual variable, and $\langle\cdot,\cdot\rangle$ is
the inner-product. Specifically, the point
$(\mathbf{r}^{*},z^{*},\mathbf{u}^{*})$ is a KKT point if and only if

$\displaystyle-\frac{1}{\sigma^{2}}u_{i}^{*}\in$ $\displaystyle\partial
f_{i}(r_{i}^{*}),$ (43) $\displaystyle r_{i}^{*}$ $\displaystyle=z^{*},$ (44)

for $i=1,\dots,N$.

######  Theorem 2.

Assume the hypotheses of Theorem 1, and that (41) has a KKT point. Then the
sequence
$\left\\{\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)\right\\}$
generated by the iterations (38)-(40) converges to a KKT point of (41).

To prove Theorem 2, we first use Proposition 1 to characterize the
subgradients of the surrogate functions. Then in Lemma 3, we state the vital
property that the distance between the iterates and a KKT point decreases at
each iteration.

######  Proposition 1.

Assume the hypotheses of Theorem 1 and let
$\left\\{\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)\right\\}$ be
generated by (38)-(40). Then

$\displaystyle-\frac{1}{\sigma^{2}}\left(r_{i}^{(k+1)}-z^{(k)}+u_{i}^{(k)}\right)$
$\displaystyle\in\partial\hat{f}_{i}^{(k)}(r_{i}^{(k+1)}),$ (45)

for $i=1,\ldots,N$.

###### Proof.

This follows from the definition of the proximal map in (38). ∎

To simplify notation, we define

$u_{i}^{\left(k+\frac{1}{2}\right)}\triangleq
u_{i}^{(k)}+\left(r_{i}^{(k+1)}-z^{(k)}\right)$ (46)

which can be thought of as an intermediate update of the dual variable after
the update of $\mathbf{r}$ in (38) but before the update of $z$ in (39). We
can then rewrite the subgradients in (45) as

$\displaystyle-\frac{1}{\sigma^{2}}u_{i}^{\left(k+\frac{1}{2}\right)}$
$\displaystyle\in\partial\hat{f}^{(k)}_{i}(r_{i}^{(k+1)}).$ (47)

######  Lemma 3.

Assume the hypotheses of Theorem 1, and that (41) has a KKT point,
$\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)$. Then the sequence
$\left\\{\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)\right\\}$
generated by the iterations (38)-(40) satisfies

$E^{(k+1)}\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)\leq
E^{(k)}\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right),$ (48)

where

$\displaystyle E^{(k)}\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)$
$\displaystyle=\sum_{i=1}^{N}\Big{(}L_{i}\left\lVert{r_{i}^{(k)}-r_{i}^{*}}\right\rVert^{2}$
(49)
$\displaystyle\quad+\frac{1}{\sigma^{2}}\left\lVert{z^{(k)}-z^{*}}\right\rVert^{2}+\frac{1}{\sigma^{2}}\left\lVert{u_{i}^{(k)}-u_{i}^{*}}\right\rVert^{2}\Big{)}.$
(50)

###### Proof.

Proof is given in the supplemental material. ∎

We are now ready to state and prove convergence of the iterations (38)-(40).

###### Proof of Theorem 2.

Let $(\mathbf{r}^{*},z^{*},\mathbf{u}^{*})$ be a KKT point of (41). By Lemma
3, we have that
$\left\\{\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)\right\\}$ is
bounded. Let
$\left\\{\left(\mathbf{r}^{(k_{j})},z^{(k_{j})},\mathbf{u}^{(k_{j})}\right)\right\\}$
be a convergent subsequence with limit point
$\left(\mathbf{r}^{\dagger},z^{\dagger},\mathbf{u}^{\dagger}\right)$. We will
show that this limit point is also a KKT point of (41).

From (98)-(103) in the proof of Lemma 3, the iterates satisfy

$\displaystyle\frac{1}{2\sigma^{2}}$
$\displaystyle\sum_{i=1}^{N}\left\lVert{z^{(k+1)}-z^{(k)}}\right\rVert^{2}+\left\lVert{u_{i}^{(k+1)}-u_{i}^{(k)}}\right\rVert^{2}$
(51)
$\displaystyle\leq\frac{1}{2}\left(E^{(k)}\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)-E^{(k+1)}\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)\right).$
(52)

Using Lemma 3, summing both sides of (51)-(52) over all $k$, and discarding
the telescoping terms, we get $\left\lVert z^{(k+1)}-z^{(k)}\right\rVert\to 0$
and $\left\lVert u_{i}^{(k+1)}-u_{i}^{(k)}\right\rVert\to 0$ as $k\to\infty$.
Since $z^{(k_{j})}\rightarrow z^{\dagger}$ as $j\rightarrow\infty$, we can set
$k=k_{j}$ and take $j\rightarrow\infty$ to get $z^{(k_{j}+1)}\rightarrow
z^{\dagger}$. Likewise, $u_{i}^{(k_{j})}\rightarrow u_{i}^{\dagger}$.

Rearranging (40), we have

$r_{i}^{(k_{j}+1)}=u_{i}^{(k_{j}+1)}-u_{i}^{(k_{j})}+z^{(k_{j}+1)},$

so $r_{i}^{(k_{j}+1)}$ converges to $z^{\dagger}$, which we denote also as
$r_{i}^{\dagger}$. Lastly, $u_{i}^{(k_{j}+\frac{1}{2})}\to u_{i}^{\dagger}$
can be seen by taking the limit of (46) along the subsequence $k=k_{j}$.

Since
$-\frac{1}{\sigma^{2}}u_{i}^{(k_{j}+\frac{1}{2})}\in\partial\hat{f}_{i}^{(k_{j})}(r_{i}^{(k_{j}+1)})$,
we have by Lemma 2 that for all $\mathbf{r}$,

$f(\mathbf{r}^{(k_{j}+1)})-f(\mathbf{r})+\frac{1}{\sigma^{2}}\sum_{i=1}^{N}\left\langle
u_{i}^{(k_{j}+\frac{1}{2})},r_{i}^{(k_{j}+1)}-r_{i}\right\rangle\leq 0.$ (53)

Since $f$ is convex and finite on ${\mathbb{R}}^{n}$ and $r_{i}^{(k_{j}+1)}\to
r_{i}^{\dagger}$, we have that $f(\mathbf{r}^{(k_{j}+1)})\to
f(\mathbf{r}^{\dagger})$, and so taking the limit of the above inequality
along the subsequence gives

$f(\mathbf{r}^{\dagger})-f(\mathbf{r})+\frac{1}{\sigma^{2}}\sum_{i=1}^{N}\left\langle
u_{i}^{\dagger},r_{i}^{\dagger}-r_{i}\right\rangle\leq 0,$ (54)

which implies that $-\frac{1}{\sigma^{2}}u_{i}^{\dagger}\in\partial
f_{i}(r_{i}^{\dagger})$ for $i=1,\dots,N$. This completes the proof that
$\left(\mathbf{r}^{\dagger},z^{\dagger},\mathbf{u}^{\dagger}\right)$ is a KKT
point.

Finally, we show that
$\left(\mathbf{r}^{\dagger},z^{\dagger},\mathbf{u}^{\dagger}\right)$ is the
only limit point of the sequence
$\left\\{\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)\right\\}$.
First, we replace $\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)$ with
$\left(\mathbf{r}^{\dagger},z^{\dagger},\mathbf{u}^{\dagger}\right)$ in Lemma
3. Then for any $k>k_{j}$, we have

$E^{(k)}\left(\mathbf{r}^{\dagger},z^{\dagger},\mathbf{u}^{\dagger}\right)\leq
E^{(k_{j})}\left(\mathbf{r}^{\dagger},z^{\dagger},\mathbf{u}^{\dagger}\right).$
(55)

Since the right hand side converges to 0 as $j\to\infty$, the definition of
$E^{(k)}$ implies that
$\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)\to\left(\mathbf{r}^{\dagger},z^{\dagger},\mathbf{u}^{\dagger}\right)$
as $k\to\infty$. ∎

To end, we prove Theorem 1 part (ii).

###### Proof of Theorem 1 part (ii).

As noted after equation (37), by introducing the variables
$u_{i}^{(k)}=\overline{\mathbf{r}}^{(k)}-w_{i}^{(k)}$, and
$z^{(k+1)}=\overline{\mathbf{r}}^{(k+1)}+\overline{\mathbf{u}}^{(k)}$, the
iterates of (23) are equivalent to the iterates (38)-(40). By Theorem 2,
$\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)$ converges to a KKT
point, $\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)$, of (41).
Resubsituting the change of variables, we get
$(\mathbf{w}^{*},\mathbf{r}^{*})$ is a fixed point of the system (23), which
by Theorem 1 part (i) is a solution to the exact MACE formulation. ∎

## References

  * [1] M. F. Spencer, “Spatial Heterodyne,” in Encyclopedia of Modern Optics, vol. 4, Elsevier, 2018.
  * [2] N. H. Farhat, “Holography, wave‐length diversity and inverse scattering,” AIP Conference Proceedings, vol. 65, no. 1, p. 627–642, 1980.
  * [3] J. C. Marron and T. J. Schulz, “Three-dimensional, fine-resolution imaging using laser frequency diversity,” Optics Letters, vol. 17, no. 4, p. 285, 1992.
  * [4] J. C. Marron and K. S. Schroeder, “Three-dimensional lensless imaging using laser frequency diversity,” Applied Optics, vol. 31, no. 2, p. 255, 1992\.
  * [5] J. W. Stafford, B. D. Duncan, and D. J. Rabb, “Phase gradient algorithm method for three-dimensional holographic ladar imaging,” Applied Optics, vol. 55, no. 17, p. 4611, 2016.
  * [6] M. T. Banet and J. R. Fienup, “Effects of speckle decorrelation on motion-compensated, multi-wavelength, 3D digital holography,” in Unconventional Imaging and Adaptive Optics 2022 (J. J. Dolne and M. F. Spencer, eds.), vol. 12239, p. 122390C, International Society for Optics and Photonics, SPIE, 2022.
  * [7] W. Farriss, Iterative Phase Estimation Algorithms in Interferometric Systems.  PhD thesis, University of Rochester, 2021.
  * [8] W. E. Farriss, J. R. Fienup, J. W. Stafford, and N. J. Miller, “Sharpness-based correction methods in holographic aperture ladar (HAL),” in Unconventional and Indirect Imaging, Image Reconstruction, and Wavefront Sensing 2018, vol. 10772, pp. 178–185, SPIE, 2018.
  * [9] M. T. Banet and J. R. Fienup, “Image sharpening on 3D intensity data in deep turbulence with scintillated illumination,” in Unconventional Imaging and Adaptive Optics 2021, vol. 11836, pp. 116–124, SPIE, 2021\.
  * [10] M. T. Banet, J. R. Fienup, J. D. Schmidt, and M. F. Spencer, “3D multi-plane sharpness metric maximization with variable corrective phase screens,” Applied Optics, vol. 60, no. 25, p. G243, 2021.
  * [11] C. J. Pellizzari, M. F. Spencer, and C. A. Bouman, “Phase-error estimation and image reconstruction from digital-holography data using a Bayesian framework,” Journal of the Optical Society of America A, vol. 34, no. 9, p. 1659, 2017.
  * [12] C. J. Pellizzari, R. Trahan, H. Zhou, S. Williams, S. E. Williams, B. Nemati, M. Shao, and C. A. Bouman, “Synthetic Aperature LADAR: A Model-Based Approach,” IEEE Transactions on Computational Imaging, vol. 3, no. 4, pp. 901–916, 2017.
  * [13] C. J. Pellizzari, M. F. Spencer, and C. A. Bouman, “Coherent Plug-and-Play: Digital Holographic Imaging Through Atmospheric Turbulence Using Model-Based Iterative Reconstruction and Convolutional Neural Networks,” IEEE Transactions on Computational Imaging, vol. 6, pp. 1607–1621, 2020.
  * [14] C. J. Pellizzari, M. F. Spencer, and C. A. Bouman, “Imaging through distributed-volume aberrations using single-shot digital holography,” Journal of the Optical Society of America A, vol. 36, no. 2, p. A20, 2019.
  * [15] C. J. Pellizzari, R. Trahan, H. Zhou, S. Williams, S. E. Williams, B. Nemati, M. Shao, and C. A. Bouman, “Optically coherent image formation and denoising using a plug and play inversion framework,” Applied Optics, vol. 56, no. 16, p. 4735, 2017.
  * [16] L. E. Baum and T. Petrie, “Statistical Inference for Probabilistic Functions of Finite State Markov Chains,” The Annals of Mathematical Statistics, vol. 37, no. 6, pp. 1554 – 1563, 1966.
  * [17] L. E. Baum, T. Petrie, G. Soules, and N. Weiss, “A Maximization Technique Occurring in the Statistical Analysis of Probabilistic Functions of Markov Chains,” The Annals of Mathematical Statistics, vol. 41, no. 3, pp. 698 – 704, 1970.
  * [18] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum Likelihood from Incomplete Data Via the _EM_ Algorithm,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 39, no. 1, pp. 1–22, 1977.
  * [19] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, “Plug-and-Play priors for model based reconstruction,” in 2013 IEEE Global Conference on Signal and Information Processing, pp. 945–948, IEEE, 2013.
  * [20] S. Sreehari, S. V. Venkatakrishnan, B. Wohlberg, L. F. Drummy, J. P. Simmons, and C. A. Bouman, “Plug-and-Play Priors for Bright Field Electron Tomography and Sparse Interpolation,” IEEE Transactions on Computational Imaging, pp. 1–1, 2016.
  * [21] C. A. Bouman, Foundations of Computational Imaging: A Model-Based Approach.  Other Titles in Applied Mathematics, Society for Industrial and Applied Mathematics, 2022.
  * [22] Y. Altmann, X. Ren, A. McCarthy, G. S. Buller, and S. McLaughlin, “Robust bayesian target detection algorithm for depth imaging from sparse single-photon data,” IEEE Transactions on Computational Imaging, vol. 2, no. 4, pp. 456–467, 2016.
  * [23] J. Rapp and V. K. Goyal, “A few photons among many: Unmixing signal and noise for photon-efficient active imaging,” IEEE Transactions on Computational Imaging, vol. 3, no. 3, pp. 445–459, 2017.
  * [24] J. Tachella, Y. Altmann, X. Ren, A. McCarthy, G. S. Buller, S. McLaughlin, and J.-Y. Tourneret, “Bayesian 3d reconstruction of complex scenes from single-photon lidar data,” SIAM Journal on Imaging Sciences, vol. 12, no. 1, pp. 521–550, 2019.
  * [25] J. Tachella, Y. Altmann, N. Mellado, A. McCarthy, R. Tobin, G. S. Buller, J.-Y. Tourneret, and S. McLaughlin, “Real-time 3d reconstruction from single-photon lidar data using plug-and-play point cloud denoisers,” Nature Communications, vol. 10, p. 4984, Nov. 2019.
  * [26] A. Ziabari, D. H. Ye, S. Srivastava, K. D. Sauer, J.-B. Thibault, and C. A. Bouman, “2.5d deep learning for ct image reconstruction using a multi-gpu implementation,” in 2018 52nd Asilomar Conference on Signals, Systems, and Computers, pp. 2044–2049, 2018.
  * [27] S. Majee, T. Balke, C. A. J. Kemp, G. T. Buzzard, and C. A. Bouman, “4d x-ray ct reconstruction using multi-slice fusion,” in 2019 IEEE International Conference on Computational Photography (ICCP), pp. 1–8, 2019.
  * [28] G. T. Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman, “Plug-and-Play Unplugged: Optimization-Free Reconstruction Using Consensus Equilibrium,” SIAM Journal on Imaging Sciences, vol. 11, no. 3, pp. 2001–2020, 2018.
  * [29] D. Tucker and L. C. Potter, “Speckle Suppression in Multi-Channel Coherent Imaging: A Tractable Bayesian Approach,” IEEE Transactions on Computational Imaging, vol. 6, pp. 1429–1439, 2020.
  * [30] S. Boyd, “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers,” Foundations and Trends® in Machine Learning, vol. 3, no. 1, pp. 1–122, 2010.
  * [31] C. Lu, J. Feng, S. Yan, and Z. Lin, “A unified alternating direction method of multipliers by majorization minimization,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 40, no. 3, pp. 527–541, 2018.
  * [32] J. Mairal, “Incremental majorization-minimization optimization with application to large-scale machine learning,” SIAM Journal on Optimization, vol. 25, no. 2, pp. 829–855, 2015.
  * [33] J. Mairal, “Optimization with first-order surrogate functions,” in Proceedings of the 30th International Conference on Machine Learning (S. Dasgupta and D. McAllester, eds.), vol. 28 of Proceedings of Machine Learning Research, (Atlanta, Georgia, USA), pp. 783–791, PMLR, 17–19 Jun 2013\.
  * [34] D. R. H. Kenneth Lange and I. Yang, “Optimization transfer using surrogate objective functions,” Journal of Computational and Graphical Statistics, vol. 9, no. 1, pp. 1–20, 2000.
  * [35] T. Allen, D. Rabb, G. T. Buzzard, and C. A. Bouman, “Multi-Agent Consensus Equilibrium for Range Compressed Holographic Surface Reconstruction,” in 21st Coherent Laser Radar Conference, 2022.
  * [36] J. Goodman, Introduction to Fourier Optics.  MaGraw-Hill, 2nd ed., 1996.
  * [37] J. W. Goodman, “Some fundamental properties of speckle*,” Journal of the Optical Society of America, vol. 66, no. 11, p. 1145, 1976.
  * [38] J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications, Second Edition.  SPIE, 2020.
  * [39] A. X. Chang, T. Funkhouser, L. Guibas, P. Hanrahan, Q. Huang, Z. Li, S. Savarese, M. Savva, S. Song, H. Su, J. Xiao, L. Yi, and F. Yu, “ShapeNet: An Information-Rich 3D Model Repository,” Tech. Rep. arXiv:1512.03012 [cs.GR], Stanford University — Princeton University — Toyota Technological Institute at Chicago, 2015.
  * [40] R. Gribonval, “Should Penalized Least Squares Regression be Interpreted as Maximum A Posteriori Estimation?,” vol. 59, no. 5, pp. 2405–2410, 2011.
  * [41] X. Xu, Y. Sun, J. Liu, B. Wohlberg, and U. S. Kamilov, “Provable Convergence of Plug-and-Play Priors With MMSE Denoisers,” vol. 27, pp. 1280–1284, 2020.
  * [42] M. W. Jacobson and J. A. Fessler, “An expanded theoretical treatment of iteration-dependent majorize-minimize algorithms,” IEEE Transactions on Image Processing, vol. 16, no. 10, pp. 2411–2422, 2007.
  * [43] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” IEEE Transactions on Signal Processing, vol. 65, no. 3, pp. 794–816, 2017.
  * [44] J. H. Friedman, J. L. Bentley, and R. A. Finkel, “An algorithm for finding best matches in logarithmic expected time,” ACM Transactions on Mathematical Software, vol. 3, p. 209–226, Sept. 1977.
  * [45] G. Harauz and M. van Heel, “Exact filters for general geometry three dimensional reconstruction,” Optik 73, pp. 146–156, 1986.
  * [46] M. v. Heel and M. Schatz, “Fourier shell correlation threshold criteria,” Journal of Structural Biology, vol. 151, no. 3, p. 250–262, 2005.
  * [47] V. Sridhar, S. J. Kisner, S. P. Midkiff, and C. A. Bouman, “Fast algorithms for model-based imaging through turbulence,” in Artificial Intelligence and Machine Learning in Defense Applications II (J. Dijk, ed.), vol. 11543, p. 1154304, International Society for Optics and Photonics, SPIE, 2020\.
  * [48] Z. Yu, J.-B. Thibault, K. Sauer, C. Bouman, and J. Hsieh, “Accelerated Line Search for Coordinate Descent Optimization,” in 2006 IEEE Nuclear Science Symposium Conference Record, vol. 5, pp. 2841–2844, 2006.
  * [49] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang, “Beyond a Gaussian Denoiser: Residual Learning of Deep CNN for Image Denoising,” IEEE Transactions on Image Processing, vol. 26, no. 7, pp. 3142–3155, 2017\.
  * [50] E. Ryu, J. Liu, S. Wang, X. Chen, Z. Wang, and W. Yin, “Plug-and-play methods provably converge with properly trained denoisers,” in Proceedings of the 36th International Conference on Machine Learning (K. Chaudhuri and R. Salakhutdinov, eds.), vol. 97 of Proceedings of Machine Learning Research, pp. 5546–5557, PMLR, 09–15 Jun 2019.

Supplemental Material for “CLAMP: Majorized Plug-and-Play for Coherent 3D
LIDAR Imaging”

## Forward Model Derivation

In this section, we derive the forward model for coherent LIDAR imaging as
presented in the main text. We begin by presenting the continuous model, and
then discretize it to obtain the discrete model (1).

We assume an off-axis digital holographic imaging system, as depicted in
Figure 2, where a linear frequency modulated (LFM) chirped plane wave
illuminates an object with a complex-valued reflectance
$\tilde{g}:\mathbb{R}^{3}\to\mathbb{C}$. We recover the wave in the pupil
plane, $\tilde{y}:\mathbb{R}^{3}\to\mathbb{C}$, with a digital holographic
imaging system. For a full account of this process, we refer the reader to
spencerSpatialHeterodyne.

By standard Fresnel diffraction theory goodmanIntroductionFourierOptics, each
transverse plane in the object’s coordinate system can be propagated to the
entrance pupil by the Fresnel diffraction integral. For the slice of the
object given by $\tilde{g}(\cdot,\cdot,z)$, the field propagated to the
entrance pupil at time $t$, $U_{z}(\xi,\nu;t)$, is given by

$\displaystyle U_{z}(\xi,\nu;z,t)$ $\displaystyle=\frac{e^{2\pi
id_{z}/\lambda_{t}}}{i\lambda_{t}z}\tilde{a}(\xi,\nu,t)\exp\left\\{\frac{2\pi
i}{2\lambda_{t}d_{z}}(\xi^{2}+\nu^{2})\right\\}$
$\displaystyle\iint_{-\infty}^{\infty}\left\\{\tilde{g}(x,y,z)\exp\left\\{\frac{2\pi
i}{2\lambda_{t}d_{z}}(x^{2}+y^{2})\right\\}\right\\}$
$\displaystyle\exp\left\\{-\frac{2\pi i}{\lambda_{t}d_{z}}(\xi x+\nu
y)\right\\}\,dx\,dy,$ (56)

where $(\xi,\nu)$ are the coordinates in the pupil plane, $(x,y,z)$ are the
coordinates of the target, $d_{z}$ is the distance from the object plane to
the pupil plane, $\tilde{a}(\xi,\nu,t)$ is the aperture, and $\lambda_{t}$ is
the instantaneous wavelength of the LFM chirp at time $t$.

To proceed, we assume a small fractional bandwidth and small extent of the
object compared to the range, which allows us to replace $\lambda_{t}d_{z}$
with $\lambda d$. These assumptions are reasonable and are satisfied in many
practical scenarios farrissIterativePhaseEstimation2021, including the one
presented in this paper. With these assumptions, we have

$\displaystyle U_{z}(\xi,\nu;t)$ $\displaystyle=\frac{e^{2\pi
id_{z}/\lambda_{t}}}{i\lambda d}\tilde{a}(\xi,\nu,t)\exp\left\\{\frac{2\pi
i}{2\lambda d}(\xi^{2}+\nu^{2})\right\\}$
$\displaystyle\iint_{-\infty}^{\infty}\left\\{\tilde{g}(x,y,z)\exp\left\\{\frac{2\pi
i}{2\lambda d}(x^{2}+y^{2})\right\\}\right\\}$
$\displaystyle\exp\left\\{-\frac{2\pi i}{\lambda d}(\xi x+\nu
y)\right\\}\,dx\,dy.$ (57)

Next, by a small-angle approximation that $(\xi^{2}+\nu^{2})/\lambda d\approx
0$, we can drop the quadratic phase term outside the integral. Moreover, as
the variable of interest is the magnitude of $\tilde{g}$, we can absorb the
quadratic phase term inside the integral into $\tilde{g}$, and write

$\displaystyle U_{z}(\xi,\nu;t)$ $\displaystyle=\frac{e^{2\pi
id_{z}/\lambda_{t}}}{i\lambda d}\tilde{a}(\xi,\nu,t)$
$\displaystyle\iint_{-\infty}^{\infty}\left\\{\tilde{g}(x,y,z)\exp\left\\{-\frac{2\pi
i}{\lambda d}(\xi x+\nu y)\right\\}\right\\}\,dx\,dy.$ (58)

Finally, the total pupil field $\tilde{y}$ at time $t$ is the combination of
the fields $U_{z}$ from all $z$. Hence, by integrating over the extent of the
object we have

$\displaystyle\tilde{y}(\xi,\nu,t)$ $\displaystyle=\frac{e^{2\pi
id/\lambda_{t}}}{i\lambda
d}\tilde{a}(\xi,\nu,t)\iiint\limits_{\mathbb{R}^{3}}\Bigg{[}\tilde{g}(x,y,z)$
(59) $\displaystyle\exp\left\\{-2\pi i\left(\frac{\xi}{\lambda
d}x+\frac{\nu}{\lambda
d}y+\frac{1}{\lambda_{t}}z\right)\right\\}\Bigg{]}\,dx\,dy\,dz.$ (60)

In other words, $\tilde{y}$ is the 3D Fourier transform of the object
$\tilde{g}$. By normalizing and discretizing the continuous model, we obtain
the discrete model (1).

## Proof of Lemma 1

Lemma 1.  Let $\hat{f}\in\mathcal{S}_{L,p}\left(f,\xi\right)$. Then $\partial
f(\xi)=\partial\hat{f}(\xi)$.

###### Proof of Lemma 1.

The smoothness condition in Definition 1 states $\nabla e(\xi)=0$, where
$e\triangleq f-\hat{f}$. Hence, the directional derivatives of $f$ and
$\hat{f}$ at $\xi$ must be equal in all directions. That is, for all
directions $h$,

$f^{\prime}(\xi;h)=\hat{f}^{\prime}(\xi;h)$ (61)

where

$f^{\prime}(\xi;h)=\lim_{t\to 0^{+}}\frac{f(\xi+th)-f(\xi)}{t},$ (62)

and

$\hat{f}^{\prime}(\xi;h)=\lim_{t\to
0^{+}}\frac{\hat{f}(\xi+th)-\hat{f}(\xi)}{t}.$ (63)

Since subgradients can be specified by the set of all directional derivatives
[Theorem 8.8 and 8.10]recht-wright, we have

$\displaystyle\partial f(\xi)$ $\displaystyle=\left\\{u\mid\langle
u,h\rangle\leq f^{\prime}(\xi;h),\forall h\right\\}$ (64)
$\displaystyle=\left\\{u\mid\langle
u,h\rangle\leq\hat{f}^{\prime}(\xi;h),\forall h\right\\}$ (65)
$\displaystyle=\partial\hat{f}(\xi).$ (66)

∎

## Proof of Lemma 2

Before proving Lemma 2, we state in Lemma 4 a useful property of the error
function, $e\triangleq f-\hat{f}$.

######  Lemma 4.

If $e$ is differentiable with an $L$-Lipschitz continuous gradient, then for
any $\xi,y$,

$\lvert e(y)-e(\xi)-\left\langle\nabla
e(\xi),y-\xi\right\rangle\rvert\leq\frac{L}{2}\left\lVert
y-\xi\right\rVert^{2}.$ (67)

###### Proof.

This is Lemma 1.2.3 in nesterovIntroductoryLecturesConvex2004. ∎

Lemma 2.  Let $\hat{f}\in\mathcal{S}_{L,p}\left(f,\xi\right)$ (see Definition
1), and let $u\in\partial\hat{f}(x)$ be any subgradient of $\hat{f}$ at $x$.
Then for any $y$,

$\displaystyle f(x)-f(y)-$ $\displaystyle\left\langle u,x-y\right\rangle$
$\displaystyle\leq\frac{1}{2}\left(L\left\lVert{y-\xi}\right\rVert^{2}-p\left\lVert{x-y}\right\rVert^{2}\right).$

###### Proof of Lemma 2.

We begin with the fact that $\hat{f}$ is $p$-strongly convex. There are
several equivalent definitions of strong convexity, see for example [Section
9.1.2]boydConvexOptimization2004. One such definition is that for all $x,y$
and subgradients $u\in\partial\hat{f}(x)$, we have

$\hat{f}(y)\geq\hat{f}(x)-\left\langle
u,x-y\right\rangle+\frac{p}{2}\left\lVert{y-x}\right\rVert^{2}.$ (68)

Using this, we have

$\displaystyle f(x)$ $\displaystyle\leq\hat{f}(x)$ (69)
$\displaystyle\leq\hat{f}(y)+\left\langle
u,x-y\right\rangle-\frac{p}{2}\left\lVert{y-x}\right\rVert^{2},$ (70)

where (69) is the majorization property of $\hat{f}$, and (70) follows from
Equation (68). Now, defining $e(y)\triangleq\hat{f}(y)-f(y)$, the fact that
$\hat{f}\in\mathcal{S}_{L,p}\left(f,\xi\right)$ implies that $e(\xi)=0$ and
$\nabla e(\xi)=0$. Hence Lemma 4 implies that
$e(y)\leq\frac{L}{2}\left\lVert{y-\xi}\right\rVert^{2}$. Therefore, continuing
the inequalities in (69)-(70), we have

$\displaystyle f(x)$ $\displaystyle\leq
f(y)+\frac{L}{2}\left\lVert{y-\xi}\right\rVert^{2}+\left\langle
u,x-y\right\rangle-\frac{p}{2}\left\lVert{y-x}\right\rVert^{2},$

completing the proof. ∎

## Proof of Lemma 3

Before we prove Lemma 3, we prove two additional useful lemmas.

######  Lemma 5.

Assume the hypotheses of Proposition 1. Then for any $(\mathbf{r},z)$, we have

$\displaystyle f(\mathbf{r}^{(k+1)})-f(\mathbf{r})$
$\displaystyle\quad+\frac{1}{\sigma^{2}}\sum\limits_{i=1}^{N}\left\langle
u_{i}^{\left(k+\frac{1}{2}\right)},\left(r_{i}^{(k+1)}-r_{i}\right)-\left(z^{(k+1)}-z\right)\right\rangle$
$\displaystyle\leq\frac{1}{2}\sum\limits_{i=1}^{N}L_{i}\left(\left\lVert{r_{i}^{(k)}-r_{i}}\right\rVert^{2}-\left\lVert{r_{i}^{(k+1)}-r_{i}}\right\rVert^{2}\right)$
$\displaystyle\quad+\frac{N}{2\sigma^{2}}\left\lVert{z^{(k)}-z}\right\rVert^{2}-\frac{N}{2\sigma^{2}}\left\lVert{z^{(k+1)}-z}\right\rVert^{2}$
$\displaystyle\quad-\frac{N}{2\sigma^{2}}\left\lVert{z^{(k+1)}-z^{(k)}}\right\rVert^{2}.$

###### Proof.

Since $f(\mathbf{r})=\sum_{i=1}^{N}f_{i}(r_{i})$, we create a bound for each
$f_{i}$, and then sum over $i=1,\ldots,N$ to get the final result. First, we
apply Lemma 2 to $f_{i}$, mapping the variables $\xi$ to $r_{i}^{(k)},$ $x$ to
$r_{i}^{(k+1)},$ $y$ to $r_{i},$ and $u$ to the subgradient in Equation (47).
This gives

$\displaystyle f_{i}(r_{i}^{(k+1)})$ $\displaystyle-
f_{i}(r_{i})+\frac{1}{\sigma^{2}}\left\langle
u_{i}^{\left(k+\frac{1}{2}\right)},r_{i}^{(k+1)}-r_{i}\right\rangle$ (71)
$\displaystyle\leq\frac{1}{2}\left(L_{i}\left\lVert{r_{i}^{(k)}-r_{i}}\right\rVert^{2}-p_{i}\left\lVert{r_{i}^{(k+1)}-r_{i}}\right\rVert^{2}\right)$
(72)
$\displaystyle\leq\frac{L_{i}}{2}\left(\left\lVert{r_{i}^{(k)}-r_{i}}\right\rVert^{2}-\left\lVert{r_{i}^{(k+1)}-r_{i}}\right\rVert^{2}\right),$
(73)

where we have used the assumption $p_{i}\geq L_{i}$ from Definition 1 to
obtain the second inequality.

Next, we obtain a similar expression with $z$. First, note that by (46) and
(39), we have

$\frac{1}{N}\sum_{i=1}^{N}u_{i}^{\left(k+\frac{1}{2}\right)}=z^{(k+1)}-z^{(k)}.$
(74)

We then use (74) to rewrite the inner product of the subgradient in (47) with
$\left(z^{(k+1)}-z\right)$ as

$\displaystyle\frac{1}{\sigma^{2}}\sum\limits_{i=1}^{N}$
$\displaystyle\left\langle
u_{i}^{\left(k+\frac{1}{2}\right)},z^{(k+1)}-z\right\rangle$ (75)
$\displaystyle=\frac{N}{\sigma^{2}}\left\langle
z^{(k+1)}-z^{(k)},z^{(k+1)}-z\right\rangle$ (76)
$\displaystyle=\frac{N}{2\sigma^{2}}\Big{(}\left\lVert{z^{(k+1)}-z}\right\rVert^{2}-\left\lVert{z^{(k)}-z}\right\rVert^{2}$
(77)
$\displaystyle\quad+\left\lVert{z^{(k+1)}-z^{(k)}}\right\rVert^{2}\Big{)},$
(78)

where we have used the identity $2\left\langle
a-b,c-d\right\rangle=\left\lVert{a-d}\right\rVert^{2}-\left\lVert{a-c}\right\rVert^{2}-\left\lVert{b-d}\right\rVert^{2}+\left\lVert{b-c}\right\rVert^{2}$.

To prove the lemma, we sum the inequalities in (71)-(73) over $i=1,\ldots,N$
and subtract (75)-(78). ∎

######  Lemma 6.

Assume the hypotheses of Proposition 1. Then for any $\mathbf{u}$,

$\displaystyle\frac{1}{\sigma^{2}}\sum_{i=1}^{N}$ $\displaystyle\left\langle
u_{i}^{(k+\frac{1}{2})}-u_{i},r_{i}^{(k+1)}-z^{(k+1)}\right\rangle$ (79)
$\displaystyle=\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\Big{(}\left\lVert{u_{i}^{(k+1)}-u_{i}}\right\rVert^{2}-\left\lVert{u_{i}^{(k)}-u_{i}}\right\rVert^{2}$
(80) $\displaystyle\mbox{\hskip
110.96556pt}+\left\lVert{u_{i}^{(k+1)}-u_{i}^{(k)}}\right\rVert^{2}\Big{)}$
(81)

###### Proof.

From (40), we have $u_{i}^{(k+1)}-u_{i}^{(k)}=r_{i}^{(k+1)}-z^{(k+1)}$. Using
this right hand side directly plus solving for $u_{i}^{(k+1)}$ to simplify
(46), we can rewrite the inner products in the lemma statement as

$\displaystyle\Bigl{<}u_{i}^{(k+\frac{1}{2})}$ $\displaystyle-
u_{i},r_{i}^{(k+1)}-z^{(k+1)}\Bigr{>}$ (82) $\displaystyle=\left\langle
u_{i}^{(k+1)}-u_{i},u_{i}^{(k+1)}-u_{i}^{(k)}\right\rangle$ (83)
$\displaystyle\quad+\left\langle
z^{(k+1)}-z^{(k)},u_{i}^{(k+1)}-u_{i}^{(k)}\right\rangle.$ (84)

As noted just before (38), $\overline{{\mathbf{u}}}^{(k)}=0$ for $k\geq 2$, so
summing (84) over $i$ as in the lemma statement yields 0, hence this term may
be dropped. The term in (83) can be rewritten as

$\displaystyle\Bigl{<}u_{i}^{(k+1)}$ $\displaystyle-
u_{i},u_{i}^{(k+1)}-u_{i}^{(k)}\Bigr{>}$ (85)
$\displaystyle=\frac{1}{2}\Big{(}\left\lVert{u_{i}^{(k+1)}-u_{i}}\right\rVert^{2}-\left\lVert{u_{i}^{(k)}-u_{i}}\right\rVert^{2}$
(86)
$\displaystyle\quad+\left\lVert{u_{i}^{(k+1)}-u_{i}^{(k)}}\right\rVert^{2}\Big{)},$
(87)

which follows from the identity $2\left\langle
a-b,c-d\right\rangle=\left\lVert{a-d}\right\rVert^{2}-\left\lVert{a-c}\right\rVert^{2}-\left\lVert{b-d}\right\rVert^{2}+\left\lVert{b-c}\right\rVert^{2}$.
Summing over $i=1,\dots,N$ and multiplying by $\frac{1}{\sigma^{2}}$ completes
the proof. ∎

Now, we prove Lemma 3.

Lemma 3. Assume the hypotheses of Theorem 1, and that (41) has a KKT point,
$\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)$. Then the sequence
$\left\\{\left(\mathbf{r}^{(k)},z^{(k)},\mathbf{u}^{(k)}\right)\right\\}$
generated by the iterations (38)-(40) satisfies

$E^{(k+1)}\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)\leq
E^{(k)}\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right),$ (88)

where

$\displaystyle E^{(k)}$
$\displaystyle\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)=$ (89)
$\displaystyle\sum_{i=1}^{N}\Big{(}L_{i}\left\lVert{r_{i}^{(k)}-r_{i}^{*}}\right\rVert^{2}+\frac{1}{\sigma^{2}}\left\lVert{z^{(k)}-z^{*}}\right\rVert^{2}$
(90) $\displaystyle\mbox{\hskip
89.626pt}+\frac{1}{\sigma^{2}}\left\lVert{u_{i}^{(k)}-u_{i}^{*}}\right\rVert^{2}\Big{)}.$
(91)

###### Proof.

Since $\left(\mathbf{r}^{*},z^{*},\mathbf{u}^{*}\right)$ is a KKT point, it
satisfies the saddle point condition

$\mathcal{L}(\mathbf{r}^{*},z^{*},\mathbf{u})\leq\mathcal{L}(\mathbf{r}^{*},z^{*},\mathbf{u}^{*})\leq\mathcal{L}(\mathbf{r},z,\mathbf{u}^{*})$
(92)

for any $(\mathbf{r},z,\mathbf{u})$, where $\mathcal{L}$ is the Lagrangian
function in (42). Specifically at
$(\mathbf{r},z)=(\mathbf{r}^{(k+1)},z^{(k+1)})$, we have

$\displaystyle f(\mathbf{r}^{(k+1)})-f(\mathbf{r}^{*})$ (93)
$\displaystyle\quad+\frac{1}{\sigma^{2}}\sum_{i=1}^{N}\left\langle
u_{i}^{*},r_{i}^{(k+1)}-z^{(k+1)}\right\rangle\geq 0.$ (94)

We now use Lemma 5 and Lemma 6 to produce an upper-bound for this quantity.
First we consider the left side of the inequalities in each lemma. By
subtracting the left side of the inequality in Lemma 6 evaluated at
$\mathbf{u}=\mathbf{u}^{*}$ from the left side of the inequality in Lemma 5
evaluated at $\left(\mathbf{r},z\right)=\left(\mathbf{r}^{*},z^{*}\right)$, we
get

$\displaystyle f(\mathbf{r}^{(k+1)})-f(\mathbf{r}^{*})$ (95)
$\displaystyle\quad+\frac{1}{\sigma^{2}}\sum\limits_{i=1}^{N}\left\langle
u_{i}^{\left(k+\frac{1}{2}\right)},\left(r_{i}^{(k+1)}-r_{i}^{*}\right)-\left(z^{(k+1)}-z^{*}\right)\right\rangle$
(96) $\displaystyle\quad-\frac{1}{\sigma^{2}}\sum_{i=1}^{N}\left\langle
u_{i}^{(k+\frac{1}{2})}-u_{i}^{*},r_{i}^{(k+1)}-z^{(k+1)}\right\rangle.$ (97)

By the feasibility condition of the KKT point, $r_{i}^{*}=z^{*}$ for each $i$,
and hence (95)–(97) simplifies to the left side of (93). Now, by subtracting
the respective right sides of the inequalities in Lemma 5 and Lemma 6, and
combining with (93), we get

$\displaystyle 0$ $\displaystyle\leq f(\mathbf{r}^{(k+1)})-f(\mathbf{r}^{*})$
(98) $\displaystyle\quad+\frac{1}{\sigma^{2}}\sum_{i=1}^{N}\left\langle
u_{i}^{*},r_{i}^{(k+1)}-z^{(k+1)}\right\rangle$ (99)
$\displaystyle\leq\frac{1}{2}\sum\limits_{i=1}^{N}\Biggl{[}L_{i}\left(\left\lVert{r_{i}^{(k)}-r_{i}^{*}}\right\rVert^{2}-\left\lVert{r_{i}^{(k+1)}-r_{i}^{*}}\right\rVert^{2}\right)$
(100)
$\displaystyle\quad+\frac{1}{\sigma^{2}}\left(\left\lVert{z^{(k)}-z^{*}}\right\rVert^{2}-\left\lVert{z^{(k+1)}-z^{*}}\right\rVert^{2}\right)$
(101)
$\displaystyle\quad+\frac{1}{\sigma^{2}}\left(\left\lVert{u_{i}^{(k)}-u_{i}^{*}}\right\rVert^{2}-\left\lVert{u_{i}^{(k+1)}-u_{i}^{*}}\right\rVert^{2}\right)$
(102)
$\displaystyle\quad-\frac{1}{\sigma^{2}}\left(\left\lVert{z^{(k+1)}-z^{(k)}}\right\rVert^{2}+\left\lVert{u_{i}^{(k+1)}-u_{i}^{(k)}}\right\rVert^{2}\right)\Biggr{]}.$
(103)

Finally, dropping the final non-positive term and rearranging the remaining
terms gives

$\displaystyle\sum_{i=1}^{N}L_{i}\left\lVert{r_{i}^{(k)}-r_{i}^{*}}\right\rVert^{2}+\frac{1}{\sigma^{2}}\left\lVert{z^{(k)}-z^{*}}\right\rVert^{2}+\frac{1}{\sigma^{2}}\left\lVert{u_{i}^{(k)}-u_{i}^{*}}\right\rVert^{2}$
(104)
$\displaystyle\geq\sum_{i=1}^{N}L_{i}\left\lVert{r_{i}^{(k+1)}-r_{i}^{*}}\right\rVert^{2}+\frac{1}{\sigma^{2}}\left\lVert{z^{(k+1)}-z^{*}}\right\rVert^{2}$
(105)
$\displaystyle\quad+\frac{1}{\sigma^{2}}\left\lVert{u_{i}^{(k+1)}-u_{i}^{*}}\right\rVert^{2},$
(106)

which completes the proof. ∎

zotero ieeetr