An accelerated splitting algorithm for radio-interferometric imaging: when natural and uniform weighting meet

TL;DR

This paper introduces an accelerated primal-dual algorithm for radio-interferometric imaging that efficiently combines natural and uniform weighting, improving convergence speed, sensitivity, and resolution in large-scale data processing.

Contribution

It proposes a novel preconditioning strategy that incorporates sampling density and noise statistics into a primal-dual algorithm, unifying natural and uniform weighting schemes.

Findings

Accelerated convergence with realistic sampling patterns.

Enhanced sensitivity and resolution over traditional methods.

Effective reconstruction of real radio galaxy data.

Abstract

Next generation radio-interferometers, like the Square Kilometre Array, will acquire tremendous amounts of data with the goal of improving the size and sensitivity of the reconstructed images by orders of magnitude. The efficient processing of large-scale data sets is of great importance. We propose an acceleration strategy for a recently proposed primal-dual distributed algorithm. A preconditioning approach can incorporate into the algorithmic structure both the sampling density of the measured visibilities and the noise statistics. Using the sampling density information greatly accelerates the convergence speed, especially for highly non-uniform sampling patterns, while relying on the correct noise statistics optimises the sensitivity of the reconstruction. In connection to CLEAN, our approach can be seen as including in the same algorithmic structure both natural and uniform weighting, thereby simultaneously optimising both the resolution and the sensitivity. The method relies on a new non-Euclidean proximity operator for the data fidelity term, that generalises the projection onto the $\ell_2$ ball where the noise lives for naturally weighted data, to the projection onto a generalised ellipsoid incorporating sampling density information through uniform weighting. Importantly, this non-Euclidean modification is only an acceleration strategy to solve the convex imaging problem with data fidelity dictated only by noise statistics. We showcase through simulations with realistic sampling patterns the acceleration obtained using the preconditioning. We also investigate the algorithm performance for the reconstruction of the 3C129 radio galaxy from real visibilities and compare with multi-scale CLEAN, showing better sensitivity and resolution. Our MATLAB code is available online on GitHub.

Equations

$$y(\mathbf{u})=\int D(\mathbf{l},\mathbf{u})x(\mathbf{l})e^{-2i\pi \mathbf{u}\cdot \mathbf{l}}{\mathrm{d}}^2\mathbf{l},$$

$$\mathbf{y}=\mathsf{\Phi }\mathbf{x}+\mathbf{n},$$

$$\boldsymbol{x}^{(t)}=\boldsymbol{x}^{(t-1)}+\boldsymbol{\mathcal{T}}\Big{(}\boldsymbol{\mathsf{\Phi}}^{\dagger}\big{(}\boldsymbol{y}-\boldsymbol{\mathsf{\Phi}}\boldsymbol{x}^{(t-1)}\big{)}\Big{)}.$$

$$\underset{\mathbf{x}}{\min }\Vert {\mathsf{\Psi }}^{†}\mathbf{x}{\Vert }_0\mathrm{subject}\mathrm{to}\Vert \mathbf{y}-\mathsf{\Phi }\mathbf{x}{\Vert }_2\le ϵ\mathrm{and}\mathbf{x}\in {\mathbb{R}}_{+}^{\mathrm{N}}.$$

$$(\forall \mathbf{z}){\iota }_{\mathcal{C}}(\mathbf{z})\stackrel{\Delta }{=}\{\begin{array}{rl}0& \mathbf{z}\in \mathcal{C}\\ +\infty & \mathbf{z}\notin \mathcal{C}.\end{array}$$

$$\underset{\mathbf{x}}{\min }f(\mathbf{x})+l({\mathsf{W}}^{†}{\mathsf{\Psi }}^{†}\mathbf{x})+h(\mathsf{\Phi }\mathbf{x}).$$

$${\mathrm{prox}}_g(\mathbf{z})\stackrel{\Delta }{=}\underset{\bar{\mathbf{z}}}{\mathrm{argmin}}g(\bar{\mathbf{z}})+\frac{1}{2}\Vert \mathbf{z}-\bar{\mathbf{z}}{\Vert }_2^2.$$

$$\underset{\mathbf{x}}{\min }f(\mathbf{x})+\gamma \sum _{i=1}^b{l}_i({\mathsf{W}}_i^{†}{\mathsf{\Psi }}_i^{†}\mathbf{x})+\sum _{j=1}^d{h}_j({\mathsf{\Phi }}_j\mathbf{x}),$$

$$\displaystyle\min_{\begin{subarray}{c}\boldsymbol{u}_{i}\\ \boldsymbol{v}_{j}\end{subarray}}f^{*}\Bigg{(}\!\!-\sum_{i=1}^{b}\boldsymbol{\mathsf{\Psi}}_{i}\boldsymbol{\mathsf{W}}_{i}\boldsymbol{u}_{i}-$$

$$+\frac{1}{\gamma }\sum _{i=1}^b{l}_i^{\ast }({\mathbf{u}}_i)+\sum _{j=1}^d{h}_j^{\ast }({\mathbf{v}}_j).$$

$$(\forall \mathbf{v})g^{\ast }(\mathbf{v})\stackrel{\Delta }{=}\underset{\mathbf{z}}{\sup }{\mathbf{z}}^{†}\mathbf{v}-g(\mathbf{z}).$$

$$\mathbf{y}=\left[\begin{array}{c}{\mathbf{y}}_1\\ ⋮\\ {\mathbf{y}}_d\end{array}\right],\mathsf{\Phi }=\left[\begin{array}{c}{\mathsf{\Phi }}_1\\ ⋮\\ {\mathsf{\Phi }}_d\end{array}\right]=\left[\begin{array}{c}{\mathsf{\Theta }}_1{\mathsf{G}}_1{\mathsf{M}}_1\\ ⋮\\ {\mathsf{\Theta }}_d{\mathsf{G}}_d{\mathsf{M}}_d\end{array}\right]\mathsf{F}\mathsf{Z},$$

$${\mathbf{y}}_j={\mathsf{\Phi }}_j\mathbf{x}+{\mathbf{n}}_j,$$

$${\mathcal{D}}_e\left({\mathsf{W}}_i^{(k)}\right)=\frac{{\omega }^{(k)}}{{\omega }^{(k)}+{\left(\left|{\mathsf{\Psi }}_i^{†}{\mathbf{x}}^{(k-1)}\right|\right)}_e},$$

$$\Big{(}\operatorname{\boldsymbol{\mathcal{P}}}_{\mathcal{C}}(\boldsymbol{z})\Big{)}_{e}\overset{\Delta}{=}\left\{\begin{aligned} \Re(z_{e})&\qquad\Re(z_{e})>0\\ 0\quad&\qquad\Re(z_{e})\leq 0\end{aligned}\right.\quad\forall e$$

$$\Big{(}\operatorname{\boldsymbol{\mathcal{S}}}_{\alpha}(\boldsymbol{z})\Big{)}_{e}\overset{\Delta}{=}\left\{\begin{aligned} \displaystyle\frac{z_{e}\big{\{}|z_{e}|-\alpha\big{\}}_{+}}{|z_{e}|}&\qquad|z_{e}|>0\\ \displaystyle 0\qquad\quad&\qquad|z_{e}|=0\\ \end{aligned}\right.\quad\forall e,$$

$$P_{{\mathcal{B}}_j}(\mathbf{z})\stackrel{\Delta }{=}\{\begin{array}{rl}{ϵ}_j\frac{\mathbf{z}-{\mathbf{y}}_j}{\Vert \mathbf{z}-{\mathbf{y}}_j{\Vert }_2}+{\mathbf{y}}_{\mathbf{j}}& \Vert \mathbf{z}-{\mathbf{y}}_j{\Vert }_2>{ϵ}_j\\ \mathbf{z}& \Vert \mathbf{z}-{\mathbf{y}}_j{\Vert }_2\le {ϵ}_j.\end{array}$$

$${\mathrm{prox}}_g^{\mathsf{U}}(\mathbf{z})\stackrel{\Delta }{=}\underset{\bar{\mathbf{z}}}{\mathrm{argmin}}g(\bar{\mathbf{z}})+\frac{1}{2}(\mathbf{z}-\bar{\mathbf{z}}{)}^{†}\mathsf{U}(\mathbf{z}-\bar{\mathbf{z}}).$$

$$\operatorname{prox}^{\boldsymbol{\mathsf{U}}^{-1}}_{\alpha g^{*}}(\boldsymbol{z})=\Big{(}\operatorname{\boldsymbol{\mathcal{I}}}-\alpha\boldsymbol{\mathsf{U}}\operatorname{prox}^{\boldsymbol{\mathsf{U}}}_{\alpha^{-1}g}\Big{)}\big{(}\alpha^{-1}\boldsymbol{\mathsf{U}}^{-1}\boldsymbol{z}\big{)},$$

$${\mathrm{prox}}_{h_j}^{{\mathsf{U}}_j}(\mathbf{z})=\underset{\bar{\mathbf{z}}}{\mathrm{argmin}}{h}_j(\bar{\mathbf{z}})+\frac{1}{2}(\mathbf{z}-\bar{\mathbf{z}}{)}^{†}{\mathsf{U}}_j(\mathbf{z}-\bar{\mathbf{z}})$$

$$\displaystyle\;\>=\operatorname*{argmin}_{\bar{\boldsymbol{z}}}h_{j}(\bar{\boldsymbol{z}})+\frac{1}{2}\Big{(}\boldsymbol{\mathsf{U}}_{j}^{\frac{1}{2}}\boldsymbol{z}-\boldsymbol{\mathsf{U}}_{j}^{\frac{1}{2}}\bar{\boldsymbol{z}}\Big{)}^{\dagger}\Big{(}\boldsymbol{\mathsf{U}}_{j}^{\frac{1}{2}}\boldsymbol{z}-\boldsymbol{\mathsf{U}}_{j}^{\frac{1}{2}}\bar{\boldsymbol{z}}\Big{)}.$$

$$\displaystyle\operatorname{prox}^{\boldsymbol{\mathsf{U}}_{j}}_{h_{j}}\Big{(}\boldsymbol{\mathsf{U}}_{j}^{-\frac{1}{2}}\boldsymbol{s}\Big{)}=\boldsymbol{\mathsf{U}}_{j}^{-\frac{1}{2}}\bigg{(}\operatorname*{argmin}_{\bar{\boldsymbol{s}}}h_{j}\Big{(}\boldsymbol{\mathsf{U}}_{j}^{-\frac{1}{2}}\bar{\boldsymbol{s}}\Big{)}+$$

$$\displaystyle\qquad\qquad\qquad+\frac{1}{2}(\boldsymbol{s}-\bar{\boldsymbol{s}})^{\dagger}(\boldsymbol{s}-\bar{\boldsymbol{s}})\bigg{)}=\boldsymbol{\mathsf{U}}_{j}^{-\frac{1}{2}}\operatorname{\boldsymbol{\mathcal{P}}}_{\mathcal{E}_{j}}(\bar{\boldsymbol{s}}).$$

$$\underset{\mathbf{x}}{\min }f(\mathbf{x})+\gamma \sum _{i=1}^b{l}_i({\mathsf{W}}_i^{†}{\mathsf{\Psi }}_i^{†}\mathbf{x})+\sum _{j=1}^d{\tilde{h}}_j({\mathsf{G}}_j{\mathsf{M}}_j\mathsf{F}\mathsf{Z}\mathbf{x}),$$

$$\displaystyle\operatorname{prox}^{\boldsymbol{\mathsf{U}}_{j}}_{\tilde{h}_{j}}\big{(}\boldsymbol{\mathsf{\Theta}}_{j}^{-1}\boldsymbol{s}\big{)}=\boldsymbol{\mathsf{\Theta}}_{j}^{-1}\operatorname*{argmin}_{\bar{\boldsymbol{s}}}\tilde{h}_{j}\big{(}\boldsymbol{\mathsf{\Theta}}_{j}^{-1}\bar{\boldsymbol{s}}\big{)}+$$

$$\displaystyle\qquad\qquad\qquad\qquad+\frac{1}{2}\big{(}\boldsymbol{s}-\bar{\boldsymbol{s}}\big{)}^{\dagger}\boldsymbol{\mathsf{\Theta}}_{j}^{-1^{\dagger}}\boldsymbol{\mathsf{U}}_{j}\boldsymbol{\mathsf{\Theta}}_{j}^{-1}\big{(}\boldsymbol{s}-\bar{\boldsymbol{s}}\big{)}.$$

$$\displaystyle\operatorname{prox}^{\boldsymbol{\mathsf{U}}_{j}}_{\tilde{h}_{j}}\big{(}\boldsymbol{\mathsf{\Theta}}_{j}^{-1}\boldsymbol{s}\big{)}=$$

$$\displaystyle\qquad\qquad\quad\boldsymbol{\mathsf{\Theta}}_{j}^{-1}\operatorname*{argmin}_{\bar{\boldsymbol{s}}}h_{j}(\bar{\boldsymbol{s}})+\frac{1}{2}\big{(}\boldsymbol{s}-\bar{\boldsymbol{s}}\big{)}^{\dagger}\boldsymbol{\mathsf{D}}\big{(}\boldsymbol{s}-\bar{\boldsymbol{s}}\big{)},$$

$${\Vert {\left[\begin{array}{cc}\zeta \mathsf{I}& 0\\ 0& \eta \mathsf{U}\end{array}\right]}^{\frac{1}{2}}\left[\begin{array}{c}{\mathsf{W}}^{†}{\mathsf{\Psi }}^{†}\\ \mathsf{\Phi }\end{array}\right]{\left[\begin{array}{c}\tau \mathsf{I}\end{array}\right]}^{\frac{1}{2}}\Vert }_{\mathrm{S}}^2\le$$

$$\le \tau \zeta {\Vert {\mathsf{W}}^{†}{\mathsf{\Psi }}^{†}\Vert }_{\mathrm{S}}^2+\tau \eta {\Vert {\mathsf{U}}^{\frac{1}{2}}\mathsf{\Phi }\Vert }_{\mathrm{S}}^2<1,$$

$$\mathrm{SNR}=20{\log }_{10}\left(\frac{\Vert {\mathbf{x}}^{\circ }{\Vert }_2}{\Vert {\mathbf{x}}^{\circ }-{\mathbf{x}}^{(t)}{\Vert }_2}\right),$$

$$\mathrm{DR}=\frac{\sqrt{N}\Vert \mathsf{\Phi }{\Vert }_{\mathrm{S}}^2}{\Vert {\mathsf{\Phi }}^{†}(\mathbf{y}-\mathsf{\Phi }{\mathbf{x}}^{(t)}){\Vert }_2}\underset{e}{\max }{x}_e^{(t)}.$$

Full text

††pubyear: 2017††pagerange: An accelerated splitting algorithm for radio-interferometric imaging: when natural and uniform weighting meet–An accelerated splitting algorithm for radio-interferometric imaging: when natural and uniform weighting meet

1 Introduction

Radio-interferometry (RI) is a technique that permits the observation of radio emissions with great sensitivity and angular resolution. It provides valuable data for many research directions in astronomy, cosmology or astrophysics (Thompson et al., 2007). The next-generation radio telescopes, like the planned Square Kilometre Array (SKA; Dewdney et al., 2009), are expected to push the sensitivity further to achieve a dynamic range of six or seven orders of magnitude and to reconstruct large, giga-pixel size, images. To achieve such a feat, the amount of data to be acquired will be huge and the signal processing techniques from RI need to be revisited and reinvented. Fast specialised algorithmic solvers are being developed (Onose et al., 2016; Deguignet et al., 2016; Ferrari et al., 2014; Yatawatta, 2015, 2016; Carrillo et al., 2014) and vigorous research is being directed towards tackling the challenges of both RI imaging and RI calibration (Rau et al., 2009; Wijnholds et al., 2014).

The SKA, whose construction is scheduled to start in 2018, will be comprised of a huge number of antennas, approximately $131~{}000$ low frequency elements and 197 dishes for medium frequency (Dewdney et al., 2009; Broekema et al., 2015). With an expected number of $65~{}000$ frequency bands of operation, the data rates estimates will be in the terabits per second range (Broekema et al., 2015) and will present a challenge for both the communication infrastructure and signal processing. The current standard algorithmic solvers, belonging to the clean family (Högbom, 1974; Schwab, 1984; Bhatnagar & Cornwell, 2004; Cornwell, 2008), do not scale well to such tremendous data sizes.

Recently, convex optimisation techniques coupled with compressive sensing models (Wiaux et al., 2009a; Li et al., 2011; Carrillo et al., 2012; Garsden et al., 2015) have been shown to potentially outperform the standard state-of-the-art clean imaging algorithms. Such methods typically approach the imaging problem by minimising a convex objective function defined as a sum of multiple terms: several data terms dependent on the measured data (the visibilities), and a number of regularisation priors usually promoting sparsity or smoothness in an appropriate domain and positivity. This is a global approach, all algorithms searching for the unique solution that minimises the convex objective function.

Besides the reconstruction quality, the processing speed is of great interest with fast and parallelizable algorithms having been recently proposed (Carrillo et al., 2014; Ferrari et al., 2014; Yatawatta, 2015; Onose et al., 2016). Such approaches come in contrast with the standard clean methods which employ local procedures and rely on greedy updates and other signal pre-processing steps, like the RI weighting used to mitigate the effects produced by an unbalanced density profile of the sampling strategy. For algorithms that work directly in image space, like clean, the type of RI weighting is very important and affects the overall image reconstruction results (Briggs, 1995; Boone, 2013; Yatawatta, 2014). Natural weighting provides controlled noise statistics with the aim of maximising the sensitivity. Uniform weighting reduces the side-lobes of the point spread function by scaling the visibilities with the inverse sampling density and provides better resolution at the cost of lowered sensitivity. Since any weighting other than natural essentially biases the data, clean is not able to maximise both resolution and sensitivity. To mitigate this, intermediate robust weighting (Briggs, 1995) or adaptive weighting schemes (Yatawatta, 2014) have also been proposed and serve as a tradeoff between resolution and sensitivity.

Convex optimisation methods (Carrillo et al., 2012, 2014) that impose constraints directly in visibility space work with naturally weighting data. Such approaches can optimise both the resolution and sensitivity, which is impossible to achieve with clean and its evolutions. An unbalanced density profile of the sampling strategy does not influence the final solution of the convex optimisation problem. It can have however a potentially significant detrimental effect on the convergence speed of the algorithmic structures.

We study herein an acceleration strategy of the primal-dual (PD) algorithmic structure recently proposed by Onose et al. (2016). It can incorporate sampling density information into the algorithmic structure to achieve faster convergence speed for non uniform visibility distributions in $u$–$v$ space. We propose the use of a preconditioning strategy that improves the convergence speed significantly, making the PD approach even more appealing for the large-scale signal processing associated with the future radio telescopes. We rely on the same convex optimisation problem from Onose et al. (2016) but introduce a non-euclidian, skewed, proximity step that uses a preconditioning matrix reminiscent of the uniform weighting used by clean and the other RI imaging methods that work in image space. Intuitively, to link with the behaviour of clean, such an approach maintains the sensitivity of the natural weighting but achieves the resolution of the uniformly weighted data.

We show through simulations the acceleration achieved using the preconditioning strategy for simulated random, SKA and VLA coverages. A study of the computational burden of the non-euclidian proximity step is also included. We also showcase the reconstruction capabilities of the algorithm using real interferometric data of the 3C129 radio galaxy and compare with clean. The observations were performed for two 50 MHz channels using the VLA in configuration B and C.

The remainder of this article is organised as follows. Section 2 introduces the RI problem and briefly reviews the current existing standard solvers. Section 3 presents the main convex optimisation problem we associate with the image reconstruction and introduces the tools used by the preconditioned PD solver. Sections 4 details the proposed preconditioned PD algorithm and the acceleration strategy. Extensive simulations and results are presented in Section 5. Section 6 presents our final remarks and future work directions.

2 Radio-interferometric imaging

In RI, the measured data, the visibilities, are produced by an array of geographically separated antennas that are paired to measure radio emissions from a given area of the sky. Under the simplifying assumptions of non-polarised monochromatic RI imaging, the measurement equation for a measured visibility point $y(\boldsymbol{u})$ can be stated as

[TABLE]

with the direction dependent effects (DDEs) that affect the measurements, modelled through $D(\boldsymbol{l},\boldsymbol{u})$. Here, we denote by $\boldsymbol{u}=(u,v)$, the projected baseline components in the orthogonal plane relative to the line of sight. The observed sky brightness is described in the same coordinate system, with coordinates $(l,m)$. We denote $\boldsymbol{l}=(l,m)$. The well-known $w$ component effect, associated with the baseline components in the line of sight, is a known DDE. Unknown DDEs related to primary beam and ionospheric effects are assumed to have been properly calibrated so that we consider here a pure imaging problem.

The reconstruction algorithms work with a discretised version of the inverse problem (1). This resolves to the linear measurement equation

[TABLE]

where $\boldsymbol{x}\in\mathbb{R}^{N}$ is the unknown intensity image of interest of which $M$ visibility measurements $\boldsymbol{y}\in\mathbb{C}^{M}$ are taken by the radio telescope array. The measurements are corrupted by additive noise $\boldsymbol{n}$, each component $n_{e}$ assumed to have a known variance $\sigma=\sigma_{e},\forall e$. The measurement operator $\boldsymbol{\mathsf{\Phi}}=\boldsymbol{\mathsf{\Theta}}\boldsymbol{\mathsf{G}}\boldsymbol{\mathsf{F}}\boldsymbol{\mathsf{Z}}$ is a linear map from the image space to the visibility domain. It is composed of the matrix $\boldsymbol{\mathsf{G}}\in\mathbb{C}^{M\times nN}$ containing compact support interpolation kernels (Fessler & Sutton, 2003) and modeling the DDEs, an $n$-oversampled Fourier operator $\boldsymbol{\mathsf{F}}\in\mathbb{C}^{nN\times nN}$ and an oversampling and scaling operator $\boldsymbol{\mathsf{Z}}\in\mathbb{R}^{nN\times N}$ that pre-compensates for the interpolation (Fessler & Sutton, 2003). If the original visibilities are affected by noise with different variances, $\sigma_{e_{1}}\neq\sigma_{e_{2}}$ for some $e_{1}$ and $e_{2}$, a diagonal matrix $\boldsymbol{\mathsf{\Theta}}$ with diagonal elements $\theta_{e,e}=\frac{1}{\sigma_{e}}$ is used to whiten the noise. This is equivalent to the natural weighting performed in RI.

2.1 The CLEAN method

The inverse problem defined by (2) has been thoroughly studied and various deconvolution methods have been proposed. The standard imaging algorithms, belonging to the clean family, perform a greedy non-linear deconvolution based on local iterative beam removal (Högbom, 1974; Schwarz, 1978; Schwab, 1984; Thompson et al., 2007). They rely on a sparsity prior on the solution implicitly introduced through the greedy, pixel by pixel, image reconstruction procedure. This resembles the matching pursuit (MP) algorithm (Mallat & Zhang, 1993). It can also be seen as a regularised gradient descent method that minimises the residual norm $\|\boldsymbol{y}-\boldsymbol{\mathsf{\Phi}}\boldsymbol{x}\|_{2}^{2}$ via a gradient descent subject to an implicit sparsity constraint on $\boldsymbol{x}$ (Rau et al., 2009),

[TABLE]

The notation † denotes the adjoint of a linear operator. Multiple versions and improvements have been suggested, multi-scale clean (Cornwell, 2008), adaptive scale clean (Bhatnagar & Cornwell, 2004). In parallel with clean, the maximum entropy method solvers (Ables, 1974; Gull & Daniell, 1978; Cornwell & Evans, 1985) have been proposed but in practice clean was favoured.

2.2 Convex optimisation algorithms

Recently, convex optimisation methods are beginning to gain traction in RI and offer improved reconstruction quality and speed over the classical clean approaches (Wiaux et al., 2009a; Wiaux et al., 2009b; Wenger et al., 2010; Li et al., 2011; Carrillo et al., 2012, 2014; Ferrari et al., 2014; Yatawatta, 2015; Garsden et al., 2015; Dabbech et al., 2015; Onose et al., 2016). They approach the imaging problem under the framework of compressed sensing (CS). Such methods add a regularisation of the ill-posed reconstruction problem in the form of a prior that assumes a low dimensional signal model (Donoho, 2006; Candès et al., 2006). Seen through the CS framework, the signal of interest $\boldsymbol{x}$ is considered to have a sparse representation, $\boldsymbol{x}=\boldsymbol{\mathsf{\Psi}}\boldsymbol{\alpha}$ with $\boldsymbol{\alpha}\in\mathbb{C}^{D}$ containing only a few nonzero elements (Fornasier & Rauhut, 2011). The dictionary $\boldsymbol{\mathsf{\Psi}}\in\mathbb{C}^{N\times D}$ is usually a collection of wavelet bases or, more generally, an over-complete frame.

An analysis-based approach (Elad et al., 2007) to recover the signal of interest $\boldsymbol{x}$ by solving the ill-posed inverse problem (2) can be formally stated as (Carrillo et al., 2012; Carrillo et al., 2013; Carrillo et al., 2014; Onose et al., 2016)

[TABLE]

The sparsity averaging reweighed analysis (SARA) sparsity prior (Carrillo et al., 2012), used as the sparsity dictionary $\boldsymbol{\mathsf{\Psi}}$, has been shown to be a good sparsity basis. Since the solution $\boldsymbol{x}$ is an intensity image, a reality and positivity prior is also assumed. Data fidelity is enforced by constraining the residual to belong to an $\ell_{2}$ ball defined given an estimate $\epsilon$ of the noise affecting the measurements. Synthesis-based approaches have also been proposed (Wiaux et al., 2009a; Wiaux et al., 2009b; McEwen & Wiaux, 2011).

The $\ell_{0}$ norm is non-convex and thus the problem defined in (4) is intractable. By replacing the $\ell_{0}$ norm with its closest convex relaxation, the $\ell_{1}$ norm, and by reformulating the constraints from (4) with the use of the indicator function111 The indicator function $\iota_{\mathcal{C}}$ of a convex set $\mathcal{C}$ is defined as

$(\forall\boldsymbol{z})\quad\iota_{\mathcal{C}}(\boldsymbol{z})\overset{\Delta}{=}\left\{\begin{aligned} 0&\quad\boldsymbol{z}\in\mathcal{C}\\ +\infty&\quad\boldsymbol{z}\notin\mathcal{C}.\end{aligned}\right.$

$\iota_{\mathcal{C}}$ we can state a basic minimisation problem as

[TABLE]

The function $f=\iota_{\mathbb{R}^{N}_{+}}$ introduces the reality and positivity requirements for the recovered solution, the function $l=\|\cdot\|_{1}$ represents the sparsity inducing prior and $h(\boldsymbol{z})=\iota_{\mathcal{B}}(\boldsymbol{z}),\mathcal{B}=\{\boldsymbol{z}\in\mathbb{C}^{M}:\|\boldsymbol{z}-\boldsymbol{y}\|_{2}\leq\epsilon\}$ is the data fidelity term constraining the residual to be situated in the $\ell_{2}$ ball $\mathcal{B}$ defined by the noise level $\epsilon$. A re-weighted $\ell_{1}$ approach (Candès et al., 2008) is generally used to approximate the $\ell_{0}$ norm by imposing the weights $\boldsymbol{\mathsf{W}}$ on the operator $\boldsymbol{\mathsf{\Psi}}$ and solving sequentially several $\ell_{1}$ problems with different $\boldsymbol{\mathsf{W}}$. This basic minimisation problem (Carrillo et al., 2012; Onose et al., 2016) has been approached using several state-of-the-art algorithmic solvers: the simultaneous direction method of multipliers (Carrillo et al., 2014), the alternating direction method of multipliers and a PD algorithm with forward-backward iterations (Onose et al., 2016).

The forward-backward iterative structure is one of the main pillars used in the algorithmic structure presented herein. We can view it as being conceptually extremely close to the major-minor cycle structure of clean. Consider one of the most basic approaches, the unconstrained version of the minimisation problem (4), namely $\min_{\boldsymbol{x}}\|\boldsymbol{\mathsf{W}}^{\dagger}\boldsymbol{\mathsf{\Psi}}^{\dagger}\boldsymbol{x}\|_{1}+\rho\|\boldsymbol{y}-\boldsymbol{\mathsf{\Phi}}\boldsymbol{x}\|^{2}_{2}$, with $\rho$ a free parameter. This can be solved using forward-backward iterations by performing a gradient step together with a proximal step (Moreau, 1965),

[TABLE]

The forward gradient step consists in doing a step in the opposite direction to the gradient of the $\ell_{2}$ norm of the residual. This is essentially equivalent to a major cycle of clean. In this particular case, the proximal step is a simple soft-thresholding operation in the given basis $\boldsymbol{\mathsf{W}}^{\dagger}\boldsymbol{\mathsf{\Psi}}^{\dagger}$ (Combettes & Pesquet, 2007). It consists in decreasing the absolute values of all the coefficients of $\boldsymbol{\mathsf{W}}^{\dagger}\boldsymbol{\mathsf{\Psi}}^{\dagger}\boldsymbol{x}$ that are above a certain threshold by the threshold value, and setting to zero those below the threshold. Such an approach is very similar to the minor cycle of clean, with the soft-threshold value being an analogous to the loop gain factor. clean iteratively builds up the signal by picking up parts of the most important coefficients until the residuals become negligible. The soft-thresholding acts globally by removing small and insignificant coefficients, on all signal locations simultaneously. As such, clean can be intuitively understood as a very specific version of the forward-backward algorithm.

3 Forward-backward PD algorithm

We continue by reviewing the minimisation problem and the randomised PD algorithm (Condat, 2013; Vũ, 2013; Pesquet & Repetti, 2015) recently proposed for RI by Onose et al. (2016), on which this work relies. It solves a primal, block wise, minimisation problem similar to (5),

[TABLE]

together with its dual formulation (Bauschke & Combettes, 2011),

[TABLE]

Here, since the $\ell_{1}$ norm is additively separable, we have split the over-complete sparsity basis into $b$ parts, $\boldsymbol{\mathsf{\Psi}}=\begin{bmatrix}\boldsymbol{\mathsf{\Psi}}_{1}&\ldots&\boldsymbol{\mathsf{\Psi}}_{b}\end{bmatrix}$. The weighting matrix $\boldsymbol{\mathsf{W}}$ is also split to produce a weight matrix $\boldsymbol{\mathsf{W}}_{i}$ for each $\boldsymbol{\mathsf{\Psi}}_{i}$. The scalar $\gamma$ is a free configuration parameter and only affects the convergence speed (Onose et al., 2016). The functions from (7) are defined block wise but similarly to (5). Thus, the functions $l_{i}=\|\cdot\|_{1}$ represent the sparsity inducing prior and $h_{j}(\boldsymbol{z})=\iota_{\mathcal{B}_{j}}(\boldsymbol{z}),\mathcal{B}_{j}=\{\boldsymbol{z}\in\mathbb{C}^{M_{j}}:\|\boldsymbol{z}-\boldsymbol{y}_{j}\|_{2}\leq\epsilon_{j}\}$ are the data fidelity terms constraining the residual to be situated in $\ell_{2}$ balls defined by the noise level $\epsilon_{j}$, for each part of the visibility data $\boldsymbol{y}_{j}$. The notation ∗ denotes the Legendre-Fenchel conjugate function.222 The Legendre-Fenchel conjugate function $g^{*}$ of a function $g$ is

$(\forall\boldsymbol{v})\qquad g^{*}(\boldsymbol{v})\overset{\Delta}{=}\sup_{\boldsymbol{z}}\boldsymbol{z}^{\dagger}\boldsymbol{v}-g(\boldsymbol{z}).$

3.1 Distributed problem formulation

We work in a setup where the visibility data are split into $d$ blocks, such that

[TABLE]

to allow for distributed and parallelised processing (Carrillo et al., 2014; Onose et al., 2016). We also rely on the fact that $\boldsymbol{\mathsf{G}}$ is composed of compact support kernels and introduce the matrices $\boldsymbol{\mathsf{M}}_{j}\in\mathbb{R}^{nN_{j}\times nN}$ to select only the parts of the discrete Fourier plane involved in computations for block $j$. Each block operator $\boldsymbol{\mathsf{G}}_{j}\in\mathbb{C}^{M_{j}\times nN_{j}}$ requires partial Fourier information, namely only $nN_{j}$ coefficients (Onose et al., 2016). The diagonal matrix $\boldsymbol{\mathsf{\Theta}}$ is also split accordingly.

The inverse problem (2) was therefore be rewritten for each data block as

[TABLE]

with $\boldsymbol{n}_{j}$ being the part of the noise associated with the measurements $\boldsymbol{y}_{j}$ and with $\boldsymbol{\mathsf{\Phi}}_{j}$ the associated linear operator.

3.2 The re-weighted $\ell_{1}$ approach

A re-weighted $\ell_{1}$ (Candès et al., 2008) serves to approximate the $\ell_{0}$ norm by solving successive $\ell_{1}$ penalised problems. The weights $\boldsymbol{\mathsf{W}}^{(k)}_{i}$, at step $k$, are computed based on the solution $\boldsymbol{x}^{(k-1)}$ from the previously solved problem from step $k-1$ such that

[TABLE]

with the operator $\boldsymbol{\mathcal{D}}_{e}$ denoting diagonal element $e$. The parameter $\omega^{(k)}$ is decreased from a preset value at each re-weight step. This ensures that, after several such steps, if the values of the $e$th coefficient $\left(\left|\boldsymbol{\mathsf{\Psi}}^{\dagger}_{i}\boldsymbol{x}^{(k)}\right|\right)_{e}$ are large, the penalty applied is decreased towards [math]. The small coefficients, smaller than $\omega^{(k)}$, are still being largely penalised. Thus, this iterative procedure removes the bias introduced by the $\ell_{1}$ relaxation of the sparsity constraint. This procedure is summarised as Algorithm 1. Note that each call to Algorithm 3, which will be detailed in the following sections, should use the past primal and dual solutions, from step $k-1$, as initialisation in order to warm start the convergence.

3.3 Proximity operators

As previously mentioned, the PD algorithm (Pesquet & Repetti, 2015) relies on forward-backward iterations (Komodakis & Pesquet, 2015) to deal with the non smooth terms present in both the primal minimisation problem (7) and its dual formulation (8). The forward step corresponds to a gradient-like step and the backward step is an implicit sub-gradient-like step performed through the use of the proximity operator (Moreau, 1965).

Using the definition (6), the proximity operator associated with the function $f$ in (7) has a closed form solution and becomes the projection

[TABLE]

onto the positive real orthant. Similarly, the proximity operator for the sparsity prior functions $l_{i}$ is the component wise soft-thresholding operator

[TABLE]

for a given threshold $\alpha$. For the data fidelity terms $h_{j}$, the proximity operator has a closed form as the projection onto an $\ell_{2}$ ball $\mathcal{B}_{j}$,

[TABLE]

More details can be found in Onose et al. (2016), which proposed the PD algorithm for solving (7) and (8) in the absence of any preconditioning strategy.

4 Accelerated forward-backward PD algorithm

The structure of the proposed algorithm, presented in Algorithm 3, is based on Pesquet & Repetti (2015). It is similar to that of the PD algorithm proposed for RI by Onose et al. (2016). As before, we solve concurrently both the primal minimisation problem (7) and its dual formulation (8). Forward-backward iterations, consisting of a gradient descent step coupled with a proximal update, are used to update both the primal and the dual variables. The key difference that accelerates the convergence speed is the use a new non-Euclidean proximity operator for the data fidelity to replaces the projection onto the $\ell_{2}$ ball, used in Onose et al. (2016), with a projection onto a generalised ellipsoid that incorporates both the noise statistics and sampling density information. By incorporating the sampling density information, the algorithm can make a larger step towards the final solution at each iteration. This acceleration strategy changes only the forward-backward step associated with the data fidelity terms, the rest of the updates remain the same as in Onose et al. (2016). In analogy with clean, the algorithm can be understood as being composed of complex clean-like forward-backward steps performed in parallel in multiple data, prior and image spaces Onose et al. (2016).

4.1 Non-euclidean proximity operator

A generalisation of the proximity operator allows us to use additional prior information about the data when performing the computations associated with the data fidelity terms $h_{j}$, in order to accelerate the convergence speed. It offers a broad flexibility in the way the data fidelity is enforced throughout the iterations.

Thus, we rely on the generalised proximity operator relative to a metric induced by a strongly positive, self-adjoint333A linear operator $\boldsymbol{\mathsf{U}}$ is said to be strongly positive and self-adjoint if $\langle\boldsymbol{x}|\boldsymbol{\mathsf{U}}\boldsymbol{x}\rangle\geq\alpha\|\boldsymbol{x}\|^{2}_{2},\forall\boldsymbol{x},\forall\alpha>0$ and $\boldsymbol{\mathsf{U}}^{\dagger}=\boldsymbol{\mathsf{U}}$, respectively. linear operator $\boldsymbol{\mathsf{U}}$ (Hiriart-Urruty & Lemarechal, 1996),

[TABLE]

The standard definition from (6) is found when $\boldsymbol{\mathsf{U}}=\boldsymbol{\mathsf{I}}$. A generalisation of the Moreau decomposition provides the link between the proximity operators of a function $g$ and that of its conjugate $g^{*}$ (Combettes & Vũ, 2014; Pesquet & Repetti, 2015) for any operator $\boldsymbol{\mathsf{U}}$,

[TABLE]

and allows for a facile way of computing the proximity operators for the conjugate functions.

We choose the preconditioning matrices $\boldsymbol{\mathsf{U}}_{j}$ to be diagonal, with positive, non-zero diagonal elements and thus positive definite and invertible. It results directly from (15) that

[TABLE]

By making the variable change $\boldsymbol{s}=\boldsymbol{\mathsf{U}}_{j}^{\frac{1}{2}}\boldsymbol{z}$ and $\bar{\boldsymbol{s}}=\boldsymbol{\mathsf{U}}_{j}^{\frac{1}{2}}\bar{\boldsymbol{z}}$ we can rewrite (17) as

[TABLE]

Here we have denoted by $\operatorname{\boldsymbol{\mathcal{P}}}_{\mathcal{E}_{j}}$ the projection onto a generalised ellipsoid $\mathcal{E}_{j}=\{\bar{\boldsymbol{s}}\in\mathbb{C}^{M_{j}}:\|\boldsymbol{\mathsf{U}}_{j}^{-\frac{1}{2}}\bar{\boldsymbol{s}}-\boldsymbol{y}_{j}\|_{2}\leq\epsilon_{j}\}$ associated with the preconditioned matrix $\boldsymbol{\mathsf{U}}_{j}$ and the data fidelity function $h_{j}$. This formulation serves as a generalisation of the way data fidelity is enforced (Carrillo et al., 2014; Onose et al., 2016). Note that the minimisation problem (7) and its dual formulation (8) do not change when the generalised proximity operator (15) is used. This only affects the way convergence is achieved. Thus, if $\boldsymbol{\mathsf{U}}_{j}=\boldsymbol{\mathsf{I}}$, the constraints that the residual should belong to the $\ell_{2}$ balls $\mathcal{B}_{j}$ is enforced such that the Euclidian distance from the starting point $\boldsymbol{\mathsf{\Phi}}_{j}\boldsymbol{x}$ and the ball $\mathcal{B}_{j}$ is minimised. This results in the simple projection onto the $\ell_{2}$ ball $\mathcal{B}_{j}$ from (14). If instead a different metric $\boldsymbol{\mathsf{U}}_{j}\neq\boldsymbol{\mathsf{I}}$ is used, the projection becomes skewed and the Euclidian distance to the ball $\mathcal{B}_{j}$ is not minimised anymore. However, the new projection point still satisfies $\|\boldsymbol{\mathsf{\Phi}}_{j}\boldsymbol{x}-\boldsymbol{y}_{j}\|_{2}\leq\epsilon_{j}$. This can be expressed as the projection onto the ellipsoid $\mathcal{E}_{j}$ with the resulting projection point moved to the $\ell_{2}$ ball by the application of $\boldsymbol{\mathsf{U}}_{j}^{-\frac{1}{2}}$ in (18).

For a generic metric $\boldsymbol{\mathsf{U}}_{j}\neq\boldsymbol{\mathsf{I}}$, an iterative procedure is required to compute the proximity operator $\operatorname{prox}^{\boldsymbol{\mathsf{U}}_{j}}_{h_{j}}$. We propose a forward-backward approach that works directly with the definition of the proximity step (17). It performs a gradient step, with step $\mu$, in the direction of the smooth term $\frac{1}{2}(\boldsymbol{z}-\bar{\boldsymbol{z}})^{\dagger}\boldsymbol{\mathsf{U}}_{j}(\boldsymbol{z}-\bar{\boldsymbol{z}})$ followed by the application of the proximity operator for the function $h_{j}$, which is the projection (14). This is formally presented as Algorithm 2. The step size $\mu$ must satisfy $\mu\leq\nicefrac{{1}}{{\|\boldsymbol{\mathsf{U}}_{j}\|^{2}_{\mathrm{S}}}}$. Since the preconditioning matrix $\boldsymbol{\mathsf{U}}_{j}$ is diagonal, we have \|\boldsymbol{\mathsf{U}}_{j}\|_{\mathrm{S}}=\max_{e}\big{(}\boldsymbol{\mathcal{D}}_{e}(\boldsymbol{\mathsf{U}}_{j})\big{)} with the operator $\boldsymbol{\mathcal{D}}_{e}$ selecting the $e$th diagonal element of $\boldsymbol{\mathsf{U}}_{j}$.

Faster converging proximal gradient algorithms for solving (15) may be employed (Tseng, 2008). However, for simplicity we limit the presentation herein to the forward-backward approach presented as Algorithm 2. Alternatively, we can compute the projection $\operatorname{\boldsymbol{\mathcal{P}}}_{\mathcal{E}_{j}}$ onto the ellipsoid $\mathcal{E}_{j}$ and then estimate $\operatorname{prox}^{\boldsymbol{\mathsf{U}}_{j}}_{h_{j}}(\boldsymbol{z})$ as in (18). A very fast iterative approach was developed by Dai (2006) for any choice of metric $\boldsymbol{\mathsf{U}}_{j}$. It requires an initial point on the feasible region, which, due to $\boldsymbol{\mathsf{U}}_{j}$ being positive definite and invertible, can be easily computed using Algorithm 2. Note that this is not the case for a general operator $\boldsymbol{\mathsf{U}}_{j}$, for which the derivations form (17) and (18) are not guaranteed to hold.

4.2 The preconditioned algorithmic structure

All the updates associated with the dual variables $\boldsymbol{v}^{(t)}_{j}$ and $\boldsymbol{u}^{(t)}_{i}$ from (8) are performed in Algorithm 3 in parallel in steps $9$–$13$ and $18$–$22$, respectively. Randomisation is supported given a probabilistic construction of the active sets $\mathcal{P}$ and $\mathcal{D}$. Thus, only a part of the dual variables are updated per iteration, the rest remaining unchanged as in steps $14$–$17$ and $23$–$26$. The forward-backward updates rely on the Moreau decomposition (16) to compute the proximity operator associated with the conjugate functions $l_{i}^{*}$ and $h_{j}^{*}$ relying on the proximity operator of the functions $l_{i}$ and $h_{j}$. The resulting updates become the soft-thresholding (13) for the prior dual variables $\boldsymbol{u}^{(t)}_{i}$ from step $19$ and the skewed projection (18) onto the ellipsoid $\mathcal{E}_{j}$ for the data fidelity dual variables $\boldsymbol{v}^{(t)}_{j}$ from step $10$. For the soft-thresholding, we perform a re-parametrisation similar to the one performed in Onose et al. (2016). Since $\gamma$ is a free parameter, we replace the resulting algorithmic soft-threshold size $\frac{\gamma}{\zeta}$ with $\kappa\|\boldsymbol{\mathsf{\Psi}}\boldsymbol{\mathsf{W}}\|^{2}_{\rm{S}}$ to produce a operator independent configuration parameter $\kappa$. The parameter $\kappa$ is only linked to the scale of the unknown image to be recovered. The application of the operators $\boldsymbol{\mathsf{G}}_{j}^{\dagger}\boldsymbol{\mathsf{\Theta}}^{\dagger}_{j}$ and $\boldsymbol{\mathsf{\Psi}}_{i}\boldsymbol{\mathsf{W}}_{i}$ is also performed in parallel, in steps $12$ and $21$. The contribution of all the dual variables is then used to update the primal variable, the image of interest $\boldsymbol{x}^{(t)}$ in steps $28$–$29$. This is a forward-backward step which, through the use of the Moreau decomposition, resumes to the projection (12) onto the positive orthant presented in step $28$.

4.3 The epiphany: when natural and uniform weighting meet

For the data fidelity terms $h_{j}$ we propose the use of a non-trivial invertible preconditioning matrix $\boldsymbol{\mathsf{U}}_{j}$ which has links to the standard weighting schemes. The weighting is used to mitigate the effects produced by the sampling strategy (Briggs, 1995; Yatawatta, 2014) and serves as an important pre-processing step for the clean family of algorithms. We aim to incorporate the sampling density information into the PD algorithmic structure, through $\boldsymbol{\mathsf{U}}_{j}$, while solving the same problems defined in (7) and (8). This does not change the overall results due to the convergence guarantees of the convex optimisation methods and increases the speed of convergence, as will be shown through simulations.

With this aim, we employ a diagonal preconditioning matrix $\boldsymbol{\mathsf{U}}_{j}$, for each visibility block $\boldsymbol{y}_{j}$. The matrix $\boldsymbol{\mathsf{U}}_{j}$ accounts for the sampling density similarly to the uniform weighting. It contains on the diagonal the inverse of the sampling density in the vicinity of each associated visibility point. This has the benefit of allowing for a facile computation of its inverse which is important to the computational complexity of the resulting strategy. Other types of preconditioning could also be supported.

To give further insight into the behaviour of this preconditioning strategy, consider the problem (7) written in an equivalent formulation

[TABLE]

by introducing the natural weighting matrix $\boldsymbol{\mathsf{\Theta}}_{j}$ in the definition of the function $\tilde{h}_{j}(\boldsymbol{z})=\iota_{\tilde{\mathcal{E}}_{j}}(\boldsymbol{z}),\tilde{\mathcal{E}}_{j}=\{\boldsymbol{z}\in\mathbb{C}^{M_{j}}:\|\boldsymbol{\mathsf{\Theta}}_{j}\boldsymbol{z}-\boldsymbol{y}_{j}\|_{2}\leq\epsilon_{j}\}$. Now, the convex set associated with $\tilde{h}_{j}$ becomes the ellipsoid $\tilde{\mathcal{E}}$ associated with the natural weight matrix $\boldsymbol{\mathsf{\Theta}}_{j}$. This does not change the definition of the minimisation problems but changes significantly how the problem is approached algorithmically. It changes the manner in which the data fidelity constraint is enforced to make it similar to the way the generalised proximity operator is used in the algorithm. As such, it allows us to provide an intuitive link between the whitening matrices $\boldsymbol{\mathsf{\Theta}}_{j}$ and the preconditioning matrices $\boldsymbol{\mathsf{U}}_{j}$ by highlighting that they enter the algorithmic structure through a similar mechanism.

Thus, based on the definition of the proximity operator (15) and by performing the variable change $\boldsymbol{s}=\boldsymbol{\mathsf{\Theta}}_{j}\boldsymbol{z}$ and $\bar{\boldsymbol{s}}=\boldsymbol{\mathsf{\Theta}}_{j}\bar{\boldsymbol{z}}$, we can write $\operatorname{prox}^{\boldsymbol{\mathsf{U}}_{j}}_{\tilde{h}_{j}}(\boldsymbol{z})$ as

[TABLE]

Since both $\boldsymbol{\mathsf{\Theta}}_{j}$ and $\boldsymbol{\mathsf{U}}_{j}$ are diagonal matrices and since \tilde{h}_{j}\big{(}\boldsymbol{\mathsf{\Theta}}_{j}^{-1}\boldsymbol{s}\big{)}=h_{j}(\boldsymbol{s}), (20) becomes

[TABLE]

with diagonal elements $d_{e,e}=\sigma_{e}^{2}\boldsymbol{\mathcal{D}}_{e}(\boldsymbol{\mathsf{U}}_{j})$. The operator $\boldsymbol{\mathcal{D}}_{e}$ selects the $e$th diagonal element from $\boldsymbol{\mathsf{U}}_{j}$. Since they affect the data fidelity term $h_{j}$ in a similar way, this provides an intuitive link between the natural weighting matrix $\boldsymbol{\mathsf{\Theta}}_{j}$ and the preconditioning matrix $\boldsymbol{\mathsf{U}}_{j}$, which is based on the inverse of the sampling density. A large value for $d_{e,e}$ corresponds to either a low sample density for the frequency vicinity of the given measurement $e$ or a large noise variance for the same measurement. Low values $d_{e,e}$ correspond to less noisy measurements or a high sampling density. Since sampling the same $u$–$v$ region multiple times can be seen as lowering the noise by averaging the data, the similitude between the effect of the noise on the measurement and the sampling density is immediate.

Let us emphasise again that only the natural weighting performed through $\boldsymbol{\mathsf{\Theta}}_{j}$ is reflected back into the definition of the minimisation problem through the application of $\boldsymbol{\mathsf{\Theta}}_{j}^{-1}$ in (21). In contrast, the preconditioning matrix is only an internal algorithmic flexibility to solve the very same problem. Thus, such an approach can be seen to incorporate all the benefits from both natural and uniform weighting in clean terms. On one hand it optimises resolution by accounting for the correct noise statistics, leveraging natural weighting in the definition of the minimisation problem for image reconstruction. On the other hand it optimises sensitivity by enabling accelerated convergence through a preconditioning strategy incorporating sampling density information à la uniform weighting.

4.4 Convergence requirements

The variables $\boldsymbol{x}^{(t)}$, $\boldsymbol{v}_{j}^{(t)}$ and $\boldsymbol{u}_{i}^{(t)}$, $\forall i,j$, are guaranteed to converge to the solution of the PD problem (7)–(8) for an adequately chosen set of configuration parameters, $\tau$, $\zeta$ and $\eta$. The convergence conditions (Pesquet & Repetti, 2015, Lemma 4.3) can be stated explicitly for Algorithm 3 as

[TABLE]

with the use of the triangle and Cauchy-Schwarz inequalities and with the diagonal matrices $\boldsymbol{\mathsf{I}}$ of a proper dimension. The matrix $\boldsymbol{\mathsf{U}}$ represents a diagonal concatenation of all the preconditioning matrices $\boldsymbol{\mathsf{U}}_{j}$ associated with the differently split operators and data blocks. A relaxation with the factor $0<\lambda\leq 1$ of the updates is also permitted. The additional parameter $\gamma>0$ imposes that $\kappa>0$ as well. For the randomised setup, the probabilities with which the active sets $\mathcal{P}$ and $\mathcal{D}$ are generated have to be nonzero and the activated variables need to be drawn in an independent and identical manner along the iterations.

The general framework of the PD with forward-backward iterations approach and its mathematical analysis are presented by Pesquet & Repetti (2015).

4.5 Computational complexity

The complexity and parallelised and distributed implementation details follow closely the study from Onose et al. (2016). The only difference is the introduction of the preconditioning matrix and the need for the iterative computation of the resulting proximity operator. The complexity class of Algorithm 2 is $\mathcal{O}(M_{j})$ per data block $j$. The computations involving the projection are to be performed in a distributed fashion similarly to the computations involving the data fidelity terms. The convergence speed of Algorithm 2 is linked to the conditioning number of the preconditioning matrix and may slow down for ill-conditioned matrices. In such case, Algorithm $4$ proposed by Dai (2006) or faster proximal gradient methods (Tseng, 2008) become preferable. Empirical evidence however suggests that the accuracy of the projection can be lowered by reducing the number of iteration performed without damaging the convergence speed of the whole algorithm. The algorithm is resilient to errors in the computations and in practice as little as $1$ iteration can be enough to achieve a significant acceleration. This can serve to control the added complexity due to the sub-iterative computation of the preconditioned proximity operator. Comparing the added total computational complexity of the preconditioning, which is $\mathcal{O}(M)$ per sub-iteration, with that of the basic non-preconditioned PD algorithm, which is of the order $\mathcal{O}(nN\log nN)+\mathcal{O}(dN)+\mathcal{O}(MN)$ per iteration, it is evident that the added cost due to the preconditioning in PPD is negligible when the number of sub-iterations is kept small.

For more details regarding the complexity, randomisation and general structure of the PD algorithm solving (7) and (8) we direct the reader to Onose et al. (2016).

5 Simulations and results

We study the acceleration for different sampling strategies of the $u$–$v$ space. To judge the efficacy of the acceleration, we compare the preconditioned algorithm PPD against the non-preconditioned PD and ADMM algorithms (Onose et al., 2016), solving the same minimisation problem. We also compare the reconstruction quality and acceleration using real interferometric measurement of the 3C129 radio galaxy. In this case, we showcase the reconstruction in comparison with clean, as implemented by the wsclean package (Offringa et al., 2014). We provide reconstruction for multi-scale clean, denoted as MS-CLEAN. We do not study the distribution and randomisation, an extensive study being performed by Onose et al. (2016).

We work with pre-calibrated measurements, for both simulated and real data. We assume the absence of DDEs and a small field of view such that the measurement operator is a Fourier operator. We used an oversampled factor $n=4$ and a matrix $\boldsymbol{\mathsf{G}}$ that performs an interpolation of the frequency data using $8\times 8$ Kaiser-Bessel interpolation kernels (Fessler & Sutton, 2003) to average nearby uniformly distributed frequency. The diagonal preconditioning matrix $\boldsymbol{\mathsf{U}}$ contains the inverse of the sampling density as diagonal elements.

Thus, we begin by performing synthetic tests with the $u$–$v$ space sampled using a zero-mean, generalised Gaussian distribution (Novey et al., 2010) with shape parameter $\beta$. This allows us to have control of the sampling densities and see how the preconditioning is able to accelerate the convergence speed for various sampling patterns. We also use realistic simulations of VLA and SKA coverages and we study, through simulations, the behaviour of the algorithms. The $u$–$v$ coverages used are included in Figure 2. For all these tests we use two test images to generate the visibilities, namely a $477\times 1024$ image of the Cygnus A radio galaxy and a $512\times 512$ simulated image of a galaxy cluster with faint extended emission, respectively. The galaxy cluster image was produced using the faraday tool (Murgia et al., 2004). The two images are presented in Figure 5. The simulated visibilities are corrupted by zero-mean complex independent Gaussian noise. We run simulations for two noise levels, to produce an input signal to noise ratio $\rm{iSNR}=30~{}\rm{dB}$ and $\rm{iSNR}=50~{}\rm{dB}$ on the visibilities, respectively. This is accomplished by choosing the appropriate noise power relative to the power of the simulated, noise free, signal. In this case, the resulting noise statistics are used to generate the weight matrix $\boldsymbol{\mathsf{\Theta}}$.

For the comparison with clean we rely on observations of the 3C129 radio galaxy: right ascension 04h 45m 31.695s, declination 44∘ 55’ 19.95”, J2000. The observations were performed using the VLA for two 50 MHz channels centred at 4.59 and 4.89 GHz on the 25th of July 1994 in configuration B and 3rd of November 1994 in configuration C, respectively. The calibration and flagging for radio frequency interference have been performed in Pratley et al. (2016) according to the casa manual. We additionally remove approximatively $20~{}000$ visibility points that contained large noise outliers, probably visibilities affected by radio frequency interference or poorly calibrated. The remaining data consist of $307~{}780$ visibilities. The normalised $u$–$v$ coverage is also included in Figure 2. All reconstructions are performed at twice the resolution of the telescope array. This is necessary to avoid tension between the band limitation of the reconstructed image and the positivity constraint introduced by our approach.

For the synthetic tests, we assess the reconstruction performance in terms of the signal to noise ratio,

[TABLE]

where $\boldsymbol{x}^{\circ}$ is the original image and $\boldsymbol{x}^{(t)}$ is the reconstructed estimate of the original. For the real data reconstructions, since we do not have access to the ground truth, we report the dynamic range obtained for the reconstruction,

[TABLE]

5.1 Choice of parameters

The PPD algorithms converge given that (22) is satisfied. To ensure this we set $\zeta=\nicefrac{{1}}{{\|\boldsymbol{\mathsf{\Psi}}\boldsymbol{\mathsf{W}}\|_{\rm{S}}^{2}}}$, $\eta=\nicefrac{{1}}{{\|\boldsymbol{\mathsf{U}}^{\frac{1}{2}}\boldsymbol{\mathsf{\Phi}}\|_{\rm{S}}^{2}}}$ and $\tau=0.49$. The relaxation parameter is set to 1. For the ADMM and PD algorithms we set the parameters as recommended by Onose et al. (2016). We do not use randomisation, all data and all sparsity priors are used at each iteration. We use the SARA collection of wavelets (Carrillo et al., 2012), namely a concatenation of a Dirac basis with the first eight Daubechies wavelets, as sparsity prior. For the simulations, we set the normalised soft-threshold values $\kappa=10^{-4}$ for all three methods, PPD, PD and ADMM. We run PPD for a number of sub-iteration $n_{\rm itr}\in\{1,5,50\}$. In all tests we impose that the square of the global bound $\epsilon^{2}$ is 2 standard deviations above the mean of the $\chi^{2}$ distribution associated with the noise (Onose et al., 2016).

For the real data reconstruction we set $\kappa=10^{-6}$, since the recovered image has the brightest pixel on the order of $10^{-2}$. In this case we also perform $10$ re-weighting steps, one every $1024$ iterations, according to Algorithm 1. We start with $\omega^{(0)}=10^{-2}$ and set $\omega^{(k)}=0.25^{k}\omega^{(0)}$ for each step $k$. In this case the global bound $\epsilon^{2}$ is set to be $1.05$ times mean of the $\chi^{2}$ distribution associated with the thermal noise affecting the visibilities. Such a bound was observed to provide good reconstruction results. MS-CLEAN was run using the wsclean software package, version $2.2.1$, with both uniform and natural weighting. For both weighting schemes, we use 6 scales, $\{0,16,24,32,48,64\}$. We set the major loop gain to $\gamma_{M}=0.6$ and the minor loop gain to $\gamma_{m}=0.08$. The stopping threshold is set to 2 standard deviations above the automatically estimated noise level on the different scales. The uniform weighting test reached the stopping threshold. The natural weighting test was stopped after $35~{}000$ iterations since, for a larger number of iterations, the method was only accumulating spurious components without improving the solution.

5.2 Simulations

To study the behaviour of PPD across a broad range of $u$–$v$ sampling strategies, we use coverages with the sampled $u$–$v$ points distributed according to a generalised Gaussian distribution with the shape parameter $\beta$. We study the acceleration for the reconstruction of the galaxy cluster test image in Figure 3 and for the reconstruction of the Cygnus A test image in Figure 4. Here, we report the evolution of the ${\rm SNR}$ as a function of the number of iterations. In both cases we have performed tests for two levels of input noise, $30~{}\rm{dB}$ and $50~{}\rm{dB}$. For all test cases we provide the distribution of the normalised $u$ and $v$ coordinates to showcase the link between the convergence speed and sampling pattern.

For sampling strategies that are farther away from uniform, the preconditioning strategy improves the convergence rate dramatically in all test cases. For a Gaussian sampling, when $\beta=2$, the converge speed of the PPD is similar to that of PD and ADMM. A decrease in $\beta$ does not affect PPD greatly. It maintains almost the same convergence speed throughout all the test cases. In the extreme case when $\beta=0.25$, the density of measurements is much greater in the centre of the $u$–$v$ space and PPD becomes one order of magnitude faster than PD and ADMM. In all test cases, the PPD algorithm remains robust to an inexact computation of the ellipsoid projection. In practice there is little difference between performing $1$ sub-iteration and performing as many as $50$. Due to this, its complexity per iteration is marginally larger than that of PD. This, coupled with the improved convergence rate, makes PPD much more suitable for the large-scale problems arising in RI. Comparing the two input noise regimes, for lower input noise, the gap between PPD and PD becomes larger. For less noisy data, the sampling density becomes the most important factor that limits the convergence speed. This is due to the high frequency data having lower power than the low frequency data. For large noise, the high frequency visibilities are below the noise level and the effective coverage can be considered to be truncated at the point where the data are overwhelmed by the noise. For the low noise setup, the algorithms can improve the reconstruction and achieve a higher $\rm SNR$ but the coverage becomes more important for the convergence speed because the effective useful visibilities cover a wider range of frequencies in the $u$–$v$ space.

To further validate the behaviour of the algorithms, we also study them for the reconstruction of the two test images using simulated, but realistic SKA and VLA coverages. The evolution of the ${\rm SNR}$ a s function of iteration number for these test cases is presented in Figure 5. In all tests PPD maintains a similar level of acceleration as observed before, for the generalised Gaussian distributed $u$–$v$ coverages. For the SKA coverages, where the conditioning number of the preconditioning matrix is large, the number of sub-iterations begins to affect the evolution of PPD. Especially of the Cygnus A image it seems that using only one sub-iteration is actually faster. This behaviour is probably due to the fact that the preconditioning matrix is not optimal. Performing only one sub-iteration can be understood as projection onto a slightly different ellipsoid.

Figures 6 and 7 contain the reconstructed images for PPD and PD at iteration $99$ for the galaxy cluster image with VLA coverage and the Cygnus A image with the SKA coverage, respectively. The reconstruction quality achieved by PPD at this iteration is evident. Such a reconstruction is possible with PD only by performing approximatively $10$ times more iterations. The figures also contain embedded an animation that cycles through the iterations and shows the solution estimates at each iteration.666The animation is only supported when the PDF file is opened using Adobe Acrobat Reader, https://get.adobe.com/reader/ The evolution of PPD resembles a behaviour that is associated with the uniform weighting used for clean while the evolution of PD resembles that associated with natural weighting. Both methods however converge towards the same global solution, the solution of the natural weighted data. The sampling density information is only incorporated into PPD to accelerate the convergence speed.

5.3 Real data reconstruction

For the real data scenario, we study the reconstruction quality of PPD in comparison with MS-CLEAN using observations of the 3C129 radio galaxy performed with the VLA. The reconstructed images are illustrated in $\log_{10}$ scale in Figure 8. We note that PPD achieves better quality in terms of both resolution and sensitivity. It is able to better recover the faint emissions towards the tail of the main emission and has very little noise incorporated in the image. In comparison, MS-CLEAN includes multiple spurious components in the model map and due to the post processing achieves a poor resolution, especially around the main bright source that generates the two emission plumes. The resolution is much worse when the natural weighting is used since the size of the clean beam used is larger. The clean model is also lower resolution than in the uniform weighting case.

To better visualise the reconstruction quality, we provide in all images enlarged sections of the main source in the two boxes on the left and of the fainter point sources, from the lower part of the recovered image, in the two boxes on the right. The faint emission showcased enlarged in the right, upper box for the PPD reconstruction is most likely the source $\rm C$ reported by Lane et al. (2002). This is the faintest emission PPD can detect without introducing noise and deconvolution artefacts. Note that this source, as well as the emission tail of 3C129 are around $2.5$ orders of magnitude fainter than the brightest source. MS-CLEAN is unable to recover these emissions well and has brighter spurious components around the main emission. Setting the deconvolution threshold lower for MS-CLEAN, in order to extract more of the signal from the measurements, greatly increases the amount of spurious components detected.

As a last figure, we present the evolution of the $\rm DR$ for PPD and PD as a function of the number of iterations in Figure 9. This serves to validate the acceleration also on real data. Here, PPD is shown to be faster than PD. The distribution of the $u$–$v$ coordinates, also reported in Figure 9, is not that extreme in this case and the speed up is small, of the order of $1.5$, which is consistent with the previous simulations. This test serves to prove that the preconditioning works on real data. For more unbalanced sampling profiles we expect a larger acceleration, as demonstrated through simulations. Also, since the number of sub-iterations performed by PPD to approximate the preconditioned proximity operator is small, $n_{\rm itr}=5$, the complexity per iteration is similar to that of PD.

6 Conclusions

We proposed an acceleration of the PD algorithmic framework for solving the RI imaging problem. Building on the highly parallelizable structure of the PD algorithm, the accelerated PPD algorithm, benefits from all the flexibility of the PD, allowing for an efficient distributed implementation, by using full splitting and randomised updates. The analogy between the clean major-minor loop and the forward-backward iterations used by the method, can portray PPD as being composed of sophisticated clean-like iterations running in parallel in multiple data, prior, and image spaces.

The proposed approach reconciles natural and uniform weighting of clean algorithms. It optimises resolution by accounting for the correct noise statistics, leveraging natural weighting in the definition of the minimisation problem for image reconstruction. It also optimises sensitivity by enabling accelerated convergence through a preconditioning strategy incorporating sampling density information à la uniform weighting.

We study the acceleration through extensive simulations with realistic $u$–$v$ coverages and using real visibilities from the observation of the 3C129 radio galaxy with the VLA. The preconditioning strategy is able to increase the convergence speed by up to an order of magnitude for highly non-uniformly sampled coverages. We also showcase the reconstruction quality in comparison with MS-CLEAN for this data, exemplifying the improved resolution and sensitivity the PPD method offers.

Our Matlab code is available online on GitHub, http://basp-group.github.io/ppd-for-ri/. In the near future we intend to provide an efficient implementation in the purify c++ package for a distributed computing infrastructure.

Acknowledgements

This work was supported by the UK Engineering and Physical Sciences Research Council (EPSRC, grants EP/M011089/1 and EP/M008843/1). We would like to thank Federica Govoni and Matteo Murgia for providing the simulated galaxy cluster image.