# Image reconstruction algorithms in radio interferometry: from handcrafted to
learned regularization denoisers

Matthieu Terris,1 Arwa Dabbech,1 Chao Tang,1 Yves Wiaux1
1Institute of Sensors, Signals and Systems, Heriot-Watt University, Edinburgh
EH14 4AS, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We introduce a new class of iterative image reconstruction algorithms for
radio interferometry, at the interface of convex optimization and deep
learning, inspired by plug-and-play methods. The approach consists in learning
a prior image model by training a deep neural network (DNN) as a denoiser, and
substituting it for the handcrafted proximal regularization operator of an
optimization algorithm. The proposed AIRI (“AI for Regularization in radio-
interferometric Imaging”) framework, for imaging complex intensity structure
with diffuse and faint emission from visibility data, inherits the robustness
and interpretability of optimization, and the learning power and speed of
networks. Our approach relies on three steps. Firstly, we design a low dynamic
range training database from optical intensity images. Secondly, we train a
DNN denoiser at a noise level inferred from the signal-to-noise ratio of the
data. We use training losses enhanced with a nonexpansiveness term ensuring
algorithm convergence, and including on-the-fly database dynamic range
enhancement via exponentiation. Thirdly, we plug the learned denoiser into the
forward-backward optimization algorithm, resulting in a simple iterative
structure alternating a denoising step with a gradient-descent data-fidelity
step. We have validated AIRI against CLEAN, optimization algorithms of the
SARA family, and a DNN trained to reconstruct the image directly from
visibility data. Simulation results show that AIRI is competitive in imaging
quality with SARA and its unconstrained forward-backward-based version uSARA,
while providing significant acceleration. CLEAN remains faster but offers
lower quality. The end-to-end DNN offers further acceleration, but with far
lower quality than AIRI.

###### keywords:

techniques: image processing – techniques: interferometric

††pubyear: 2021††pagerange: Image reconstruction algorithms in radio
interferometry: from handcrafted to learned regularization denoisers–B

##  1 Introduction

Data acquisition by interferometry in radio astronomy relies on an array of
antennas measuring an incomplete coverage of the spatial Fourier domain of the
image of interest, yielding an ill-posed inverse problem towards image
formation. Modern and upcoming radio telescopes are designed to bring
unprecedented resolution and sensitivity. In this context, the algorithms to
be deployed to solve the inverse imaging problem face an ever increasing
requirement to jointly deliver precision (_i.e._ high resolution and dynamic
range), robustness (_i.e._ endowed with calibration and uncertainty
quantification functionalities), and scalability (_i.e._ the capability to
process sheer data volumes).

Most RI imaging pipelines utilized by astronomers are based on the CLEAN
algorithm, originally proposed by Högbom (1974). Assuming a sparse sky model,
this greedy algorithm iteratively removes the contribution of each point
source in the dirty image (more generally known as the backprojected image,
defined as the inverse Fourier transform of the data). Numerous extensions of
CLEAN have been devised over the last fifty years (_e.g._ Schwarz, 1978;
Schwab & Cotton, 1983; Bhatnagar & Cornwell, 2004; Cornwell et al., 2008;
Thompson et al., 2017). Albeit simple and computationally efficient, CLEAN-
based methods often require extensive manual intervention to ensure their
stability. Furthermore, CLEAN restored images often exhibit sub-optimal
imaging quality. On the one hand, their resolution is limited to that of the
observations due to the convolution with a smoothing beam. On the other hand,
their sensitivity is limited to the noise level, through the addition of the
residual image.

Since their first inception by Wiaux et al. (2009), convex and nonconvex
optimization algorithms underpinned by sparsity priors have emerged in radio
interferometric (RI) imaging (Li et al., 2011; Dabbech et al., 2015; Garsden
et al., 2015; Onose et al., 2016; Onose et al., 2017; Repetti et al., 2017;
Pratley et al., 2018; Dabbech et al., 2018; Repetti et al., 2018, 2019; Birdi
et al., 2020; Thouvenin et al., 2022a, b) and have led to a much higher image
precision than that offered by advanced versions of CLEAN. In this framework,
a so-called objective function is defined, typically as the sum of a data
fidelity term and a regularization term injecting a prior image model to
compensate for data incompleteness. The image estimate is defined as the
minimizer of this objective and is reached via provably convergent algorithms.
The obtained solution can also be understood in a Bayesian framework as a
_maximum a posteriori_ (MAP) estimate with respect to a posterior
distribution, the negative logarithm of which is the objective. Uncertainty
quantification approaches fully powered by optimization algorithms have also
been proposed (Repetti et al., 2018, 2019). The versatility of the framework
has progressively imposed optimization algorithms as a cornerstone of current
state-of-the-art techniques in RI imaging. We note that they have also proven
their worth for image reconstruction in the related fields of optical
interferometry (Thiébaut & Young, 2017) and very long baseline interferometry
(VLBI) (Akiyama et al., 2019, 2022). However, albeit robust and interpretable,
optimization algorithms come at the expense of a significant increase of the
computational cost over CLEAN, ultimately affecting their scalability to large
data volumes.

Bayesian inference approaches have also been proposed to address the RI image
formation problem (Cai et al., 2018; Arras et al., 2019), also for VLBI
(Akiyama et al., 2022), naturally enabling uncertainty quantification, but
remaining computationally very demanding for the data volumes expected by
modern radio telescopes.

Recent works have contemplated the use of deep neural networks (DNNs) for end-
to-end reconstruction or postprocessing in RI imaging (Terris et al., 2019;
Connor et al., 2022; Gheller & Vazza, 2021). Despite promising results and
tremendous scalability capabilities, these approaches have neither been
thoroughly validated for high resolution high dynamic range imaging of complex
structure involving diffuse and faint emission, nor compared with CLEAN or
state-of-the-art optimization approaches. Moreover, end-to-end DNNs are known
to be subject to robustness issues (Goodfellow et al., 2015; Nguyen et al.,
2015; Pang et al., 2018). Their exploitation raises questions, not only with
regards to the interpretability of the obtained solution, specifically in a
Bayesian context, but also their generalizability, _i.e._ their ability to
provide accurate reconstruction quality in the case of acquisition conditions
unseen during training.

In this work, we propose a new RI image reconstruction approach based on the
Plug-and-Play (PnP) framework, whereby a prior image model is learned by
training a DNN as a denoiser, and substituted for the so-called proximal
operator enforcing regularization at the heart of proximal optimization
algorithms (Venkatakrishnan et al., 2013; Chan et al., 2016; Zhang et al.,
2017, 2019). PnP algorithms were shown to deliver outstanding image
reconstruction quality in applications such as image restoration (Zhang et
al., 2021) and medical imaging (Ahmad et al., 2020). Importantly, as the DNN
is trained as a denoiser on an image database to serve as a simple
regularization operator, it is by construction applicable for any sensing
procedure, _de facto_ avoiding the generalizability issue affecting end-to-end
approaches. Leveraging monotone operator theory, recent works have shown that
these algorithms provably converge to well characterized solutions under
nonexpansiveness conditions on the denoiser, yielding similar interpretations
to that of MAP estimators (Cohen et al., 2021; Pesquet et al., 2021). In
summary, this versatile framework simultaneously inherits the robustness and
interpretability of optimization approaches, and the learning power and
scalability capabilities of DNNs.

Building on Pesquet et al. (2021), we propose a PnP framework dubbed AIRI,
standing for “AI for Regularization in radio-interferometric Imaging”.
Focusing on monochromatic intensity imaging, our approach for developing a
first AIRI algorithm capable of imaging complex structure with diffuse and
faint emission can be summarized in three steps. The first step is to design,
from publicly available optical images, a realistic but low dynamic range
database of intensity images for supervised training. The second step is to
train a DNN denoiser with basic architecture ensuring positivity of the
reconstructed images, at the signal-to-noise ratio of the data. The training
relies on a simple $\ell_{2}$ or $\ell_{1}$ loss enhanced with a firm
nonexpansiveness term ensuring algorithm convergence, and on an on-the-fly
exponentiation procedure for dynamic range enhancement. The third step resides
in plugging the DNN into the FB optimization algorithm, resulting in a simple
iterative structure alternating between a denoising step and a gradient-
descent data-fidelity step.

We finally validate the resulting AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$
algorithms, implemented in Matlab, on high resolution high dynamic range
simulations utilizing intensity images containing diffuse and faint emission
across the field of view. Our test images are 3c353, Hercules A, Centaurus A,
and Cygnus A, with size $512\times 512$. Our benchmark algorithms are (i) the
multi-scale CLEAN version implemented in the C++ WSClean software package
(Offringa & Smirnov, 2017), (ii) two optimization algorithms from the SARA
family leveraging the handcrafted “average sparsity” proximal regularization
operator: SARA itself (Onose et al., 2017) and its unconstrained FB-based
version uSARA that we introduce here, and (iii) a UNet trained in an end-to-
end fashion to reconstruct target images from dirty images. While SARA relies
on an advanced primal-dual FB (PDFB) algorithm (Pesquet & Repetti, 2015)
enabling non-differentiable data-fidelity constraints, uSARA only differs from
AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ by the use of the average sparsity
proximal operator in lieu of a learned denoiser for regularization. It in fact
corresponds to the standalone version of the imaging module of the joint
calibration and imaging approach proposed in Repetti et al. (2017) and Dabbech
et al. (2021). Results show that AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ are
already competitive with SARA and uSARA in terms of imaging quality, while
providing a significant acceleration. In a nutshell, AIRI does indeed inherit
the robustness and interpretability of optimization approaches, and the
learning power and speed of DNNs. The WSClean code remains significantly
faster but offers lower reconstruction quality. The UNet offers further
acceleration, but with far lower quality than AIRI.

The remainder of the paper is organized as follows. In Section 2, we review
imaging approaches from proximal optimization, recall SARA, and introduce its
unconstrained version uSARA. In Section 3, we review PnP algorithms through
the prism of proximal algorithms and introduce AIRI. In Section 4, we study
the performance of AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ through extensive
simulation in comparison with the benchmark algorithms. In Section 5, we draw
our conclusions and discuss necessary future work.

##  2 State-of-art optimization algorithms & RI

###  2.1 RI imaging problem

Aperture synthesis in radio astronomy probes the sky by measuring an
incomplete coverage of the spatial Fourier domain of the image of interest
through an array of antennas. Focusing on monochromatic intensity imaging,
assuming a narrow field of view, and in the absence of atmospheric and
instrumental perturbations, each pair of antennas acquires a noisy Fourier
component of the intensity image to be formed, called a visibility. The
associated Fourier mode (also called $uv$-point) is given by the projection of
the corresponding baseline, expressed in units of the observation wavelength,
onto the plane perpendicular to the line of sight (Thompson et al., 2017). A
discrete formulation of the resulting linear RI image formation problem,
aiming to restore a target intensity image
$\overline{\bm{x}}\in\mathbb{R}^{n}$ from the measured complex visibilities
$\bm{y}\in\mathbb{C}^{m}$, reads (Onose et al., 2016)

$\bm{y}=\bm{\Phi}\overline{\bm{x}}+\bm{e},$ (1)

where $\bm{\Phi}=\bm{\mathrm{GFZ}}\in\mathbb{C}^{m\times n}$ is the
measurement operator, $\bm{\mathrm{G}}\in\mathbb{C}^{m\times d}$ is a sparse
interpolation matrix, encoding the non-uniform Fourier transform,
$\bm{\mathrm{F}}\in\mathbb{C}^{d\times d}$ is the 2D Discrete Fourier
Transform, $\bm{\mathrm{Z}}\in\mathbb{R}^{d\times n}$ is a zero-padding
operator, incorporating the correction for the convolution performed through
the operator $\bm{\mathrm{G}}$, and $\bm{e}\in\mathbb{C}^{m}$ is a realization
of some i.i.d. Gaussian random noise, with zero mean and standard deviation
$\tau>0$. We refer to Appendix A for considerations regarding more general
noise distributions, and the inclusion of so-called direction-dependent
effects (DDEs) in $\bm{\mathrm{G}}$.

We also note that, backprojecting problem (1) into the image domain gives

$\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{y}\\}=\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}\overline{\bm{x}}+\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{e}\\},$
(2)

where $(\cdot)^{\dagger}$ denotes the complex conjugate transpose and
$\operatorname{Re}\\{\cdot\\}$ denotes the real part.
$\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{y}\\}$ is known as the dirty
image, $\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}$ is the operator
representing the convolution of $\overline{\bm{x}}$ by the point spread
function, also known as the dirty beam, and
$\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{e}\\}$ is the noise backprojected
in the image domain.

We finally emphasise that in the optical interferometry and VLBI contexts,
imaging is often performed from closure quantities derived from the
visibilities, bringing a nonlinear inverse problem affected by non-Gaussian
noise (Akiyama et al., 2019, 2022). Imaging in the presence of DDEs or from
closure quantities lies beyond the scope of the present work.

###  2.2 Proximal algorithms & FB

Leaving the nature of $\bm{\Phi}$ aside, problems of the likes of (1) are
ubiquitous in imaging sciences, and arise for instance in image restoration
(Levin et al., 2009; Yang et al., 2010; Bredies & Holler, 2020), hyperspectral
imaging (Wang et al., 2015; Xie et al., 2019) and medical imaging (Gupta et
al., 2018; Zbontar et al., 2018; Fessler, 2020), to name a few. A widespread
approach to solve such a problem is to reformulate it as a convex minimization
problem

$\underset{\bm{x}\in\mathbb{R}^{n}}{\text{minimize}}\,f(\bm{x})+\lambda
r(\bm{x}),$ (3)

where $f(\bm{x})+\lambda r(\bm{x})$ is called the objective function,
$f\in\Gamma_{0}(\mathbb{R}^{n})$111$\Gamma_{0}(\mathbb{R}^{n})$ denotes the
set of convex, proper and lower-semicontinuous functions from $\mathbb{R}^{n}$
to $(-\infty,\infty]$. is the term enforcing data-fidelity, and
$r\in\Gamma_{0}(\mathbb{R}^{n})$ is the regularization term introduced to
address the ill-posedness of (1) by enforcing a prior image model.

The theory of optimization offers a myriad of algorithms to solve such
minimization problems, ranging from the simple Forward-Backward (FB) algorithm
(Bauschke & Combettes, 2017), applicable if one of the two terms is
differentiable, to ADMM (Boyd et al., 2011), or the Douglas-Rachford algorithm
(Eckstein & Bertsekas, 1992), applicable even if both $f$ and $r$ are non-
differentiable, and more evolved structures, such as PDFB (Pesquet & Repetti,
2015), that can handle multi-term objective functions with as many non-
differentiable data-fidelity and regularization terms processed in parallel.
All these algorithms are part of the same class of proximal algorithms, where
differentiable and non-differentiable terms respectively translate in the
algorithm into gradient and proximal operators.

Considering $r\in\Gamma_{0}(\mathbb{R}^{n})$, its proximal operator
$\operatorname{prox}_{r}$ is defined as

$(\forall\bm{z}\in\mathbb{R}^{n}),\qquad\operatorname{prox}_{r}(\bm{z})=\underset{\bm{u}\in\mathbb{R}^{n}}{\text{argmin}}\frac{1}{2}\|\bm{z}-\bm{u}\|^{2}+r(\bm{u}),$
(4)

and can be interpreted as a generalization of a projection operator, or
interestingly, as a denoising operator. Indeed, problem (4) can be interpreted
as the minimization problem to be solved for a simple denoising problem of the
form $\bm{z}=\bm{u}+\bm{w}$ where the data $\bm{z}$ result from adding some
i.i.d. Gaussian noise $\bm{w}$ to the unknown $\bm{u}$. From this perspective,
the first term of the objective function in (4) would be the standard data-
fidelity term, given the Gaussian nature of the noise, and $r$ the
regularization term. A proximal operator is a denoiser!

The choice of $f$ and $r$ in (3) is of paramount importance as it influences
the final solution to the minimization task (Mallat, 1999; Selesnick et al.,
2005; Bredies & Holler, 2020). Interpreting the reconstructed image as a MAP
estimate leads naturally to choosing $f$ as the negative log-likelihood
associated with a given statistical model of the noise (typically a squared
$\ell_{2}$ norm for i.i.d. additive Gaussian noise), and $r$ as the negative
logarithm associated with a statistical prior image model (typically,
$\ell_{1}$ norm for a Laplace prior). However, optimization algorithms are not
tied to this statistical interpretation. In the context of compressive sensing
theory (Candès et al., 2006; Donoho, 2006; Baraniuk, 2007), where sparsity is
the postulated signal model, the choice of the $\ell_{1}$ norm as a
regularization term simply emanates from it being the closest convex
relaxation of the $\ell_{0}$ norm that is the natural sparsity measure.
$\ell_{1}$ regularization was shown to yield state-of-the-art image
reconstruction methods beyond the context of compressive sensing. For
instance, sparsity of the gradient is enforced as the $\ell_{1}$ norm of the
magnitude of the image gradient, that is the so-called Total Variation semi-
norm (TV) (Rudin et al., 1992; Bredies et al., 2010; Bredies & Holler, 2020),
whereas sparsity of the sought image in a sparsifying domain
$\bm{\mathrm{A}}\in\mathbb{R}^{n\times p}$ can be enforced by choosing
$r(\bm{x})=\|\bm{\mathrm{A}}^{\dagger}\bm{x}\|_{1}$ where $\bm{\mathrm{A}}$
can be a wavelet dictionary (Mallat, 1999), an x-let transform (Candes &
Demanet, 2003; Do & Vetterli, 2003), or a learned dictionary (Mairal et al.,
2009), to name a few. Generally, the choice of the data-fidelity term is often
driven by the statistical nature of the measurement noise, while
regularization terms are carefully handcrafted to meet the specificity of the
applications, with sparsity models defining the state-of-the-art. In summary,
the MAP interpretation has its limitation, and in particular, setting the
regularization parameter $\lambda$ according to the proper normalization of
the statistical models of which $f$ and $r$ would be the negative logarithms,
is known to be suboptimal.

The FB algorithm is a proximal algorithm designed to solve problems of the
form (3), where $f,r\in\Gamma_{0}(\mathbb{R}^{n})$ and $f$ is differentiable.
It is defined by the iterative sequence

$(\forall
k\in\mathbb{N}),\qquad\bm{x}_{k+1}=\operatorname{prox}_{\gamma\lambda
r}(\bm{x}_{k}-\gamma\nabla f(\bm{x}_{k})),$ (5)

which is proven to converge towards a minimizer of the objective function (3),
provided that $0<\gamma<2/L$ where $L$ is the Lipschitz constant222The
Lipschitz constant of an operator $h:\mathbb{R}^{n}\to\mathbb{R}^{p}$ is
defined as $\text{sup}_{x\neq y}\|h(x)-h(y)\|/\|x-y\|$. of $\nabla f$
(Bauschke & Combettes, 2017). Accelerated versions of FB, leveraging inertial
terms, preconditioning, and stochastic approaches, also exist. We also
emphasize that FB can be used to solve problems involving nonconvex
regularization terms $r$ (Beck & Teboulle, 2009; Attouch et al., 2013;
Chouzenoux et al., 2014; Combettes & Pesquet, 2015; Repetti & Wiaux, 2021).

###  2.3 SARA and its unconstrained version

####  2.3.1 SARA

The SARA family represent state-of-the-art optimization algorithms for RI
imaging (Carrillo et al., 2012, 2014; Abdulaziz et al., 2016; Onose et al.,
2016; Onose et al., 2017; Birdi et al., 2018; Pratley et al., 2018; Dabbech et
al., 2018; Abdulaziz et al., 2019; Thouvenin et al., 2022a). In its
monochromatic intensity imaging version, SARA consists in solving the
constrained minimization problem given by

$\underset{\bm{x}\in\mathbb{R}^{n}}{\text{minimize}}\,\,\iota_{\mathcal{B}(\bm{y},\epsilon)}(\bm{\Phi}\bm{x})+\lambda
r(\bm{x}),$ (6)

where data fidelity is imposed via a constraint encoded by the indicator
function333The indicator function of a non-empty, closed, and convex set
$\mathcal{S}$ is defined as $\iota_{\mathcal{S}}(\bm{x})=0$ if
$\bm{x}\in\mathcal{S}$, else $\iota_{\mathcal{S}}(\bm{x})=+\infty$.
$\iota_{\mathcal{B}(\bm{y},\epsilon)}$ of the $\ell_{2}$-ball centred at the
data $\bm{y}$ and of radius $\epsilon>0$:
$\mathcal{B}(\bm{y},\epsilon)=\left\\{\bm{u}\in\mathbb{C}^{m}\,|\,\|\bm{y}-\bm{u}\|\leq\epsilon\right\\}$.
The $\ell_{2}$-bound $\epsilon$ is derived from the assumed i.i.d. Gaussian
random noise statistics as $\epsilon^{2}=(2m+4\sqrt{m})\tau^{2}/2$ (Carrillo
et al., 2012). The function $r$ defines the SARA prior, also called the
“average sparsity” prior, which consists of a positivity constraint and a log-
sum prior promoting average sparsity in an over-complete dictionary
$\bm{\Psi}\in\mathbb{R}^{n\times bn}$. More specifically, the dictionary
$\bm{\Psi}$ is defined as the concatenation of $b=9$ orthogonal bases (the
first eight Daubechies wavelets and the Dirac basis), re-normalized by
$\sqrt{b}=3$ to ensure that $\bm{\Psi}\bm{\Psi}^{\dagger}=\mathbf{I}$, where
$\mathbf{I}$ denotes the identity operator. The SARA prior explicitly reads

$r(\bm{x})=\rho\sum_{j=1}^{bn}\log\left(\rho^{-1}\left|\left(\bm{\mathrm{\Psi}}^{\dagger}\bm{x}\right)_{j}\right|+1\right)+\iota_{\mathbb{R}^{n}_{+}}(\bm{x}),$
(7)

where $\left(.\right)_{j}$ denotes the $j^{\text{th}}$ coefficient of its
argument vector, and $\mathbb{R}^{n}_{+}$ denotes the $n$-dimensional real
positive orthant. The parameter $\lambda>0$ is a regularization parameter,
while $\rho>0$ is used to avoid reaching zero values in the argument of the
logarithmic terms. Thouvenin et al. (2022a) suggest to set both parameters to
an estimate of the standard deviation of the measurement noise in the wavelet
domain.

Minimizing (6) should typically involve the proximal operator
$\operatorname{prox}_{r}$, dubbed the “average sparsity” proximal operator.
However, if the log-sum prior enforces a stronger sparsity of the solution
than a typical $\ell_{1}$ prior, it is also nonconvex. To address the
resulting nonconvex minimization task, a reweighting procedure is adopted,
where a sequence of convex surrogate minimization problems is solved
iteratively (Candes et al., 2008; Carrillo et al., 2012), each involving a
weighted-$\ell_{1}$ prior $g$, given by

$g(\bm{x},\bm{\mathrm{W}})=\|\bm{\mathrm{W}}\bm{\mathrm{\Psi}}^{\dagger}\bm{x}\|_{1}+\iota_{\mathbb{R}^{n}_{+}}(\bm{x}),$
(8)

where ${\bm{\mathrm{W}}\in\mathbb{R}^{bn\times bn}}$ is a diagonal weighting
matrix that needs to be updated after each resolution of the surrogate
problem. The full SARA algorithm for solving (6) can be summarized as

$\begin{array}[]{l}\text{for}\;i=0,1,\ldots\\\
\left\lfloor\begin{array}[]{ll}\bm{\tilde{x}}_{i+1}&=\underset{\bm{x}\in\mathbb{R}^{n}}{\text{argmin}}\,\,\iota_{\mathcal{B}(\bm{y},\epsilon)}(\bm{\Phi}\bm{x})+g(\bm{x},\bm{\mathrm{W}}_{i})\\\
\bm{\mathrm{W}}_{i+1}&=\mathrm{Diag}\left(\rho\left/\left(\rho+\bm{\mathrm{\Psi}}^{\dagger}\bm{\tilde{x}}_{i+1}\right)\right.\right),\\\
\end{array}\right.\\\ \end{array}$ (9)

with $\bm{\mathrm{W}}_{0}=\mathbf{I}$, where $\text{Diag}(\cdot)$ denotes the
diagonal matrix containing its (vector) argument on the diagonal. Minimizing
each of these problems will now involve the weighted-$\ell_{1}$ prior $g$ via
its proximal operator $\operatorname{prox}_{g}$. One should acknowledge that
the reweighting strategy also introduces significant complexity to the
algorithm since a sequence of convex minimization problems needs to be solved.

Given the non-differentiability of the data-fidelity term in (6), the PDFB
algorithm (Pesquet & Repetti, 2015) is leveraged to solve each minimization
task of the form (9). We note that because $\bm{\Psi}^{\dagger}$ is over-
complete and due to the presence of the positivity constraint in (8), the
proximal operator $\operatorname{prox}_{g}$ does not admit a closed-form
solution, thus _a priori_ requiring a sub-iterative structure. However, the
full splitting functionalities of PDFB enable the decomposition of $g$ into
functions admitting simple proximal operators (a component-wise projection on
the real positive orthant and a component-wise thresholding operation). These
proximal operators are handled in parallel by PDFB without the need for sub-
iterations (Onose et al., 2016). Interestingly, the same functionalities
enable handling large data volumes via parallel processing of data blocks
(Onose et al., 2016), and large image sizes via parallel processing of image
facets (Thouvenin et al., 2022a), providing significant scalability to SARA.

####  2.3.2 Unconstrained SARA (uSARA)

Building on the works of Repetti & Wiaux (2020, 2021), we propose to focus on
the unconstrained formulation of the SARA problem, that we dub unconstrained
SARA (uSARA), again looked at from the prism of monochromatic intensity RI
imaging only. More precisely, we consider the problem (3), with a
differentiable data-fidelity term $f$ chosen as the mean squared error loss,
corresponding to the negative log-likelihood of the data under the i.i.d.
Gaussian random noise assumption, given by

$f(\bm{x})=\frac{1}{2}\|\bm{\Phi}\bm{x}-\bm{y}\|^{2}.$ (10)

Adopting the prior model $r$ described in (7), the resulting nonconvex
minimization task writes

$\underset{\bm{x}\in\mathbb{R}^{n}}{\text{minimize}}\,\frac{1}{2}\|\bm{\Phi}\bm{x}-\bm{y}\|^{2}+\lambda
r(\bm{x}),$ (11)

where $\lambda>0$ is a regularization parameter. The nonconvexity of the prior
$r$ from (7) is handled through a similar reweighting procedure as in (9):

$\begin{array}[]{l}\text{for}\;i=0,1,\ldots\\\
\left\lfloor\begin{array}[]{ll}\bm{\tilde{x}}_{i+1}&=\underset{\bm{x}\in\mathbb{R}^{n}}{\text{argmin}}\,\,\frac{1}{2}\|\bm{\Phi}\bm{x}-\bm{y}\|^{2}+\lambda
g(\bm{x},\bm{\mathrm{W}}_{i})\\\
\bm{\mathrm{W}}_{i+1}&=\mathrm{Diag}\left(\rho\left/\left(\rho+\bm{\mathrm{\Psi}}^{\dagger}\bm{\tilde{x}}_{i+1}\right)\right.\right),\\\
\end{array}\right.\\\ \end{array}$ (12)

with $\bm{\mathrm{W}}_{0}=\mathbf{I}$. At each reweighting iteration, problem
(12) can be solved with a standard FB algorithm. The full reweighting
procedure for solving the uSARA problem (11) is summarized in Algorithm 1,
where $\xi_{1}>0$ is a relative variation convergence criterion, and where the
gradient of $f$ at $\bm{x}$ reads as

$\nabla
f(\bm{x})=\text{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}\bm{x}-\text{Re}\\{\bm{\Phi}^{\dagger}\bm{y}\\},$
(13)

and its Lipschitz constant is the spectral norm444The spectral norm of a
linear operator is its maximum singular value, which can be computed via the
power method (see _e.g._ Golub & Van Loan (2013)). of
$\text{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}$, _i.e._

$L=\|\text{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}\|_{\rm{S}}.$ (14)

We recall that the stepsize $\gamma$ in (5) is to be upper-bounded by $2/L$.

As appears in Algorithm 1 solving (12) (steps 5-7), the weighted-$\ell_{1}$
prior $g$ is involved via $\operatorname{prox}_{\gamma\lambda g}$. The
algorithmic structure does not exhibit the same full splitting functionalities
as PDFB, and $\operatorname{prox}_{\gamma\lambda g}$ must be computed
iteratively. This can typically be achieved via a so-called dual FB algorithm
(Combettes & Pesquet, 2011), detailed in Algorithm 2, where $\xi_{2}>0$ is a
relative variation convergence criterion. This algorithm alternates between a
projection on the real positive orthant (Step 3), denoted by
$\Pi_{\mathbb{R}^{n}_{+}}$, and the proximal operator of the dual of the
weighted-$\ell_{1}$ prior, which involves a component-wise soft-thresholding
operator
$\operatorname{prox}_{\gamma\lambda\|\bm{\mathrm{W}}\cdot\|_{1}}$555This
operator, applied to $\bm{z}$ boils down component-wise to
$(\operatorname{prox}_{\gamma\lambda\|\bm{\mathrm{W}}\cdot\|_{1}}(\bm{z}))_{j}=\operatorname{prox}_{\eta_{j}|\cdot|}(z_{j})=\operatorname{sign}(z_{j})\operatorname{max}\\{|z_{j}|-\eta_{j},0\\}$
where $\eta_{j}=\gamma\lambda\mathrm{W}_{j,j}$ is called the soft-thresholding
parameter. All values below the threshold are set to $0$ while those above are
reduced by the value of the threshold (in absolute value). (Step 4).

Formally, Algorithms 1 and 2 lead to a triply sub-iterative structure: a
proximal operator loop (Step 2), inside a weighted-$\ell_{1}$ loop (Step 5),
inside a reweighting loop (Step 3). This suggests a very computationally
expensive structure. Nonetheless, firstly, Repetti & Wiaux (2021) have shown
that only a fixed number of FB iterations $K$ are required (Step 2),
corresponding to the approximate minimization of (12), while preserving the
convergence of the overall algorithm to a minimizer of the nonconvex objective
(11). This contrasts with SARA, where the sub-problems (9) need to be solved
to convergence. The value of $K$ can in fact be optimized to provide
significant acceleration of Algorithm 1, effectively removing one iteration
layer. Secondly, the number of iterations in Algorithm 2 is moderate in
practice when using an appropriate initialization strategy for the dual
variable $\bm{v}_{0}$, leading to mild computational cost of
$\operatorname{prox}_{\gamma\lambda g}$.

We conclude this section by underlining that uSARA corresponds to the imaging
module of the joint calibration imaging approach described in Repetti et al.
(2017) and Dabbech et al. (2021).

Algorithm 1 Re-weighted FB algorithm for uSARA

1:given $0<\gamma<2/L$, $\lambda>0$, $\rho>0$,
$\bm{\tilde{x}}_{0}\in\mathbb{R}^{n}$, $\xi_{1}>0$

2:Set $\bm{\mathrm{W}}_{0}=\mathbf{I}$

3:repeat for $i=0,1,\ldots$

4: $\bm{x}_{0}=\bm{\tilde{x}}_{i}$

5: repeat for $k=0,\ldots,K$

6: $\bm{x}_{k+1}=\operatorname{prox}_{\gamma\lambda
g(\cdot,\bm{\mathrm{W}}_{i})}(\bm{x}_{k}-\gamma\nabla f(\bm{x}_{k}))$

7: end for

8: ${\bm{\tilde{x}}}_{i+1}=\bm{x}_{K}$

9:
$\bm{\mathrm{W}}_{i+1}=\mathrm{Diag}\left(\rho\left/\left(\rho+\bm{\mathrm{\Psi}}^{\dagger}\bm{\tilde{x}}_{i+1}\right)\right.\right)$

10:until $\|x_{k+1}-x_{k}\|/\|x_{k+1}\|<\xi_{1}$

11:return ${\bm{\tilde{x}}}_{i+1}$

Algorithm 2 Dual FB algorithm for computing
$\operatorname{prox}_{\gamma\lambda g(\cdot,\bm{\mathrm{W}})}(\bm{z})$ in
Algorithm 1

1:given $\gamma>0$, $\lambda>0$, $\bm{\mathrm{W}}$, $\bm{z}\in\mathbb{R}^{n}$,
$\bm{v}_{0}\in\text{Span}(\bm{\mathrm{\Psi}}^{\dagger}$), $\xi_{2}>0$

2:repeat for $l=0,1,\ldots$

3:
$\bm{x}_{l+1}=\Pi_{\mathbb{R}_{+}^{n}}(\bm{z}-\bm{\mathrm{\Psi}}\bm{v}_{l})$

4:
$\bm{v}_{l+1}=\left(\mathbf{I}-\operatorname{prox}_{\gamma\lambda\|\bm{\mathrm{W}}\cdot\|_{1}}\right)\left(\bm{v}_{l}+\bm{\mathrm{\Psi}}^{\dagger}\bm{x}_{l+1}\right)$

5:until $\|x_{k+1}-x_{k}\|/\|x_{k+1}\|<\xi_{2}$

6:return $\bm{x}_{l+1}$

##  3 AIRI: AI for Regularization in RI Imaging

###  3.1 PnP-FB

When solving an inverse imaging problem from an optimization theory viewpoint,
one defines an objective function, of which the sought image would be a
minimizer and which can be obtained via an iterative algorithm. The PnP
approach (Venkatakrishnan et al., 2013) looks at the problem directly through
the lens of the algorithm. As discussed in Section 2.2, it follows from the
definition of a proximal operator that it can be interpreted as a denoiser.
Borrowing a proximal optimization algorithm, the PnP approach proposes to
replace the proximal regularization operator with a more general denoiser. We
here note that this procedure can be applied to a wide class of algorithms,
ranging from FB to PDFB. Specifically, assuming a minimization problem of the
form (3) with the differentiable data fidelity term $f$, recall that the FB
algorithm reads as in (5), with the proximal operator (4). Its PnP counterpart
(PnP-FB) then simply follows as

$(\forall
k\in\mathbb{N}),\qquad\bm{x}_{k+1}=\operatorname{D}(\bm{x}_{k}-\gamma\nabla
f(\bm{x}_{k})),$ (15)

where $\operatorname{D}$ is a denoising operator, _i.e._ an operator
specifically designed to remove i.i.d. Gaussian random noise from an image.

Denoisers are ubiquitous in imaging sciences, and a large variety of denoisers
have been defined and studied in the signal processing literature, including
Gaussian filters, BM3D (Dabov et al., 2007), non local means (NLM) (Buades et
al., 2011), with, more recently, deep neural networks (DNNs) leading the state
of the art (Zhang et al., 2017). It is therefore very tempting to plug
powerful denoisers in lieu of proximal operators… and play. In this context,
we aim at taking advantage of the learning capabilities of DNNs in order to
learn an appropriate regularization denoiser. Just as handcrafting the
regularization term $r$ (from which the proximal regularization operator
$\operatorname{prox}_{r}$ results) is crucial to the reconstruction quality
for a pure optimization approach, learning a powerful denoiser is cornerstone
for the reconstruction quality of PnP algorithms. Recent advances in DNNs for
image denoising tasks have shown new state of the art results over traditional
denoisers (Zhang et al., 2017; Wang et al., 2018; Zhou et al., 2020).
Simultaneously, PnP algorithms with DNNs as denoisers have become the new
state-of-the-art in image reconstruction (Zhang et al., 2021), significantly
improving over their traditional optimization counterparts.

Yet, replacing the proximal operator with an off-the-shelf denoising operator
is not inconsequential. In particular, a general non-proximal denoiser can _a
priori_ not be related with a regularization term in some overarching
objective function. As a consequence, the characterization of the PnP solution
(_e.g._ as a MAP estimator), and worse, the convergence properties of the
iterative structure, are not ensured anymore, thus questioning the robustness
and interpretability of PnP solutions. To overcome that issue, we follow the
approach of Pesquet et al. (2021) allowing to enforce the so-called “firm
nonexpansiveness constraint” on the denoiser $\operatorname{D}$ during
training. This constraint ensures that $\operatorname{D}$ “contracts
distances” and yields both the convergence of Algorithm 3 and the
characterization of the limit point (see Appendix B for more details). The
PnP-FB algorithm with established convergence guarantees is summarized in
Algorithm 3, where $\xi_{3}>0$ is a relative variation convergence criterion.
In the RI imaging case of interest, the PnP counterpart to uSARA follows from
defining $f$ as in (10), with the measurement operator $\bm{\Phi}$ from (1).
The gradient of $f$ at $\bm{x}$ therefore reads as in (13), with Lipschitz
constant given in (14).

Algorithm 3 AIRI (PnP-FB) algorithm

1:given $0<\gamma<2/L$, denoiser $\operatorname{D}$,
$\bm{x}_{0}\in\mathbb{R}^{n}$, $\xi_{3}>0$

2:repeat for $k=0,1,\ldots$

3: $\bm{x}_{k+1}=\operatorname{D}(\bm{x}_{k}-\gamma\nabla f(\bm{x}_{k}))$

4:until $\|x_{k+1}-x_{k}\|/\|x_{k+1}\|<\xi_{3}$

5:return $x_{k+1}$

###  3.2 Training noise level vs. uSARA regularization parameter

####  3.2.1 Denoising i.i.d. Gaussian noise of specific standard deviation

Building from Pesquet et al. (2021), our AIRI denoiser $\operatorname{D}$ will
be trained in a supervised manner to remove i.i.d. Gaussian random noise of
specific standard deviation $\sigma$. In other words, $\operatorname{D}$ is
set to tackle the denoising problem

$\bm{z}=\bm{u}+\sigma\bm{w},$ (16)

where $\bm{z}$ denotes the noisy data input to the DNN, $\bm{u}$ is the
unknown image corresponding to the target output of the DNN, and
$\bm{w}\sim\mathcal{N}(0,\mathbf{I})$. As discussed in the next section, the
input signal-to-noise ratio of the RI data provides a reliable handle on
$\sigma$.

####  3.2.2 Equating training & target dynamic ranges

As already emphasized, AIRI will emerge from PnP-FB with $f$ from (10) and
$\bm{\Phi}$ from (1). We propose to set the training noise level $\sigma$ for
$\operatorname{D}$ to an appropriate estimate of the standard deviation of
some effective image-domain noise induced by the original i.i.d. Gaussian
random noise of standard deviation $\tau$ on the Fourier data in (1). We
propose below a procedure to estimate the standard deviation of this image-
domain noise, which, for simplicity, we also approximate to be i.i.d. Gaussian
random noise.

Firstly, we resort to the image-domain formulation (2) of (1), normalized by
the Lipschitz constant $L$ in (14):
$\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{y}\\}/L=\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}\overline{\bm{x}}/L+\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{e}\\}/L$.
By construction, the chosen normalization ensures that the convolution
operator $\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}/L$ has unit
spectral norm, so that the dirty image is at the same scale as the original
image. The covariance matrix of the image-domain noise
$\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{e}\\}/L$ thus writes
$\tau^{2}\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}/2L^{2}$.
Secondly, as proposed in Thouvenin et al. (2022a), we discard the correlation
structure of the image-domain noise via the approximation
$\operatorname{Re}\\{\bm{\Phi}^{\dagger}\bm{\Phi}\\}\simeq L\mathbf{I}$,
leading to a noise covariance matrix approximation $\tau^{2}\mathbf{I}/2L$. As
a result the effective image-domain noise is assumed i.i.d. Gaussian, and the
original standard deviation $\tau$ in the measurement domain is rescaled to
$\tau/\sqrt{2L}$ in the transfer to the image domain. Equating the standard
deviations of the training noise for $\operatorname{D}$ and the image-domain
noise leads to the following heuristic:

$\sigma=\frac{\tau}{\sqrt{2L}}.$ (17)

We note that this heuristic is algorithm agnostic, _i.e. a priori_ independent
of the algorithmic structure in which the denoiser is plugged. In this sense,
while similar, it is simpler and more general than the one proposed in Pesquet
et al. (2021), which was derived specifically by estimating $\sigma$ around
the fixed point of (15). As will be shown through simulations in Section 4,
the proposed heuristic provides a very precise estimate of the optimal
training noise level.

Assuming the validity of this heuristic, one must acknowledge a dependency of
the training noise level $\sigma$ on the measurement model, via $\tau$ and
$L$, suggesting that a denoiser should be trained independently for each
target reconstruction. Also, we consider a normalized database of images whose
maximum intensity values are lower or equal to 1 (see Section 3.3), so that
the denoising performance cannot be ascertained on non-normalized images.
These two generalizability questions can be solved by a simple rescaling of
the original inverse problem (1) to simultaneously fit the normalization
constraint (an upper bound on the peak intensity of the image being accessible
from the dirty image) and rescale the noise level in the data to the training
noise level of some already available denoiser. Dabbech et al. (2022) and
Wilber et al. (2022a) build from these considerations and propose both a
“single denoiser” approach and a more flexible “denoiser shelf” approach,
where denoisers are respectively trained at a single predefined noise level or
on a limited number of such noise levels, to which the target inverse problems
can be rescaled at will. It also appears from the above rescaling
considerations that the dimension of $\sigma$ is that of an inverse input
image-domain peak signal-to-noise ratio rather than an absolute noise level.
As the AIRI denoiser consistently categorizes $\sigma$-level values as noise
to be removed, the input image-domain signal-to-noise ratio also naturally
sets the target output peak signal-to-noise ratio or dynamic range of the
reconstruction. In summary, our heuristic simply stipulates that the training
dynamic range should be set equal to the target reconstruction dynamic range.

####  3.2.3 uSARA regularization parameter

We further note the very similar roles of the noise level $\sigma$ of the
regularization denoiser in PnP-FB and the regularization parameter $\lambda$
in the FB algorithm (5) solving (3). At the algorithmic level, on the one
hand, $\operatorname{D}$ is a denoiser and we have proposed to set $\sigma$ to
the estimate (17) of the effective image-domain noise level, or more generally
to the inverse target reconstruction dynamic range value. On the other hand,
for the traditional FB algorithm, $\lambda$ is involved in (5) via
$\operatorname{prox}_{\gamma\lambda r}$. Focusing on the specific uSARA
instance described in Section 2.3.2, and summarized on Algorithms 1 and 2,
$\gamma\lambda$ typically acts as a soft-thresholding parameter, thresholding
out small values in the wavelet domain defined by the average sparsity
dictionary. Setting its value to the standard deviation of the effective
wavelet-domain noise, which according to Thouvenin et al. (2022a) is the same
as in the image domain given the unit spectral norm of the average sparsity
dictionary $\bm{\Psi}$, yields the heuristic:

$\gamma\lambda=\frac{\tau}{\sqrt{2L}}.$ (18)

We readily note that, in contrast with the reasoning leading to (17), the
above reasoning for $\lambda$ is specific to the FB structure underpinning
uSARA. We also recall the parameter $\rho>0$ in the definition of the average
sparsity prior $r$ in (7), used to avoid reaching zero values in the argument
of the logarithmic terms. This parameter represents a floor level for the
wavelet coefficients (Carrillo et al., 2012; Thouvenin et al., 2022a), and
should naturally be set to the same effective wavelet-domain noise level
estimate as $\lambda$: $\rho=\lambda$. We further note that a variation of the
proposed heuristic for the wavelet-domain noise level arises from a slightly
different approximation in the transfer of the noise level from the image
domain to the wavelet domain. Instead of assuming noise level conservation due
to the normalization of $\bm{\Psi}$, one can analyze the noise in each of the
$9$ bases of the dictionary separately, discarding the noise correlation
structure between each pair of bases. In this case, the orthonormality of each
basis and its normalization by $3$ suggests a correction factor $1/3$ in the
heuristic values: $\rho=\gamma\lambda\simeq\tau/3\sqrt{2L}$.

As will be shown in Section 4, relation (18) provides a useful reference to
set $\lambda$ in uSARA. Even though the $1/3$ correction will be shown to
bring a closer-to-optimal value, it is still not as precise an estimate as
(17) is for the training noise level of AIRI denoisers. We emphasize that a
significant body of work has been dedicated to the question of setting the
regularization parameter in problems of the form (3), with no simple and
accurate solution to our knowledge (Donoho & Johnstone, 1995; Luisier et al.,
2007; Vidal et al., 2020). In this context, the fact that (18) provides an
appropriate reference value for $\lambda$ in uSARA is an interesting result.
However, the fact that the PnP-FB admits a very accurate heuristic to set the
training noise level, avoiding the requirement for fine-tuning an arbitrary
regularization parameter, represents a major advantage of AIRI over uSARA.

###  3.3 Training database with adaptive dynamic range

####  3.3.1 State-of-the-art and proposed approach

Supervised training of DNNs for imaging requires the definition of a database
of images, which remains a main difficulty in translating deep learning
methods to RI imaging due to the absence of groundtruths. The few works that
have leveraged deep learning techniques for RI imaging so far resorted to
various approaches. In Terris et al. (2019), we used a high dynamic range
simulated radio map from Bonaldi & Braun (2018) of size $32000\times 32000$
containing mainly elliptical Gaussian sources, broken down into a database of
1100 images of size $512\times 512$. Gheller & Vazza (2021) relied on 1000
radio images of size $2000\times 2000$ simulated with ENZO (Bryan et al.,
2014) and augmented by random rotations, resulting in the appearance of
realistic extended emission in the field of view. Connor et al. (2022)
generated a database of 900 synthetic images of size $2000\times 2000$
containing elliptical Gaussian sources.

In this work, we are targeting high dynamic range imaging of complex structure
with diffuse and faint emission across the field of view. The database to be
designed for the training of the regularization denoiser thus needs to reflect
this target. We propose (i) to source publicly available optical astronomy
images containing a rich combination of compact and diffuse emission, (ii) to
preprocess these images, aiming to remove existing noise and artefacts via a
preliminary denoising procedure, and (iii) to exponentiate their intensity
values to create high dynamic range images. The procedure results in the
creation of a rich and clean training database of $2235$ training images of
size $512\times 512$. More technical details are provided below.

####  3.3.2 Building a rich low dynamic range database

Firstly, we constituted a set of 32 large grayscale images
$(\bm{u}^{\text{raw}}_{k})_{1\leq k\leq 32}$ from optical astronomical
intensity images available online, containing complex structure with both
compact and diffuse emission. Each image is normalized with maximum intensity
value equal to $1$, and when the image contained colour channels, those were
averaged. The typical size of each of these images is $3500\times 3500$. The
sample standard deviation of the residual noise in image regions with no
signal, averaged over the 32 images reads $\sigma_{0}\simeq 1/64$, leading to
a typical dynamic range of $\sigma_{0}^{-1}\simeq 64$, orders of magnitude
below the values of interest in modern RI observations.

Secondly, the images are preprocessed with the aim of removing the residual
noise, which includes various compression artefacts, and would otherwise be
learned as part of the image model during training. To that effect, for all
$1\leq k\leq 32$, a simple denoising problem is formulated for the recovery of
a clean image $\overline{\bm{u}}_{k}$ from $\bm{u}_{k}^{\text{raw}}$:
$\bm{u}_{k}^{\text{raw}}=\overline{\bm{u}}_{k}+\bm{e}_{k}$, assuming i.i.d.
Gaussian random noise $\bm{e}_{k}$ with standard deviation
$\widehat{\sigma}_{k}$. The denoising procedure applied consists in solving
the optimization problem underpinning SARA in (6), with $\bm{\Phi}=\mathbf{I}$
and $\bm{y}=\bm{u}_{k}^{\text{raw}}$, resulting in a denoised database
$(\bm{u}^{\text{low}}_{k})_{1\leq k\leq 32}$. We note that this procedure does
not affect the (low) dynamic range of the images. In a nutshell, noise and
artefacts are removed, but the floor signal level and maximum intensity remain
essentially unchanged.

Finally, the cleaned images are split to create a rich and clean database
$(\bm{u}^{\text{low}}_{s})_{1\leq s\leq S}$ with $S=2235$ images of size
$512\times 512$ containing both compact and diffuse emission. When the image
dimensions are not multiples of $512$, a symmetric padding is used before
splitting. The resulting images exhibit a whole distribution of maximum
intensity values bounded by $1$, and a corresponding distribution of dynamic
range values, with $\sigma_{0}^{-1}\simeq 64$ representing the “nominal”
dynamic range value of the database.

####  3.3.3 Enhancing database dynamic range to target

In order to simulate the high dynamic range of interest in RI imaging, we
propose to exponentiate pixel-wise the intensity value of all groundtruth
images of the low dynamic range database. The exponentiation parameter will be
set to achieve a final nominal dynamic range of the same order as the training
dynamic range, which according to (17) is also equivalent to the target
dynamic range of the reconstruction. Technically, the low dynamic range images
$\bm{u}^{\text{low}}_{s}$ are exponentiated through pixel-wise as:

$\bm{u}_{s}=h(\bm{u}^{\text{low}}_{s},a)=\frac{a^{\bm{u}^{\text{low}}_{s}}-1}{a},$
(19)

for some parameter $a\gg 1$. Recalling that the maximum span of the pixel
values in any $\bm{u}^{\text{low}}_{s}$ is $[\sigma_{0},1]$, the resulting
maximum span of the exponentiated images is approximately
$[a^{-1}(a^{\sigma_{0}}-1),1]$. The exponentiation thus preserves the database
normalization with maximum intensity value across images bounded by $1$, while
providing control on the final nominal dynamic range
$a(a^{\sigma_{0}}-1)^{-1}$. Imposing this value to be equal to the training
dynamic range leads to the following equation for $a$:

$a=\left(1+a\sigma\right)^{\sigma_{0}^{-1}}.$ (20)

Parametrizing $a$ as $a=b\sigma^{-1}$ leads to the following equation for the
scaling factor $b$: $b=\sigma\left(1+b\right)^{\sigma_{0}^{-1}}$. This
equation can be solved numerically to estimate the exponentiation factor $a$
as a function of $\sigma_{0}$ and $\sigma$. In summary, the high dynamic range
groundtruth database $(\bm{u}_{s})_{1\leq s\leq S}$ results from the low
dynamic range groundtruth database $(\bm{u}^{\text{low}}_{s})_{1\leq s\leq S}$
via (19), with exponentiation parameter $a$ set in (20).

Figure 1 summarizes the preprocessing pipeline for the creation of the low and
high dynamic range databases. Notice the diversity of images in the training
samples displayed in Figure 1 (c).

Figure 1: Summary of the proposed strategy for building a training database of
postprocessed radio groundtruths. We first gather a set of 32 astronomical
(optical) images of high quality but low dynamic range: (a) shows a subset of
four images. (b) illustrates the effect of preprocessing on a small part of
the image from (a). Eventually, we split the images into $2235$ images of size
$512\times 512$ to create the low dynamic range training database. Samples of
our training database, after application of the exponentiation (19) with
$a=10^{3}$, are shown in (c) in logarithmic scale.

(a) 3c353

(b) Hercules A

(c) Centaurus A

(d) Cygnus A

Figure 2: The four test RI images, displayed in logarithmic scale.

####  3.3.4 Generating noisy & groundtruth pairs on-the-fly

The supervised training approach considered relies on image pairs: the
groundtruth images $(\bm{u}_{s})_{1\leq s\leq S}$ corresponding to the target
outputs of the DNN, and the noisy images $(\bm{z}_{s})_{1\leq s\leq S}$
corresponding to the DNN input. Given our choice, described in Section 3.2.1,
to build a denoiser tailored to tackle i.i.d. Gaussian random noise of
specific standard deviation $\sigma$, the noisy high dynamic range images read

$\bm{z}_{s}=\bm{u}_{s}+\sigma\bm{w}_{s},$ (21)

with $\bm{w}_{s}\sim\mathcal{N}(0,\bm{\mathrm{I}})$ and $\sigma$ given by
(17).

Importantly, as the high dynamic range database of noisy and groundtruth pairs
depends on the target dynamic range through (19) and (21), one version of the
database would need to be computed and stored for each dynamic range of
interest, which is a highly impractical and unnecessary approach. Instead,
$\bm{z}_{s}$ and $\bm{u}_{s}$ can simply be seen as functions of
$\bm{u}^{\text{low}}_{s}$ and computed on the fly during training.

###  3.4 Training loss function

In this section, we introduce the training loss for the supervised training of
the denoiser $\operatorname{D}$. Building on Pesquet et al. (2021) and
following technical arguments detailed in Appendix B, we regularize the
training loss with an appropriate Jacobian term to enforce the firm
nonexpansiveness of $\operatorname{D}$. By denoting
$\operatorname{Q}=2\operatorname{D}-\operatorname{I}$ where $\operatorname{I}$
denotes the identity operator, and denoting $\bm{\theta}\in\mathbb{R}^{c}$ the
learnable parameters of $\operatorname{D}$, the resulting training loss reads:

$\displaystyle\underset{\bm{\theta}\in\mathbb{R}^{c}}{\text{minimize}}\,\frac{1}{S}\sum_{s=1}^{S}\bigg{(}\mathcal{L}(\operatorname{D}_{\bm{\theta}}(\bm{z}_{s})-\bm{u}_{s})+\kappa\operatorname{max}\\{\|\bm{\nabla}\operatorname{Q}_{\bm{\theta}}(\bm{z}_{s})\|_{\rm{S}},1-\varepsilon\\}\bigg{)},$
(22)

with $\bm{z}_{s}$ and $\bm{u}_{s}$ given in (19) and (21) respectively, as
functions of $\bm{u}^{\text{low}}_{s}$ computed on the fly during training.
The second term in (22) is the firm nonexpansiveness regularization term, with
$\kappa>0$ the associated regularization parameter, and $\varepsilon>0$ a
safety margin parameter. The first term $\mathcal{L}$ is the standard part of
the loss function, for which practical choices involve the $\ell_{2}$ loss
_i.e._ $\mathcal{L}(\cdot)=1/2\|\cdot\|^{2}$, and the $\ell_{1}$ loss, _i.e._
$\mathcal{L}(\cdot)=\|\cdot\|_{1}$; both will be studied in Section 4.

Figure 3: Architecture of the proposed denoising DNN, inspired from Zhang et
al. (2017) and Pesquet et al. (2021), containing 20 convolutional layers in
total. “Conv” denotes convolutional layers. The number of channels is given
above the arrows. Notice the ReLU layer at the output, ensuring that our
reconstructed RI images are positive.

##  4 Simulations and results

###  4.1 Overall simulation set up

####  4.1.1 Test images

Our test image database consists of four radio images, namely 3c353, Hercules
A, Centaurus A and Cygnus A, all resized to a size of $n=512\times 512$, shown
in Figure 2. These images all have a peak value equal to 1 and dynamic ranges
above $10^{4}$. The simulations involving these images will target
reconstruction dynamic ranges of the order of $10^{4}$. We emphasise that this
is orders of magnitude higher than typical values in the PnP literature (Zhang
et al., 2017; Zhang et al., 2021; Pesquet et al., 2021), and therefore
significantly more challenging.

####  4.1.2 Benchmark methods

We benchmark AIRI with (i) uSARA (Algorithm 1) solving the unconstrained
problem (11), (ii) SARA (Onose et al., 2017; Thouvenin et al., 2022a), solving
the constrained problem (6), (iii) the multiscale variant of CLEAN implemented
in the WSClean software (Offringa & Smirnov, 2017), and (iv) a UNet (Jin et
al., 2017) trained end-to-end on a large database to recover groundtruths from
a database of dirty images.

All simulations involving uSARA rely on a value for $\rho$ using (18)
including the $1/3$ correction factor, while the values of $\lambda$ are fine-
tuned around (18) to optimize the reconstruction quality. As for SARA, we rely
on the prescriptions of Thouvenin et al. (2022a) for setting its
hyperparameters. While all other methods use natural weighting of the data,
CLEAN makes use of uniform weighting, which relies on the density of the
Fourier sampling for enhanced resolution (see Section 2.1 and Appendix
A)666The specific command used is: wsclean -multiscale -niter 30000 -weight
uniform -gain 0.1 -mgain 0.6 -threshold 0.001 -no-reorder -minuvw-m 0.001
-auto-mask 0.001 -mem 1.. We underline that the CLEAN restored image is in
units of Jansky per beam area, in contrast with SARA-based methods, AIRI, and
the UNet, where the associated estimated images are in the units of Jansky per
pixel. Therefore, in what follows, CLEAN images are rescaled by the area of
the CLEAN beam, for visualization purposes. Finally, we note that we always
show the dirty image computed with a uniform weighting scheme, and rescaled in
the $[0,1]$ range for visualization purposes.

Regarding the end-to-end deep learning approach, we use a Deep Residual UNet
from Zhang et al. (2021) and train it to reconstruct a target image in (1)
from its dirty version in (2). We generate a database
$(\bm{u}_{s},\operatorname{Re}\\{\bm{\Phi}_{s}^{\dagger}\bm{y}_{s}\\})_{1\leq
s\leq S}$, where for every $s$, $\bm{y}_{s}$ represents the visibility vector
computed from the groundtruth $\bm{u}_{s}$ and for a measurement operator
$\bm{\Phi}_{s}=\bm{\mathrm{G}}_{s}\bm{\mathrm{FZ}}$ generated from a random
pointing direction in the sky. The randomization ensures that the UNet
training database contains a rich and diverse set of sampling patterns. We
borrow the procedure detailed in Section 3.3.3 to enhance the nominal dynamic
range of the groundtruth database $(\bm{u}_{s})_{1\leq s\leq S}$, and equate
it to the target dynamic range of reconstruction. The UNet is then trained to
minimize the $\ell_{1}$ loss between the groundtruth $\bm{u}_{s}$ and the
reconstruction
$\operatorname{UNet}(\operatorname{Re}\\{\bm{\Phi}_{s}^{\dagger}\bm{y}_{s}\\})$.
The UNet is trained for 800 epochs with the Adam algorithm.

####  4.1.3 Evaluation metrics

For AIRI-$\ell_{2}$, AIRI-$\ell_{1}$, uSARA, SARA, and the UNet, the image
estimates $\widehat{\bm{x}}$ are the reconstructed model images, while the
image estimate considered for CLEAN is the restored image. The performance of
all algorithms is assessed qualitatively via inspection of $\widehat{\bm{x}}$,
alongside the residual images, defined as
$\widehat{\bm{r}}=\beta\bm{\Phi}^{\dagger}(\bm{y}-\bm{\Phi}\widehat{\bm{x}})$,
where $\beta$ is a normalization factor777The normalization factor is given by
$\beta=1/\max_{i}(\bm{\Phi}^{\dagger}\bm{\Phi}{\bm{\delta}})_{i}$, where
$\bm{\delta}$ is the image with value 1 at the phase centre and 0 otherwise.
Hence, the dirty beam given by
$\beta\bm{\Phi}^{\dagger}\bm{\Phi}{\bm{\delta}}$, has a peak value equal to
1.. Given the high dynamic ranges of interest, the estimated images are shown
in logarithmic scale, using the mapping
$\operatorname{rlog}=\bm{x}\mapsto\operatorname{log}_{10}(10^{3}\bm{x}+1)/3$.
The residual images are displayed in linear scale.

We also consider two quantitative evaluation metrics, namely the SNR and
logSNR of the estimated images $\widehat{\bm{x}}$ with respect to the
groundtruth $\overline{\bm{x}}$. The standard SNR metric (in dB) is defined as
$\text{SNR}(\widehat{\bm{x}},\overline{\bm{x}})=20\operatorname{log}_{10}(\|\overline{\bm{x}}\|/\|\overline{\bm{x}}-\widehat{\bm{x}}\|)$.
It simply corresponds to a logarithmic version of the norm of the error image
$\overline{\bm{x}}-\widehat{\bm{x}}$. The logSNR (in dB) is the SNR evaluated
on a logarithmic version of the intensity images. This metric is introduced
given the high dynamic ranges of interest, and assigns more relative weight to
low intensities than the SNR. Our definition is
$\text{logSNR}(\widehat{\bm{x}},\overline{\bm{x}})=\text{SNR}(\operatorname{rlog}(\widehat{\bm{x}}),\operatorname{rlog}(\overline{\bm{x}}))$.

####  4.1.4 DNN architecture & training particulars

Firstly, as highlighted already, the denoising losses that we use for
$\mathcal{L}$ in our training loss (22) are chosen to be either the $\ell_{2}$
loss or the $\ell_{1}$ loss, giving rise to AIRI-$\ell_{2}$ and
AIRI-$\ell_{1}$.

Secondly, we choose a basic DNN architecture for the denoiser
$\operatorname{D}$, in the form of a modified DnCNN (Zhang et al., 2017) where
the batch normalization layers have been removed, as in Pesquet et al. (2021).
But in this work, a ReLU has been added as the last layer. This is a simple
way to enforce positivity at the output of the denoiser, and therefore of the
final reconstruction, much necessary as the developed AIRI algorithm targets
intensity (versus polarization) imaging (see Figure 3).

Thirdly, the DNN is trained for $10^{5}$ epochs with the Adam algorithm
(Kingma & Ba, 2015) and learning rate $10^{-4}$, the learning rate being
divided by 2 every $2\times 10^{4}$ epochs. $(\bm{z}_{s},\bm{u}_{s})$ in (22)
are of size $46\times 46$, generated according to (21) and (19) from patches
randomly extracted from the low dynamic range database of $512\times 512$
images $(\bm{u}^{\text{low}}_{s})_{1\leq s\leq S}$, itself augmented with
random rotations and zooms. The Jacobian regularization in (22) is computed
with a power method with a number of iterations increased to 10 in the later
stages of the training to improve the precision of the computation. As
observed in Pesquet et al. (2021), starting the optimization of the network
from a pretrained state with $\kappa=0$ in (22) enhances the results
significantly.

Finally, we have observed that performing, at each iteration of the PnP
algorithm, random flips and $90$ degrees rotations of the image before
applying the DNN, and inverting the transform after the denoising step, can
significantly boost the reconstruction quality. We therefore adopt this
strategy, inspired by Zhang et al. (2021).

###  4.2 Experiment 1: validating the denoiser for denoising

####  4.2.1 Simulation setup

The PnP approach suggests that the best denoisers should provide the best
regularizers. In this first experiment, we compare the denoising capabilities
of the AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ denoisers on one side, and the
average sparsity proximal operator propelling uSARA and SARA on the other
side. To this end, we focus on a denoising problem

$\bm{y}=\overline{\bm{x}}+\bm{e},$ (23)

where $\bm{y}\in\mathbb{R}^{n}$ is the data,
$\overline{\bm{x}}\in\mathbb{R}^{n}$ is the groundtruth image, and
$\bm{e}\in\mathbb{R}^{n}$ is the realization of an i.i.d. Gaussian random
noise with zero mean and standard deviation $\tau$. In other words, (23) is a
degenerate case of problem (1) with $\bm{\Phi}=\textbf{I}$. We consider the
3c353 image as a groundtruth $\overline{\bm{x}}$. As the following experiments
for the full validation of AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ operate at a
target dynamic range around $10^{4}$, we design this preliminary pure
denoising experiment with a similarly challenging target dynamic range, and
set the observation noise level to $\tau=10^{-4}$.

Firstly, problem (23) can be solved using AIRI via Algorithm 3. As the image
peaks at $1$, no rescaling of the inverse problem in (1) is needed (see
Section 3.2). In this simple case, the gradient in (13) reads $\nabla
f(\bm{x})=\bm{x}-\bm{y}$ with Lipschitz constant in (14) $L=1$. Taking
$\gamma=1/L=1$, Algorithm 3 trivially boils down to a simple non-iterative
application of $\operatorname{D}$ to the data, leading to the denoised image

$\widehat{\bm{x}}=\operatorname{D}(\bm{y}).$ (24)

The AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ denoisers were trained with a noise
level $\sigma=\tau=10^{-4}$ as per (17), on the proposed synthetic training
database, and following the procedure explained in Section 3.4. Following the
procedure detailed in Section 3.3.3, the chosen exponentiation parameter for
on-the-fly dynamic range enhancement of the database is $a=10^{3}$.
Furthermore, we set the regularization parameter in the training loss (22) to
$\kappa=10^{-9}$ and $\kappa=10^{-5}$ for AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$
denoisers respectively. The safety margin parameter is set to
$\varepsilon=5\times 10^{-2}$. These values are chosen here as they are the
ones to ensure the stability of Algorithm 3 in our next experiments.

Secondly, problem (23) can also be solved in a proximal approach with uSARA,
setting $\bm{\Phi}=\textbf{I}$. Interestingly, in this case, the uSARA
objective boils down to the objective defining the average sparsity proximal
operator (see definition (4)) of $\lambda r$, for the average sparsity prior
$r$ in (7), and $\lambda=\tau$ following the reasoning around equation (18).
Therefore, applying Algorithm 1 results in computing the denoised image by
simply applying the average sparsity proximal operator with appropriate
$\lambda$:

$\widehat{\bm{x}}=\operatorname{prox}_{\lambda r}(\bm{y}).$ (25)

We note, given $\nabla f(\bm{x})=\bm{x}-\bm{y}$ and $L=1$, and choosing $K=1$
and $\gamma=1/L=1$, that Algorithm 1 simplifies to the sequential application
of $\operatorname{prox}_{\gamma\lambda g}(\bm{y},\bm{\mathrm{W}}_{i})$ with
only the weights being updated at each iteration. The algorithm was run with
convergence criteria to $\xi_{1}=6\times 10^{-6}$ Algorithm 1, and
$\xi_{2}=10^{-4}$ in Algorithm 2.

####  4.2.2 Experimental results

Denoising results are displayed in Figure 4. All three denoisers appear to be
highly effective, both in terms of SNR and logSNR and visual denoising
quality. The SNR values are very comparable. In terms of logSNR the
AIRI-$\ell_{1}$ denoiser achieves better performance than the AIRI-$\ell_{2}$
denoiser, itself superior to the average sparsity proximal operator. Visually,
the AIRI-$\ell_{1}$ denoiser also achieves the reconstruction with the best
resolution and least amount of artefacts (see zooms).

These results suggests that the AIRI-$\ell_{1}$ and AIRI-$\ell_{2}$ denoisers
do encapsulate a superior prior model to the average sparsity model
encapsulated in the corresponding proximal operator.

(a) Groundtruth

(b) Noisy image
$(47.6\,\text{dB},24.1\,\text{dB})$

(c) Average sparsity prox. operator
$(53.8\,\text{dB},33.1\,\text{dB})$

(d) AIRI-$\ell_{2}$ denoiser
$(53.7\,\text{dB},35.1\,\text{dB})$

(e) AIRI-$\ell_{1}$ denoiser
$(53.9\,\text{dB},35.8\,\text{dB})$

Figure 4: Experiment 1 results: Denoising results obtained with the proposed
AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ denoisers in comparison with the averaged
sparsity proximal operator. The respective groundtruth and noisy images are
displayed in 4(a) and 4(b). Denoising results obtained with the averaged
sparsity proximal operator are shown in 4(c), with AIRI-$\ell_{2}$ in 4(d) and
with AIRI-$\ell_{1}$ in 4(e). All images are displayed in logarithmic scale.
The colourbar is saturated at $10^{-1}$ for clarity. Values of the obtained
evaluation metrics (SNR, logSNR) are indicated below their associated images.

###  4.3 Experiment 2: validating the training noise level heuristic

####  4.3.1 Simulation setup

Using our four $512\times 512$ test images (see Figure 2), we study the impact
of the noise level $\sigma$ involved in the training of the AIRI-$\ell_{2}$
and AIRI-$\ell_{1}$ denoisers according to the procedure set in Section 3,
with the aim to validate the heuristic (17). We also study the optimal value
of the regularization parameter $\lambda$ of uSARA around the heuristic (18).
The performance of AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ in terms of
reconstruction quality is compared to that of uSARA.

We simulate RI observations utilizing simulated $uv$-coverages of the radio
telescope MeerKAT. The $uv$-patterns correspond to five randomly selected
pointing directions displayed in Figures 5 and 6(c), and two of the four total
observation durations displayed in Figure 6: $\Delta
T\in\\{4\,\text{h},8\,\text{h}\\}$. The observed wavelength is $0.3\text{m}$
and the simulated sampling rate of the telescope is fixed to 100 points per
hour. The data size thus increases linearly with $\Delta T$. Data are
simulated following the model described in (1), with input SNR defined as
$\text{iSNR}=20\log_{10}\left(\|\bm{\Phi}\overline{\mathbf{x}}\|/\tau\right)=30~{}\text{dB}$.
We finally note that the $uv$-patterns extend to the edge of the considered
Fourier domain as shown in Figures 5 and 6. In other words, we target mild, if
any, super-resolution with respect to the nominal resolution of the
observation, as set by the largest baseline.

(a)

(b)

(c)

(d)

Figure 5: Experiments 1-3: MeerKAT $uv$-coverages simulated for four different
telescope pointing directions and fixed total observation duration $\Delta
T=4\,\text{h}$. The components are measured in units of the wavelength of the
observation. Each $uv$-pattern contains $m=806400$ points.

(a) $\Delta T=1\,\text{h}$

(b) $\Delta T=2\,\text{h}$

(c) $\Delta T=4\,\text{h}$

(d) $\Delta T=8\,\text{h}$

Figure 6: Experiments 1-3: MeerKAT $uv$-coverages associated with the four
observation durations and fixed telescope pointing direction (complementary
the four shown in Figure 5). The numbers of points in the $uv$-patterns are
(a) $\Delta T=1\text{h}:201600$, (b) $\Delta T=2\text{h}:403200$, (c) $\Delta
T=4\,\text{h}:806400$, and (d) $\Delta T=8\,\text{h}:1612800$.

Since our groundtruth images have peak values at $1$, no rescaling of the
inverse problem in (1) is needed for AIRI, and a common heuristic value is
reached for $\sigma$ in (17) and $\gamma\lambda$ in (18). Obviously, $\tau$
and $L$ are functions of $\bm{\Phi}$, and therefore of both the pointing
direction and total duration of observation. However, the heuristic values are
within a $25\%$ variation range across pointing directions for each $\Delta
T$. For simplicity, a single heuristic value is considered for all five
$uv$-patterns associated with the same $\Delta T$, defined as $\nu_{\Delta
T}$, taken to be the average of the values resulting for each pointing
direction. The two resulting values are as announced associated with a target
dynamic range around $10^{4}$: $\nu_{4}=1.4\times 10^{-4}$ for $\Delta
T=4\,\text{h}$, and $\nu_{8}=9.3\times 10^{-5}$ for $\Delta T=8\,\text{h}$.
The values probed for $\gamma\lambda$ and $\sigma$ are sampled a factor
$\sqrt{2}$ around the heuristic and by steps of $2$ thereafter.

Following the study by Pesquet et al. (2021), we choose the value of
$\gamma=1.98/L$ in Algorithms 1 and 3, which saturates the theoretical bound.
The convergence criteria of Algorithms 1 and 3 are set to
$\xi_{1}=\xi_{3}=5\times 10^{-6}$, with a maximum number of iterations of
$6\times 10^{3}$, and we set $\xi_{2}=10^{-5}$ in Algorithm 2.

We train our AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ denoisers with optimized
values $\kappa=10^{-9}$ and $\kappa=10^{-5}$ respectively in the training loss
(22), with $\varepsilon=5\times 10^{-2}$, and observe that this ensures the
stability of Algorithm 3 despite a spectral norm of the Jacobian slightly
above $1$. Following the procedure detailed in Section 3.3.3, the
exponentiation parameter for on-the-fly dynamic range enhancement of the
database is set to $a=10^{3}$.

####  4.3.2 Experimental results

Considering the image of 3c353 as groundtruth and the $uv$-pattern from Figure
5(a) ($\Delta T=4\,\text{h}$) for simulating the measurements, reconstructions
for different values of $\sigma$ in AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$, and
$\lambda$ in uSARA are displayed in Figures 7 and 8. More precisely, model
images, displayed in logarithmic scale, and their associated SNR and logSNR
values are provided in Figure 7, and residual images are displayed in Figure 8
in linear scale. We firstly notice that $\sigma$ and $\lambda$ play a similar
role in the reconstruction quality: the higher their value, the smoother the
recovered image; the lower their value, the stronger the artefacts. Overall,
reconstructions with uSARA exhibit some wavelet artefacts, particularly
noticeable around bright compact sources overlaying faint extended emission.
AIRI-$\ell_{2}$ reconstructions typically exhibit less artefacts, but appear
to be generally smoother. AIRI-$\ell_{1}$ reconstructions contain even less
artefacts, with an achieved resolution similar to uSARA. We conclude these
observations by highlighting the low-amplitude residual images obtained with
low values of $\sigma$ and $\lambda$. However, such results are not a token of
a good quality reconstruction. They instead reflect the over-fitting of the
data in such settings.

(a) Groundtruth

(b) Dirty image

(c) uSARA, $\gamma\lambda=\nu_{4}/4\sqrt{2}$
$(25.71\,\text{dB},23.47\,\text{dB})$

(d) uSARA, $\gamma\lambda=\nu_{4}/2\sqrt{2}$
$(24.31\,\text{dB},25.12\,\text{dB})$

(e) uSARA, $\gamma\lambda=\nu_{4}/\sqrt{2}$
$(22.77\,\text{dB},24.48\,\text{dB})$

(f) uSARA, $\gamma\lambda=2\sqrt{2}\nu_{4}$
$(18.90\,\text{dB},21.54\,\text{dB})$

(g) AIRI-$\ell_{2}$, $\sigma=\nu_{4}/4\sqrt{2}$
$(25.14\,\text{dB},19.03\,\text{dB})$

(h) AIRI-$\ell_{2}$, $\sigma=\nu_{4}/2\sqrt{2}$
$(26.05\,\text{dB},24.03\,\text{dB})$

(i) AIRI-$\ell_{2}$, $\sigma=\nu_{4}$
$(25.93\,\text{dB},24.62\,\text{dB})$

(j) AIRI-$\ell_{2}$, $\sigma=2\sqrt{2}\nu_{4}$
$(24.57\,\text{dB},21.09\,\text{dB})$

(k) AIRI-$\ell_{1}$, $\sigma=\nu_{4}/4\sqrt{2}$
$(26.64\,\text{dB},21.27\,\text{dB})$

(l) AIRI-$\ell_{1}$, $\sigma=\nu_{4}/2\sqrt{2}$
$(28.31\,\text{dB},26.02\,\text{dB})$

(m) AIRI-$\ell_{1}$, $\sigma=\nu_{4}$
$(28.85\,\text{dB},26.83\,\text{dB})$

(n) AIRI-$\ell_{1}$, $\sigma=2\sqrt{2}\nu_{4}$
$(27.13\,\text{dB},24.86\,\text{dB})$

Figure 7: Experiment 2 results: Impact of the training noise level $\sigma$
(resp. regularization parameter $\lambda$) in imaging with AIRI (resp. uSARA)
on simulated RI data using 3c353 image as a groundtruth and the $uv$-pattern
from Figure 5(a) ($\Delta T=4\,\text{h}$) in comparison with the reference
value $\nu_{4}$ suggested by the AIRI and uSARA heuristics. Top row: the
groundtruth image 7(a) and the simulated dirty image 7(b). Second row:
estimated model images obtained with uSARA, with a regularization parameter
$\lambda$ increasing with the column index. Third (resp. fourth) row:
estimated model images of AIRI-$\ell_{2}$ (resp. AIRI-$\ell_{1}$), with a
training noise level $\sigma$ increasing with the column index. Below each
image we indicate the reconstruction metrics as (SNR, logSNR). Images are
shown in logarithmic scale.

(a) uSARA, $\gamma\lambda=\nu_{4}/4\sqrt{2}$

(b) uSARA, $\gamma\lambda=\nu_{4}/2\sqrt{2}$

(c) uSARA, $\gamma\lambda=\nu_{4}/\sqrt{2}$

(d) uSARA, $\gamma\lambda=2\sqrt{2}\nu_{4}$

(e) AIRI-$\ell_{2}$, $\gamma\lambda=\nu_{4}/4\sqrt{2}$

(f) AIRI-$\ell_{2}$, $\sigma=\nu_{4}/2\sqrt{2}$

(g) AIRI-$\ell_{2}$, $\sigma=\nu_{4}$

(h) AIRI-$\ell_{2}$, $\sigma=2\sqrt{2}\nu_{4}$

(i) AIRI-$\ell_{1}$, $\gamma\lambda=\nu_{4}/4\sqrt{2}$

(j) AIRI-$\ell_{1}$, $\sigma=\nu_{4}/2\sqrt{2}$

(k) AIRI-$\ell_{1}$, $\sigma=\nu_{4}$

(l) AIRI-$\ell_{1}$, $\sigma=2\sqrt{2}\nu_{4}$

Figure 8: Experiment 2 results: Residual images associated with the
reconstructions of Figure 7. Top row: residual images for the estimated model
images obtained with uSARA, with a regularization parameter $\lambda$
increasing with the column index. Second (resp. third) row: respective
residual images for the estimated model images obtained with AIRI-$\ell_{2}$
(resp. AIRI-$\ell_{1}$), with a training noise level $\sigma$ increasing with
the column index. All images are shown in linear scale.

Figure 9 shows the reconstruction SNR and logSNR graphs obtained when varying
$\sigma$ and $\gamma\lambda$ for reconstructions with AIRI and uSARA
respectively, for both the total observation durations $\Delta T=4\,\text{h}$
and $\Delta T=8\,\text{h}$. Each point on each of the four graphs gives the
average value and $95\%$ confidence interval over 20 reconstructions arising
when considering the four test images in Figure 2 and five different pointing
directions from Figures 5 and 6(c).

Firstly, for both AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$, we observe that the SNR
and logSNR metrics clearly peak at the proposed heuristic (17) (shown as the
dashed vertical black line) for both $\Delta T=4\,\text{h}$ and $\Delta
T=8\,\text{h}$, suggesting that equating the training dynamic range to the
target dynamic range is an accurate procedure. This result is also confirmed
by the visual analysis in Figure 7. For uSARA, the SNR and logSNR results and
the visual analysis in Figure 7 suggest that the relation (18) is a useful
reference, but rather provides an upper-bound on the optimal value for the
soft-thresholding parameter, roughly 6 to 11 times smaller when looked at from
an SNR standpoint, or 3 to 6 times smaller from a logSNR perspective. In
general, we find that the logSNR provides better coherence with the visual
analysis, suggesting a optimal soft-thresholding value 3 to 6 times below the
reference. Interestingly, the $1/3$ correction factor brings a closer-to-
optimal value, but does not solve the discrepancy between the SNR and logSNR
metrics. In contrast to AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$, uSARA thus
exhibits a regularization parameter of which the value cannot be set
simultaneously automatically and optimally. This represents a significant
advantage of AIRI over uSARA.

Secondly, considering each method at its best, SNR and logSNR values for both
$\Delta T=4\,\text{h}$ and $\Delta T=8\,\text{h}$ confirm a superior
performance of AIRI-$\ell_{1}$ over both AIRI-$\ell_{2}$ and uSARA, the latter
two achieving similar reconstruction qualities, with virtually same logSNR,
and a slight SNR advantage for AIRI-$\ell_{2}$.

Thirdly, these results are in line with the visual analysis when considering
each method at its best (see Figures 7(d), 7(i), 7(m), and Figures 8(b), 8(g),
8(k)). AIRI-$\ell_{2}$ and uSARA exhibit similar reconstruction qualities with
the former, smoother but with less artefacts than the latter. AIRI-$\ell_{1}$
exhibits less artefacts than AIRI-$\ell_{2}$, while preserving the resolution
offered by uSARA. AIRI-$\ell_{1}$ provides the best residuals, uSARA residuals
clearly contain sources discarded by the model and relating to the artefacts,
while AIRI-$\ell_{2}$’s higher residuals correlate with the sub-optimal
resolution.

Fourthly, we recall that AIRI-$\ell_{2}$, AIRI-$\ell_{1}$, and uSARA are
underpinned by the same FB structure, and only differ by their choice of
denoiser. From this perspective, the results of the current experiment are in
line with those of the previous experiment, with a better regularization
denoiser naturally leading to better overall reconstruction.

Finally, we observe that the confidence intervals around the mean SNR and
logSNR values in Figure 9 are rather small for both AIRI-$\ell_{2}$ and
AIRI-$\ell_{1}$, and similar to those of uSARA. This is a first illustration
that in the PnP approach there is no issue of generalizability to measurement
conditions not considered during training as, by construction, the denoisers
are blind to the measurement setting (with regards to the noise level, see
discussion around (17)), in contrast to end-to-end networks.

(a) Case $\Delta T=4\,\text{h}$

(b) Case $\Delta T=8\,\text{h}$

(c) Case $\Delta T=4\,\text{h}$

(d) Case $\Delta T=8\,\text{h}$

Figure 9: Experiment 2 results: Reconstruction metrics as a function of the
training noise level $\sigma$ in the case of AIRI-$\ell_{2}$ and
AIRI-$\ell_{1}$, or the thresholding parameter $\gamma\lambda$ in the case of
uSARA. Top row: reconstruction SNR in the case of $\Delta T=4\,\text{h}$ (a)
and $\Delta T=8\,\text{h}$ (b). Bottom row: reconstruction logSNR in the case
of $\Delta T=4\,\text{h}$ (c) and $\Delta T=8\,\text{h}$ (d). Each point is an
average over the 20 simulated observations built from the four different
groundtruths from Figure 2 and the five different $uv$-patterns shown in
Figures 5 and 6(c). Error bars show the 95% confidence interval. On each
graph, the black vertical dotted line indicates the common heuristic value for
$\sigma$ and $\gamma\lambda$ (_i.e._ $\nu_{4}$ for $\Delta T=4\,\text{h}$ and
$\nu_{8}$ for $\Delta T=8\,\text{h}$). Values for $\sigma$ and $\gamma\lambda$
are sampled a factor $\sqrt{2}$ around the heuristic and by steps of $2$
thereafter.

###  4.4 Experiment 3: validating AIRI in precision & cost

####  4.4.1 Simulation setup

In this section, we validate AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ more
extensively, both in reconstruction quality and computational cost, across a
wider range of observation durations, and in comparison with, uSARA, SARA,
CLEAN, as well as UNets trained in an end-to-end fashion. We consider the same
four groundtruth images and experimental setup as for the previous experiment
(see Section 4.3.1), with the same set of five pointing directions (see
$uv$-coverages in Figures 5 and 6(c)). In order to further assess robustness
to varying measurement conditions, the set of observation durations is
enlarged to $\Delta T\in\\{1\,\text{h},2\,\text{h},4\,\text{h},8\,\text{h}\\}$
(see Figure 6). We note that a different UNet is trained for each observation
duration to alleviate the challenge of generalizability with respect to the
sampling pattern. The exponentiation parameter for the UNet training database,
resulting from the procedure in Section 3.3.3, is $a=10^{3}$.

We recall that our groundtruth images have peak values at $1$, yielding the
same heuristic value for $\sigma$ and $\gamma\lambda$. These values are still
within a $25\%$ variation range across pointing directions for all four
$\Delta T$ values, and a single heuristic value is considered for all five
$uv$-patterns associated with the same $\Delta T$. The resulting four average
$\nu_{\Delta T}$ are: $\nu_{1}=2.9\times 10^{-4}$, $\nu_{2}=2.2\times
10^{-4}$, $\nu_{4}=1.4\times 10^{-4}$, and $\nu_{8}=9.3\times 10^{-5}$. Given
the result of the previous experiment the values of $\sigma$ for
AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ are taken exactly at the heuristic, while
the value of $\gamma\lambda$ for uSARA is fine-tuned manually to optimize the
logSNR of each reconstruction, again resulting in soft-thresholding values
ranging between $3$ and $6$ times below the heuristic, depending on $\Delta
T$. Moreover, setting $\kappa=10^{-9}$ (resp. $\kappa=10^{-5}$) for the
AIRI-$\ell_{2}$ (resp. AIRI-$\ell_{1}$) denoiser in (22), with
$\varepsilon=5\times 10^{-2}$, ensured the stability of Algorithm 3. Following
the procedure detailed in Section 3.3.3, the exponentiation parameter for on-
the-fly dynamic range enhancement of the database is $a=10^{3}$.

At the level of Algorithms 1, 2, and 3, the stepsize $\gamma$, convergence
criteria $\left(\xi_{i}\right)_{1\leq i\leq 3}$ and maximum number of
iterations are set to the same values as in the previous experiment.

Finally, AIRI-$\ell_{2}$, AIRI-$\ell_{1}$, uSARA, and SARA are run in Matlab,
while CLEAN is implemented in the highly optimized C++ WSClean package, and
the UNets are implemented in Python using PyTorch. uSARA, SARA, and CLEAN
implementations utilize 10 CPU cores (Intel Xeon E5-2695 2.1 GHz). The
parallel PDFB algorithmic structure underpinning SARA (see Section 2.3.1)
optimizes the distribution of the operations involved at each iteration across
those 10 CPU cores: 1 for the data-fidelity term, and 1 for each of the 9
bases of the average sparsity dictionary $\bm{\Psi}$ (Onose et al., 2016). For
uSARA, the average sparsity dictionary is parallelized via a faceting
procedure using 9 facets (Pruša, 2012; Thouvenin et al., 2022a).
AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ utilize 10 CPU cores and 1 GPU (NVIDIA
Tesla V100). The latter is for the application of the denoiser, as commonly
done with DNNs in deep learning. The DNN can be run on CPU, but this slows
down the reconstructions significantly, as will be evident from the
experimental results. Implementations of uSARA, SARA, and CLEAN leveraging
GPUs are not investigated here. The UNets are run utilizing 1 GPU (NVIDIA
Tesla V100).

####  4.4.2 Reconstruction quality

Figure 10 shows the average reconstruction SNR and logSNR as a function of the
observation duration $\Delta T$. Each point, on each of the two graphs,
represents an average over 20 reconstructions arising when considering the
four test images in Figure 2 and five pointing directions, while error bars
give the $95\%$ confidence intervals. As expected, the reconstruction SNR
increases with $\Delta T$, as the data size increases with $\Delta T$. The
results extend the conclusion of the previous experiment of a superior
performance of AIRI-$\ell_{1}$ over both AIRI-$\ell_{2}$ and uSARA across
$\Delta T\in\\{1\,\text{h},2\,\text{h},4\,\text{h},8\,\text{h}\\}$. In SNR,
AIRI-$\ell_{1}$ and AIRI-$\ell_{2}$ are respectively around 3dB and 2dB above
uSARA. In logSNR, AIRI-$\ell_{1}$ is between 1dB and 2dB above both
AIRI-$\ell_{2}$ and uSARA. The comparison with SARA reveals that
AIRI-$\ell_{1}$ is on par with SARA in SNR and slightly superior in logSNR
across the range of measurement settings considered. As previously, we observe
that the confidence interval around the mean SNR and logSNR values for both
AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ is similar to that of uSARA and SARA,
reflecting the robustness of the proposed method. The reconstruction quality
of both CLEAN and the UNets trained end-to-end is much inferior to those of
AIRI and SARA approaches, for both SNR and logSNR and all measurement
conditions considered. Interestingly, notice that the UNets significantly
improve over CLEAN in linear scale, but perform significantly worse in
logarithmic scale, which confirms the difficulty of the UNets to recover very
faint emissions. This phenomenon was also noticed during training.

We recall that SARA and uSARA leverage the same prior model, encapsulated in
the averaged sparsity proximal operator. However, SARA relies on a constrained
data-fidelity term (see objective (6)), while uSARA uses an unconstrained
formulation enforcing less abruptly data fidelity, in the sense that it
imposes a soft penalty rather than a hard constraint (see objective (11)). A
by-product of our analysis is to confirm the superiority of the former over
the latter.

Figure 11 shows the reconstructed Hercules A images for the considered methods
and the $uv$-coverages shown in Figures 6(a), 6(b), and 6(c). The visual
analysis confirms the quantitative results from Figure 10. We acknowledge that
the AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ algorithms can create artefacts by
over-emphasising small structures (see the red zoom box in Figure 11(o)). SARA
exhibits very similar quality to AIRI-$\ell_{1}$, but not always with the same
resolution-to-artefacts tradeoff. The images suggest that the poor
reconstruction of CLEAN stems from a strong offset in intensity values (see in
particular the difference between the estimated higher intensity values in the
upper rectangular zoom from Figures 11(e),11(k), and 11(q), and that of the
groundtruth). Also note the important residual noise hiding low intensity
features of the image, as well as remaining artefacts for $\Delta
T=1\,\text{h}$ and $\Delta T=2\,\text{h}$. This figure also illustrates
clearly the strengths and shortcomings of the UNets trained in an end-to-end
manner. The relatively good SNR and low logSNR, together with the visual
results, confirm that the UNets manage to recover quite accurately the high
intensity emission but struggle to recover faint emission. Visually, the
results of the UNets present some hallucinated point-like sources (Muckley et
al., 2021), which we believe is the consequence of our synthetic training
database containing a large number of such sources.

Figure 10: Experiment 3 results: Reconstruction SNR (left) and logSNR (right)
as a function of the observation duration for AIRI-$\ell_{2}$,
AIRI-$\ell_{1}$, uSARA, SARA, CLEAN and the end-to-end UNets. Each point is an
average over the 20 simulated observations built from the four different
groundtruths from Figure 2 and the five different $uv$-patterns shown in
Figures 5 and 6(c). For AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$, $\sigma$ is
chosen at the heuristic. For uSARA, $\gamma\lambda$ is fine-tuned manually to
optimize the logSNR, resulting in values between $3$ and $6$ times below the
heuristic (in line with Figure 9). Error bars show the $95\%$ confidence
interval.

####  4.4.3 Computational cost and reconstruction times

The computational costs for the various approaches considered are reported in
Figure 12(b). Each point on the graphs is an average over the 20 simulated
observations built from the four different groundtruths from Figure 2 and the
five different $uv$-patterns shown in Figures 5 and 6(c). Error bars show the
95% confidence intervals.

Figure 12(a) displays the average time per iteration, independently for each
of the forward and backward steps of both uSARA and AIRI, as a function of the
observation duration. We recall once more that uSARA relies on the same
algorithmic structure as AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$. In both cases,
the underpinning FB algorithm alternates between a (forward) gradient-descent
step enforcing data fidelity and a (backward) regularization step enforcing a
prior model via the application of either the average sparsity proximal
operator (uSARA, Algorithm 1) or a learned denoiser (AIRI, Algorithm 3). We
note that, in such an algorithmic structure, the computational cost of the
regularization step scales linearly with the image size, while that of the
gradient-descent step scales linearly with the data size.

Firstly, we observe from the figure that the computation time of the gradient
step naturally scales linearly with the observation duration, which is indeed
proportional to data size. Secondly, we acknowledge that the computation of
the proximal operator within uSARA is slow, which is due to its sub-iterative
nature. For the image and data sizes considered in our simulations, and given
our CPU implementation, it always dominates the computation time of the
gradient, on average per iteration of Algorithm 1. This conclusion holds
despite the fact that we have carefully optimized the external parameters of
Algorithm 1 and Algorithm 2 to optimize the reconstruction time ($K$ in
Algorithm 1, and the tolerance criterion $\xi_{2}$ and dual variable
$\bm{v}_{0}$ in Algorithm 2). In that regard, we note that an efficient
initialization for the dual variable $\bm{v}_{0}$ consists in setting it as
the dual variable from the previous use of Algorithm 2. Using such an
initialization, the number of sub-iterations for Algorithm 2 progressively
decreases with the iteration count in Algorithm 1. Thirdly, using a standard
GPU implementation of the AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ denoisers, their
application is two orders of magnitude faster than the computation of the
proximal operator. A CPU implementation of the denoisers makes them even
slower than the proximal operator in uSARA. For the image and data sizes at
stake, the AIRI denoisers on GPU are also significantly faster than the
computation of the gradient, on average per iteration of Algorithm 1.

Figure 12(b) displays the total reconstruction times for uSARA and AIRI as a
function of the observation duration, and in comparison with SARA and CLEAN.
The AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ reconstruction times with denoiser on
GPU scale almost linearly with the observation duration, due to the domination
of the gradient step. For uSARA, and for AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$
with denoiser on CPU, the dependency on observation duration is less severe
given the domination of the regularization step at each iteration. Overall,
given the speed up in the regularization step, the total reconstruction time
of AIRI with denoiser on GPU is significantly ($5$ to $10$ times) lower than
the one of uSARA. AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ with denoiser on CPU are
much slower, in fact slightly slower than uSARA itself. SARA offers
reconstruction time between those of AIRI and uSARA. CLEAN remains
significantly faster than AIRI ($4$ to $15$ times), but the fastest method
among all is by far the UNet trained end-to-end. We underline that this is to
be expected since, not only UNets were run on GPU, but the reconstruction with
the UNet only requires a single inference step, as opposed to the other
methods that are iterative.

To summarize, the application of an AIRI denoiser (on GPU) is orders of
magnitude faster than that of the average sparsity proximal operator of uSARA
(on CPU). Substituting the former by the latter does not only provide superior
imaging quality, but significantly reduces the computational cost of the
regularization step. For image and data sizes where the uSARA regularization
step dominates over the gradient-descent data-fidelity step, a cost reduction
in the former, substituting the proximal operator for a DNN denoiser leads to
significant speed up of the algorithm.

(a) Groundtruth
(SNR, logSNR)

(b) Dirty image, $\Delta T=1\,\text{h}$
(SNR, logSNR)

(c) Dirty image, $\Delta T=2\,\text{h}$
(SNR, logSNR)

(d) Dirty image, $\Delta T=4\,\text{h}$
(SNR, logSNR)

(e) CLEAN, $\Delta T=1\,\text{h}$
$(4.90\,\text{dB},9.26\,\text{dB},4\,\text{min})$

(f) UNet, $\Delta T=1\,\text{h}$
$(11.05\,\text{dB},1.24\,\text{dB},7\,\text{ms})$

(g) SARA, $\Delta T=1\,\text{h}$
$(23.78\,\text{dB},22.99\,\text{dB},38\text{min})$

(h) uSARA, $\Delta T=1\,\text{h}$
$(19.73\,\text{dB},21.47\,\text{dB},3.5\text{h})$

(i) AIRI-$\ell_{2}$, $\Delta T=1\,\text{h}$
$(22.40\,\text{dB},21.22\,\text{dB},12\text{min})$

(j) AIRI-$\ell_{1}$, $\Delta T=1\,\text{h}$
$(24.13\,\text{dB},21.80\,\text{dB},12\,\text{min})$

(k) CLEAN, $\Delta T=2\,\text{h}$
$(5.05\,\text{dB},8.96\,\text{dB},3\,\text{min})$

(l) UNet, $\Delta T=2\,\text{h}$
$(12.22\,\text{dB},1.19\,\text{dB},7\,\text{ms})$

(m) SARA, $\Delta T=2\,\text{h}$
$(25.44\,\text{dB},23.63\,\text{dB},49\,\text{min})$

(n) uSARA, $\Delta T=2\,\text{h}$
$(20.53\,\text{dB},22.61\,\text{dB},2.9\,\text{h})$

(o) AIRI-$\ell_{2}$, $\Delta T=2\,\text{h}$
$(23.54\,\text{dB},23.15\,\text{dB},20\,\text{min})$

(p) AIRI-$\ell_{1}$, $\Delta T=2\,\text{h}$
$(25.30\,\text{dB},24.76\,\text{dB},20\,\text{min})$

(q) CLEAN, $\Delta T=4\,\text{h}$
$(5.19\,\text{dB},8.95\,\text{dB},4\,\text{min})$

(r) UNet, $\Delta T=4\,\text{h}$
$(13.81\,\text{dB},1.0\,\text{dB},7\,\text{ms})$

(s) SARA, $\Delta T=4\,\text{h}$
$(26.34\,\text{dB},24.69\,\text{dB},1.3\,\text{h})$

(t) uSARA, $\Delta T=4\,\text{h}$
$(22.42\,\text{dB},23.97\,\text{dB},3.4\,\text{h})$

(u) AIRI-$\ell_{2}$, $\Delta T=4\,\text{h}$
$(24.41\,\text{dB},23.77\,\text{dB},39\,\text{min})$

(v) AIRI-$\ell_{1}$, $\Delta T=4\,\text{h}$
$(26.41\,\text{dB},25.81\,\text{dB},39\,\text{min})$

Figure 11: Experiment 3 results: Influence of the observation duration $\Delta
T$ on the reconstruction quality for the different algorithms with the
pointing position from Figure 6; the associated $uv$-patterns are displayed in
Figures 6(a), 6(b), and 6(c). Second, third and last row show reconstruction
results for $\Delta T=1\,\text{h}$, $\Delta T=2\,\text{h}$ and $\Delta
T=4\,\text{h}$, respectively. For these three rows, each column from left to
right shows estimated model images obtained with CLEAN, UNets, SARA, uSARA,
AIRI-$\ell_{2}$, and AIRI-$\ell_{1}$, respectively. Below each reconstruction
we indicate the reconstruction metrics and time (SNR, logSNR, reconstruction
time). All images are shown in logarithmic scale.

(a) Average time per iteration.

(b) Total timing for each method.

Figure 12: Experiment 3 results: reconstruction times for AIRI-$\ell_{2}$,
AIRI-$\ell_{1}$, uSARA, SARA, CLEAN, and the end-to-end UNets (on GPU)

. AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ exhibit identical behaviours, with only
AIRI-$\ell_{2}$ results are reported and referred to as AIRI on the graphs.
Each point is an average over the 20 simulated observations. Error bars show
the (very small and virtually invisible) 95% confidence interval. For
AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$, all reconstruction times are reported for
both GPU (solid red lines) and CPU implementation (dashed red lines) of
$\operatorname{D}$. Left: average time per iteration in seconds (s.) for the
forward and backward operators involved in Algorithm 1 ($\nabla f$ and
$\operatorname{prox}_{\gamma\lambda g}$) and Algorithm 3 ($\nabla f$ and
denoiser $\operatorname{D}$), as a function of the observation duration.
Right: average total reconstruction times for AIRI-$\ell_{2}$,
AIRI-$\ell_{1}$, uSARA, SARA, CLEAN, and the end-to-end UNets, as a function
of the observation duration.

##  5 Conclusion & Future Work

###  5.1 Conclusion for both AIRI & uSARA

In this work, we have proposed a new algorithmic framework for scalable
precision RI imaging, dubbed AIRI. The approach consists in encapsulating a
prior image model in a DNN denoiser, and plugging it in lieu of the proximal
regularization operator of an optimization algorithm of the SARA family for
image reconstruction. AIRI inherits the robustness and interpretability of
optimization approaches and the learning power and speed of networks, while
avoiding generalizability limitations of end-to-end networks, for which
training involves the details of the measurement model.

More specifically, we have firstly designed a realistic low dynamic range
training database from optical images. We have also established a method to
train firmly nonexpansive denoisers at a noise level inferred from the target
dynamic ranges of reconstruction, with either an $\ell_{2}$ or $\ell_{1}$
loss, including a procedure for on-the-fly database dynamic range enhancement.
Finally, the AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ algorithms resulting from
plugging the denoisers into the FB algorithm were validated in simulation for
reconstruction of images involving complex structure with diffuse and faint
emission across the field of view (at a dynamic range of $10^{4}$). The
benchmark algorithms included CLEAN, the state-of-the-art SARA algorithm, as
well as its new unconstrained version uSARA, which only differs from
AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ by the use of the average sparsity
proximal regularization operator instead of a learned denoiser in the FB
algorithm.

Our results show that the first incarnations of AIRI are already competitive
with, if not superior to, the advanced optimization algorithms of the SARA
family in imaging quality, well beyond the capabilities of CLEAN and of the
UNets trained in an end-to-end fashion. In particular, while the AIRI
denoisers are trained on the same database as the UNets, the PnP solutions do
not exhibit hallucinated artefacts, which illustrates the robustness of AIRI
compared to the end-to-end UNets. With regards to computation time, the AIRI
approach was shown to provide a significant acceleration potential over uSARA
and SARA, thanks to the inference speed of DNN denoisers over the averaged
sparsity proximal operator, but remain significantly slower than CLEAN, itself
orders of magnitude slower than the UNets. The results also confirm the
validity of our heuristic to set the training noise level $\sigma$, consisting
in enforcing equality of the training and target dynamic range of the
reconstruction, and the validity of the database exponentiation approach (with
parameter $a$), consisting in enforcing equality of the final nominal dynamic
range of the database and target dynamic range of reconstruction. This
provides a significant advantage over uSARA, for which the regularization
parameter $\lambda$ requires manual fine-tuning, even though the fact that the
corresponding heuristic (and its $1/3$-corrected version) provides an
appropriate reference upper-bound for the optimal value is an interesting
result.

###  5.2 Current limitations & future work

In what follows, we discuss current limitations of the AIRI framework, and
future work. This also includes considerations for further developments of
uSARA and end-to-end DNNs.

Firstly, the modified DnCNN denoiser architecture used, the $\ell_{2}$ and
$\ell_{1}$ losses considered, as well as the set of optical images used to
build a training database constitute the most basic choices leading to our
first AIRI-$\ell_{2}$ and AIRI-$\ell_{1}$ incarnations of the framework. We
anticipate that future work, building a richer database from RI observation,
and considering more advanced losses, such as adversarial losses (Wang et al.,
2018), and more advanced network architectures, such as UNets (Hurault et al.,
2022), can provide a quantum jump in imaging precision with AIRI. We also
acknowledge that our approach to enforce the firm nonexpansiveness of the
denoiser is only approximate, and relies on fine-tuning the regularization
parameter $\kappa$ in the training loss. Simpler and more precise approaches
to ensure firm nonexpansiveness should be investigated in the future (_e.g._
leveraging 1-Lipschitz layers). We note however that, in the context of our
experiments, the proposed approach has, in practice, provided a robust way to
ensure convergence of the resulting PnP algorithms.

Secondly, fully optimized and parallelized implementations of AIRI algorithms,
but also SARA, or uSARA, should be investigated, leveraging GPUs not only for
DNN denoisers, but also for the implementation of the linear (measurement or
regularization) operators involved. This should lead to further acceleration
of these approaches, and possibly participate to closing the gap with CLEAN in
terms of reconstruction time. We also note that accelerated versions of uSARA
are achievable, in particular leveraging preconditioning strategies (Repetti &
Wiaux, 2021).

Thirdly, with regards to the details of the target imaging modality, we have
investigated monochromatic intensity imaging on small fields of view only,
with mild, if any, super-resolution factor with respect to the nominal
resolution of the observation, as set by the largest baseline. Extensions of
the framework to further super-resolution functionality and to wideband
polarization imaging should be contemplated. Also, a by-product of the present
analysis is to confirm the superiority of SARA over uSARA, due to the
constrained approach to data fidelity in SARA, versus the unconstrained
approach underpinning uSARA. In this context, the superior results of AIRI
over uSARA are entirely related to the capability to encapsulate a better
prior model than the average sparsity model in learned DNN denoisers. On the
one hand, a first possible enhancement of the AIRI approach would consist in
developing a PnP version of SARA (as opposed to uSARA), simultaneously taking
advantage of networks to learn prior models and of the constrained data-
fidelity formulation. On the other hand, to date, the unconstrained
formulation underpinning uSARA, AIRI-$\ell_{2}$, and AIRI-$\ell_{1}$ is a
critical building block of the only optimization-based algorithm proposed for
RI imaging capable of handling jointly DDE calibration and imaging (Repetti et
al., 2017; Dabbech et al., 2021). This is a strong justification for the study
of unconstrained versus constrained formulations of the imaging module.

Fourthly, the proposed UNet end-to-end implementation showed interesting
results, albeit significantly outperformed by the proposed AIRI algorithms. A
more thorough study should be undertaken, leveraging more advanced
architectures, such as unfolded networks that would enable to incorporate
knowledge about the measurement operator within the architecture.

Finally, the present work is simulation-based only. Dabbech et al. (2022),
Wilber et al. (2022a), and Wilber et al. (2022b) validate AIRI and uSARA on
real wide-field high-resolution high-dynamic range data from MeerKAT and
ASKAP, relying on a highly parallelised implementation of the measurement
operator tailored to the wide-field setting.

## Acknowledgements

The research of M. T. was funded by Heriot-Watt under the James Watt
Scholarship scheme. The work of M. T., A. D., and Y. W. was supported by EPSRC
under grants EP/T028270/1 and ST/W000970/1. The computing resources came from
the Cirrus UK National Tier-2 HPC Service at EPCC (http://www.cirrus.ac.uk)
funded by the University of Edinburgh and EPSRC (EP/P020267/1), partly through
time allocation under the SUSA project, and partly through GPU resources
directly provided by EPCC (Adrian Jackson). Credits for the 32 images used in
our training database go to NOIRLab/NSF/AURA/H.Schweiker/WIYN/T.A.Rector
(University of Alaska Anchorage). The radio images where taken from: for
Hercules A: R. Perley and W. Cotton (NRAO/AUI/NSF), for 3c353: NRAO/VLA, for
Centaurus A: NRAO/AUI/NSF/Univ.Hertfordshire/M.Hardcastle.

## Data Availability

AIRI code will be made available as part of a later release of and joint
parallel toolbox for AIRI and uSARA.

## References

  * Abdulaziz et al. (2016) Abdulaziz A., Dabbech A., Onose A., Wiaux Y., 2016, in European Signal Processing Conf.. pp 388–392
  * Abdulaziz et al. (2019) Abdulaziz A., Dabbech A., Wiaux Y., 2019, MNRAS, 489, 1230
  * Ahmad et al. (2020) Ahmad R., Bouman C. A., Buzzard G. T., Chan S., Liu S., Reehorst E. T., Schniter P., 2020, IEEE Signal Processing Magazine, 37, 105
  * Akiyama et al. (2019) Akiyama K., et al., 2019, The Astrophysical Journal Letters, 875, L4
  * Akiyama et al. (2022) Akiyama K., et al., 2022, The Astrophysical Journal Letters, 930, L14
  * Arras et al. (2019) Arras P., Frank P., Leike R., Westermann R., Enßlin T. A., 2019, A&A, 627, A134
  * Attouch et al. (2013) Attouch H., Bolte J., Svaiter B. F., 2013, Math. Programming, 137, 91
  * Balestriero & Baraniuk (2021) Balestriero R., Baraniuk R., 2021, arXiv preprint arXiv:2104.00219
  * Baraniuk (2007) Baraniuk R. G., 2007, IEEE Signal Processing Magazine, 24, 118
  * Bauschke & Combettes (2017) Bauschke H. H., Combettes P. L., 2017, Convex analysis and monotone operator theory in Hilbert spaces. Springer
  * Beck & Teboulle (2009) Beck A., Teboulle M., 2009, SIAM J. on Imaging Sci., 2, 183
  * Bhatnagar & Cornwell (2004) Bhatnagar S., Cornwell T., 2004, A&A, 426, 747
  * Birdi et al. (2018) Birdi J., Repetti A., Wiaux Y., 2018, MNRAS, 478, 4442
  * Birdi et al. (2020) Birdi J., Repetti A., Wiaux Y., 2020, MNRAS, 492, 3509
  * Bonaldi & Braun (2018) Bonaldi A., Braun R., 2018, arXiv preprint arXiv:1811.10454
  * Boyd et al. (2011) Boyd S., Parikh N., Chu E., 2011, Distributed optimization and statistical learning via the alternating direction method of multipliers. Now Publishers Inc
  * Bredies & Holler (2020) Bredies K., Holler M., 2020, Inverse Problems, 36, 123001
  * Bredies et al. (2010) Bredies K., Kunisch K., Pock T., 2010, SIAM J. on Imaging Sci., 3, 492
  * Briggs (1995) Briggs D. S., 1995, in American Astron. Soc. Meeting Abstracts. p. 112.02
  * Bryan et al. (2014) Bryan G. L., et al., 2014, ApJS, 211, 19
  * Buades et al. (2011) Buades A., Coll B., Morel J.-M., 2011, Image Processing On Line, 1, 208
  * Cai et al. (2018) Cai X., Pereyra M., McEwen J. D., 2018, MNRAS, 480, 4154
  * Candes & Demanet (2003) Candes E., Demanet L., 2003, Comptes Rendus Mathematique, 336, 395
  * Candès et al. (2006) Candès E. J., Romberg J., Tao T., 2006, IEEE Trans. on Inf. Theor., 52, 489
  * Candes et al. (2008) Candes E. J., Wakin M. B., Boyd S. P., 2008, J. of Fourier analysis and applications, 14, 877
  * Carrillo et al. (2012) Carrillo R. E., McEwen J. D., Wiaux Y., 2012, MNRAS, 426, 1223
  * Carrillo et al. (2014) Carrillo R. E., McEwen J. D., Wiaux Y., 2014, MNRAS, 439, 3591
  * Chan et al. (2016) Chan S. H., Wang X., Elgendy O. A., 2016, IEEE Trans. on Comput. Imaging, 3, 84
  * Chouzenoux et al. (2014) Chouzenoux E., Pesquet J.-C., Repetti A., 2014, J. of Optimization Theor. and Applications, 162, 107
  * Cohen et al. (2021) Cohen R., Elad M., Milanfar P., 2021, SIAM J. on Imaging Sci., 14, 1374
  * Combettes & Pesquet (2011) Combettes P. L., Pesquet J.-C., 2011, in , Fixed-point algorithms for inverse problems in science and engineering. Springer, pp 185–212
  * Combettes & Pesquet (2015) Combettes P. L., Pesquet J.-C., 2015, SIAM Journal on Optimization, 25, 1221
  * Connor et al. (2022) Connor L., Bouman K. L., Ravi V., Hallinan G., 2022, MNRAS, 514, 2614
  * Cornwell et al. (2008) Cornwell T. J., Golap K., Bhatnagar S., 2008, IEEE J. of Selected Topics in Signal Processing, 2, 647
  * Dabbech et al. (2015) Dabbech A., Ferrari C., Mary D., Slezak E., Smirnov O., Kenyon J. S., 2015, A&A, 576, A7
  * Dabbech et al. (2017) Dabbech A., Wolz L., Pratley L., McEwen J. D., Wiaux Y., 2017, MNRAS, 471, 4300
  * Dabbech et al. (2018) Dabbech A., Onose A., Abdulaziz A., Perley R. A., Smirnov O. M., Wiaux Y., 2018, MNRAS, 476, 2853
  * Dabbech et al. (2021) Dabbech A., Repetti A., Perley R. A., Smirnov O. M., Wiaux Y., 2021, MNRAS, 506, 4855
  * Dabbech et al. (2022) Dabbech A., Terris M., Jackson A., Ramatsoku M., Smirnov O., Wiaux Y., 2022, arXiv preprint arXiv:2207.11336
  * Dabov et al. (2007) Dabov K., Foi A., Katkovnik V., Egiazarian K., 2007, IEEE Trans. on Image Processing, 16, 2080
  * Do & Vetterli (2003) Do M. N., Vetterli M., 2003, in , Vol. 10, Studies in Comp. Math.. Elsevier, pp 83–105
  * Donoho (2006) Donoho D. L., 2006, IEEE Trans. on Inf. Theor., 52, 1289
  * Donoho & Johnstone (1995) Donoho D. L., Johnstone I. M., 1995, J. of the American Statistical Association, 90, 1200
  * Eckstein & Bertsekas (1992) Eckstein J., Bertsekas D. P., 1992, Math. Programming, 55, 293
  * Fessler (2020) Fessler J. A., 2020, IEEE Signal Processing Magazine, 37, 33
  * Garsden et al. (2015) Garsden H., et al., 2015, A&A, 575, A90
  * Gheller & Vazza (2021) Gheller C., Vazza F., 2021, MNRAS
  * Golub & Van Loan (2013) Golub G. H., Van Loan C. F., 2013, Matrix computations. JHU press
  * Goodfellow et al. (2015) Goodfellow I., Shlens J., Szegedy C., 2015, in Int. Conf. on Learning Representations.
  * Gupta et al. (2018) Gupta H., Jin K. H., Nguyen H. Q., McCann M. T., Unser M., 2018, IEEE Trans. on Medical Imaging, 37, 1440
  * Hertrich et al. (2021) Hertrich J., Neumayer S., Steidl G., 2021, Linear Algebra and its Applications, 631, 203
  * Högbom (1974) Högbom J., 1974, A&AS, 15, 417
  * Hurault et al. (2022) Hurault S., Leclaire A., Papadakis N., 2022, arXiv preprint arXiv:2201.13256
  * Jin et al. (2017) Jin K. H., McCann M. T., Froustey E., Unser M., 2017, IEEE Trans. on Image Processing, 26, 4509
  * Kingma & Ba (2015) Kingma D. P., Ba J., 2015, in Int. Conf. on Learning Representations.
  * Levin et al. (2009) Levin A., Weiss Y., Durand F., Freeman W. T., 2009, in 2009 IEEE Conf. on Comput. Vision and Pattern Recognition. pp 1964–1971
  * Li et al. (2011) Li F., Cornwell T. J., de Hoog F., 2011, A&A, 528, A31
  * Luisier et al. (2007) Luisier F., Blu T., Unser M., 2007, IEEE Trans. on Image Processing, 16, 593
  * Maddox et al. (2021) Maddox W., Tang S., Moreno P., Wilson A. G., Damianou A., 2021, in International Conference on Artificial Intelligence and Statistics. pp 2737–2745
  * Mairal et al. (2009) Mairal J., Bach F., Ponce J., Sapiro G., 2009, in Int. Conf. on Machine Learning. pp 689–696
  * Mallat (1999) Mallat S., 1999, A wavelet tour of signal processing. Elsevier
  * Muckley et al. (2021) Muckley M. J., et al., 2021, IEEE transactions on medical imaging, 40, 2306
  * Nguyen et al. (2015) Nguyen A., Yosinski J., Clune J., 2015, in IEEE Conf. on Comput. Vision and Pattern Recognition. pp 427–436
  * Offringa & Smirnov (2017) Offringa A. R., Smirnov O., 2017, MNRAS, 471, 301
  * Onose et al. (2016) Onose A., Carrillo R. E., Repetti A., McEwen J. D., Thiran J.-P., Pesquet J.-C., Wiaux Y., 2016, MNRAS, 462, 4314
  * Onose et al. (2017) Onose A., Dabbech A., Wiaux Y., 2017, MNRAS, 469, 938
  * Pang et al. (2018) Pang T., Du C., Dong Y., Zhu J., 2018, Advances in Neural Inf. Processing Syst., 31
  * Paszke et al. (2017) Paszke A., et al., 2017, in Advances in Neural Inf. Processing Syst. Workshop.
  * Pesquet & Repetti (2015) Pesquet J.-C., Repetti A., 2015, J. Nonlinear Convex Analysis, 16, 2353
  * Pesquet et al. (2021) Pesquet J.-C., Repetti A., Terris M., Wiaux Y., 2021, SIAM J. on Imaging Sci., 14, 1206
  * Pratley et al. (2018) Pratley L., McEwen J. D., d’Avezac M., Carrillo R. E., Onose A., Wiaux Y., 2018, MNRAS, 473, 1038
  * Pruša (2012) Pruša Z., 2012, PhD thesis, Brno university of technology
  * Repetti & Wiaux (2020) Repetti A., Wiaux Y., 2020, in IEEE Int. Conf. on Acoustics, Speech and Signal Processing. pp 1434–1438
  * Repetti & Wiaux (2021) Repetti A., Wiaux Y., 2021, SIAM J. on Optimization, 31, 1215
  * Repetti et al. (2017) Repetti A., Birdi J., Dabbech A., Wiaux Y., 2017, MNRAS, 470, 3981
  * Repetti et al. (2018) Repetti A., Pereyra M., Wiaux Y., 2018, in European Signal Processing Conf.. pp 2668–2672
  * Repetti et al. (2019) Repetti A., Pereyra M., Wiaux Y., 2019, SIAM J. on Imaging Sci., 12, 87
  * Romano et al. (2017) Romano Y., Elad M., Milanfar P., 2017, SIAM J. on Imaging Sci., 10, 1804
  * Rudin et al. (1992) Rudin L. I., Osher S., Fatemi E., 1992, Phys. D: nonlinear phenomena, 60, 259
  * Ryu et al. (2019) Ryu E., Liu J., Wang S., Chen X., Wang Z., Yin W., 2019, in Int. Conf. on Machine Learning. pp 5546–5557
  * Schwab & Cotton (1983) Schwab F. R., Cotton W. D., 1983, AJ, 88, 688
  * Schwarz (1978) Schwarz U., 1978, A&A, 65, 345
  * Selesnick et al. (2005) Selesnick I. W., Baraniuk R. G., Kingsbury N. C., 2005, IEEE Signal Processing Magazine, 22, 123
  * Terris et al. (2019) Terris M., Abdulaziz A., Dabbech A., Jiang M., Repetti A., Pesquet J.-C., Wiaux Y., 2019, in SPARS workshop.
  * Terris et al. (2020) Terris M., Repetti A., Pesquet J.-C., Wiaux Y., 2020, in IEEE Int. Conf. on Acoustics, Speech and Signal Processing. pp 8658–8662
  * Thiébaut & Young (2017) Thiébaut É., Young J., 2017, JOSA A, 34, 904
  * Thompson et al. (2017) Thompson A. R., Moran J. M., Swenson Jr G. W., 2017, Interferometry and synthesis in radio astronomy. Springer Nat.
  * Thouvenin et al. (2022a) Thouvenin P.-A., Abdulaziz A., Dabbech A., Repetti A., Wiaux Y., 2022a, accepted for publication in MNRAS, HWU research portal ID: 25219759
  * Thouvenin et al. (2022b) Thouvenin P.-A., Dabbech A., Jiang M., Abdulaziz A., Thiran J.-P., Jackson A., Wiaux Y., 2022b, submitted, HWU research portal ID: 52577083
  * Venkatakrishnan et al. (2013) Venkatakrishnan S. V., Bouman C. A., Wohlberg B., 2013, in 2013 IEEE Global Conf. on Signal and Inf. Processing. pp 945–948
  * Vidal et al. (2020) Vidal A. F., De Bortoli V., Pereyra M., Durmus A., 2020, SIAM J. on Imaging Sci., 13, 1945
  * Wang et al. (2015) Wang L., Xiong Z., Gao D., Shi G., Wu F., 2015, Applied optics, 54, 848
  * Wang et al. (2018) Wang X., Yu K., Wu S., Gu J., Liu Y., Dong C., Qiao Y., Change Loy C., 2018, in European Conf. on Comput. Vision Workshops. pp 0–0
  * Wiaux et al. (2009) Wiaux Y., Jacques L., Puy G., Scaife A. M., Vandergheynst P., 2009, MNRAS, 395, 1733
  * Wilber et al. (2022a) Wilber A., Dabbech A., Jackson A., Wiaux Y., 2022a, in preparation, HWU research portal ID: 52551784
  * Wilber et al. (2022b) Wilber A., Dabbech A., Terris M., Jackson A., Wiaux Y., 2022b, in preparation, HWU research portal ID: 54234902
  * Xie et al. (2019) Xie Q., Zhou M., Zhao Q., Meng D., Zuo W., Xu Z., 2019, in IEEE/CVF Conf. on Comput. Vision and Pattern Recognition. pp 1585–1594
  * Yang et al. (2010) Yang J., Wright J., Huang T. S., Ma Y., 2010, IEEE Trans. on Image Processing, 19, 2861
  * Zbontar et al. (2018) Zbontar J., et al., 2018, arXiv preprint arXiv:1811.08839
  * Zhang et al. (2017) Zhang K., Zuo W., Chen Y., Meng D., Zhang L., 2017, IEEE Trans. on Image Processing, 26, 3142
  * Zhang et al. (2019) Zhang K., Zuo W., Zhang L., 2019, in IEEE/CVF Conf. on Comput. Vision and Pattern Recognition. pp 1671–1681
  * Zhang et al. (2021) Zhang K., Li Y., Zuo W., Zhang L., Van Gool L., Timofte R., 2021, IEEE Trans. on Pattern Analysis and Machine Intelligence
  * Zhou et al. (2020) Zhou Y., Jiao J., Huang H., Wang Y., Wang J., Shi H., Huang T., 2020, in AAAI Conf. on Artificial Intelligence. pp 13074–13081

##  Appendix A On the chosen RI imaging model

We discuss two points regarding the generality of the measurement model in
(1). Firstly, at the dynamic range of interest for current and future
observations, direction-dependent effects (DDEs) of atmospheric and
instrumental origin complicate the RI measurement equation. More precisely, at
each time instant, the visibility associated with a pair of antennas is formed
from the modulation of the radio sky with the product of a DDE pattern
specific to each antenna. The DDEs are typically unknown and need to be
calibrated, leading to a blind deconvolution problem, which requires the
design of joint calibration and imaging approaches, typically alternating
between imaging and calibration modules (Repetti et al., 2017; Dabbech et al.,
2021). One exception to this, is the so-called $w$-term induced by the
projection of the baseline associated with each visibility on the line of
sight, which becomes non-negligible on wide fields of view, and thus acts as a
known DDE. In practice, DDEs can be integrated into extended convolution
kernels in $\bm{\mathrm{G}}$ (Dabbech et al., 2017). In this work, we assume a
measurement operator $\bm{\Phi}$ not affected by any known or unknown DDE,
considering a pure imaging problem on small fields of view.

Secondly, the measurement model (1) assumes noise with constant variance
across visibilities. In full generality, the noise is not white, and a
diagonal whitening matrix $\bm{\Theta}\in\mathbb{R}^{m\times m}$ with entries
equal to the inverse noise standard deviation per visibility can be applied to
the original measurement vector, ensuring that the resulting visibility vector
is indeed affected by white noise. In this case problem (1) still holds, with
a measurement operator $\bm{\Phi}=\bm{\Theta\mathrm{GFZ}}$, and $\bm{y}$ and
$\bm{e}$ the visibility and noise vectors after application of the whitening
matrix. This operation is known as natural weighting. This scheme is used by
the algorithms of the SARA family to ensure the negative log-likelihood
interpretation of its data-fidelity terms (Carrillo et al., 2012). More
complex weighting schemes, such as uniform and Briggs weighting, are used by
CLEAN (Briggs, 1995), which involve additional multiplicative terms in
$\bm{\Theta}$. Both approaches add to natural weighting by downweighting
visibilities in regions of the $uv$-plane with high sampling density, with the
aim to reduce the sidelobes of the dirty beam and improve the resolution of
the CLEAN reconstruction, though at the expense of sensitivity. Interestingly,
Onose et al. (2017) showed, for algorithms of the SARA family specifically,
that, while keeping to natural weighting, using sampling density information
in acceleration strategies such as preconditioning are very efficient at
bringing jointly optimal resolution and dynamic range. In other words, the
resolution-versus-sensitivity trade-off disappears asymptotically at
convergence. In this work, without loss of generality, our simulations
consider i.i.d. Gaussian random measurement noise so that
$\bm{\Theta}=\mathbf{I}$ for natural weighting, where $\mathbf{I}$ denotes the
identity matrix. All algorithms discussed (SARA, uSARA, AIRI) use natural
weighting, while CLEAN is implemented with uniform weighting (see Section
4.1.2).

##  Appendix B On the PnP convergence

Several works have recently focused on restoring the convergence of PnP
algorithms using DNN denoisers. The majority of these works assume a
nonexpansiveness constraint on the denoiser (Romano et al., 2017; Ryu et al.,
2019; Terris et al., 2020; Cohen et al., 2021; Hertrich et al., 2021), but
this constraint often comes with restrictive assumptions either on the
algorithm, on the DNN architecture or on the operator $\nabla f$.

In Pesquet et al. (2021), casting the problem in the framework of monotone
operator theory, we have shown that the convergence of a PnP algorithm can be
ensured by introducing a well-chosen regularization in the training loss of
the denoiser $\operatorname{D}$ enforcing its firm nonexpansiveness, without
any form of limitation on the DNN architecture itself. Importantly, the
approach also offers a generalized notion of characterization of its limit as
the solution to monotone inclusion problems, more general than pure
optimization problems, hence opening a path towards generalizing the usual
Bayesian MAP interpretation of optimization solutions. We summarize those
results in the remainder of this section.

Once again concentrating on a differentiable $f$ such as (10), minimization
problems like (3) can be reformulated (through their optimality conditions) as

$\text{find $\widehat{\bm{x}}\in\mathbb{R}^{n}$ s.t. }\qquad 0\in\nabla
f(\widehat{\bm{x}})+\lambda\partial r(\widehat{\bm{x}}),$ (26)

where $\partial r$ denotes the subdifferential of $r$, the subdifferential of
$r\in\Gamma_{0}(\mathbb{R}^{n})$ being the set-valued operator defined as
$\partial
r(\bm{x})=\\{\bm{u}\in\mathbb{R}^{n}|(\forall\bm{z}\in\mathbb{R}^{n})\langle\bm{z}-\bm{x},\bm{u}\rangle+r(\bm{x})\leq
r(\bm{z})\\}$. When $r$ is also differentiable, the subdifferential coincides
with the gradient (_i.e._ $\partial r=\nabla r$). Yet, the subdifferential is
also defined for non-differentiable functions such as
$\iota_{\mathcal{B}(\bm{y},\epsilon)}$; we refer the reader to Bauschke &
Combettes (2017) for more information. Problems as (26) are part of the class
of so-called monotone inclusion problems and encompass unconstrained
minimization problems such as uSARA (11). They can however take a more general
flavour and write as

$\text{find $\widehat{\bm{x}}\in\mathbb{R}^{n}$ s.t. }\qquad 0\in\nabla
f(\widehat{\bm{x}})+\operatorname{A}(\widehat{\bm{x}}),$ (27)

where $\operatorname{A}$ is an operator satisfying some maximal monotonicity
assumptions (Bauschke & Combettes, 2017). It is important to stress that this
class of operators includes, but is not restricted to, subdifferentials.
Problem (27) is thus in general not equivalent to a minimization problem.

Solving such inclusion problems can be done using a generalized version of the
proximal FB algorithm (5) with

$(\forall
k\in\mathbb{N}),\qquad\bm{x}_{k+1}=\operatorname{J}_{\gamma\\!\operatorname{A}}(\bm{x}_{k}-\gamma\nabla
f(\bm{x}_{k})),$ (28)

where $\operatorname{J}_{\gamma\\!\operatorname{A}}$ is the so-called
resolvent of $\gamma\\!\operatorname{A}$, yielding the convergence of the
sequence $(\bm{x}_{k})_{k\in\mathbb{N}}$ to $\widehat{\bm{x}}$ satisfying
(27). As previously, the stepsize $\gamma$ must satisfy $0<\gamma<2/L$ where
$L$ is the Lipschitz constant of $\nabla f$.

The resolvent operator can be seen as a generalized version of the proximal
operator (4) (Bauschke & Combettes, 2017). In fact, (5) is a special case of
(28), since when $\operatorname{A}=\lambda\partial r$ for
$r\in\Gamma_{0}(\mathbb{R}^{n})$, one has
$\operatorname{J}_{\gamma\lambda\partial r}=\operatorname{prox}_{\gamma\lambda
r}$. A key characterization of resolvent operators is the following: the
operator $\operatorname{J}$ is the resolvent of a maximally monotone operator
$\operatorname{A}$ if and only if it is firmly nonexpansive. Technically, by
definition, an operator $\operatorname{J}$ is firmly nonexpansive if there
exists an operator $\operatorname{Q}$ with Lipschitz constant smaller or equal
to 1 such that
$\operatorname{J}=\left(\operatorname{I}+\operatorname{Q}\right)/2$, where
$\operatorname{I}$ is the identity operator. Given the identical structure of
PnP-FB and (28), it follows that $\operatorname{D}$ in (15) should be firmly
nonexpansive for the convergence of the PnP algorithm to be ensured, with the
PnP solutions characterized as solutions of a monotone inclusion problem (27).

Thus, in order to impose the firm nonexpansiveness of $\operatorname{D}$, it
is sufficient to impose that the Lipschitz constant of
$\operatorname{Q}=2\operatorname{D}-\operatorname{I}$ is less than 1, or
equivalently, that its Jacobian satisfies
$\|\bm{\nabla}\operatorname{Q}(\bm{x})\|_{\rm{S}}\leq 1$ for all $\bm{x}$.
Since most deep learning frameworks rely on Jacobian-vector products (Paszke
et al., 2017; Balestriero & Baraniuk, 2021; Maddox et al., 2021), one can,
given some $\bm{x}\in\mathbb{R}^{n}$, compute
$\|\bm{\nabla}\operatorname{Q}(\bm{x})\|_{\rm{S}}$ via the power method. Thus,
$\operatorname{D}$ would ideally be trained as a denoiser with a standard
training loss, but under the constraint that
$\|\bm{\nabla}\operatorname{Q}(\bm{x})\|_{\rm{S}}\leq 1$ for all $\bm{x}$.
Such a constraint is not trivial to impose in practice and, as proposed in
Pesquet et al. (2021), we relax it and introduce a regularization term in the
training loss that penalizes non firmly nonexpansive networks.

Given a database of noisy / groundtruth images $(\bm{z}_{s},\bm{u}_{s})_{1\leq
s\leq S}$, and denoting $\bm{\theta}\in\mathbb{R}^{c}$ the learnable
parameters of $\operatorname{D}$, the training loss ensuring the firm
nonexpansiveness of $\operatorname{D}$ reads:

$\displaystyle\underset{\bm{\theta}\in\mathbb{R}^{c}}{\text{minimize}}\,\frac{1}{S}\sum_{s=1}^{S}\bigg{(}\mathcal{L}(\operatorname{D}_{\bm{\theta}}(\bm{z}_{s})-\bm{u}_{s})+\kappa\operatorname{max}\\{\|\bm{\nabla}\operatorname{Q}_{\bm{\theta}}(\bm{z}_{s})\|_{\rm{S}},1-\varepsilon\\}\bigg{)},$
(29)

where $\mathcal{L}$ is a loss function and $\varepsilon>0$ is a safety margin
parameter. As underlined in Pesquet et al. (2021), this penalty does not
explicitly guarantee that the solution of this problem satisfies the
constraint $\|\bm{\nabla}\operatorname{Q}(\bm{z})\|_{\rm{S}}\leq 1$ on the
database, where it is effectively evaluated, let alone for all
$\bm{z}\in\mathbb{R}^{n}$. Various countermeasures are however deployed, which
have been demonstrated to be effective at ensuring PnP convergence in practice
(see Pesquet et al. (2021)). Firstly, the safety margin parameter
$\varepsilon>0$ is a mean to strengthen the penalty. Secondly, the term
$\|\bm{\nabla}\operatorname{Q}\|_{\rm{S}}$ is not computed exactly at
$\bm{z}_{s}$, but instead at a point sampled uniformly at random on the
segment $[\bm{z}_{s},\bm{u}_{s}]$, with the aim to explore a richer set of
images than those of the database only. Last but not least, the value of
$\kappa$ can, and should, be fine-tuned manually to optimize the balance
between the regularization term and the standard loss $\mathcal{L}$.

We conclude this section by underlining that, contrary to other works (Ryu et
al., 2019; Terris et al., 2020), the proposed approach does not require to
impose layer-wise Lipschitz constraints nor architectural restrictions on the
DNN.