A parallel & automatically tuned algorithm for multispectral image deconvolution

TL;DR

This paper introduces MUFFIN, a scalable, automatically tuned parallel algorithm for multispectral radio image deconvolution, designed to handle large data sets like those from the future SKA telescope.

Contribution

The paper presents MUFFIN, a novel parallel and automatically tuned wideband deconvolution algorithm for radio interferometry, suitable for big data applications.

Findings

MUFFIN outperforms existing algorithms in reconstruction quality.

MUFFIN is scalable and suitable for large-scale data like SKA.

Automatic tuning enhances usability and performance.

Abstract

In the era of big data, radio astronomical image reconstruction algorithms are challenged to estimate clean images given limited computing resources and time. This article is driven by the need for large scale image reconstruction for the future Square Kilometre Array (SKA), which will become in the next decades the largest low and intermediate frequency radio telescope in the world. This work proposes a scalable wideband deconvolution algorithm called MUFFIN, which stands for "MUlti Frequency image reconstruction For radio INterferometry". MUFFIN estimates the sky images in various frequency bands given the corresponding dirty images and point spread functions. The reconstruction is achieved by minimizing a data fidelity term and joint spatial and spectral sparse analysis regularization terms. It is consequently non-parametric w.r.t. the spectral behaviour of radio sources. MUFFIN algorithm is endowed with a parallel implementation and an automatic tuning of the regularization parameters, making it scalable and well suited for big data applications such as SKA. Comparisons between MUFFIN and the state-of-the-art wideband reconstruction algorithm are provided.

Figures

Equations

$${\mathbf{y}}_{\ell }={\mathbf{H}}_{\ell }{\mathbf{x}}_{\ell }^{\star }+{\mathbf{w}}_{\ell },\ell =1,\dots ,L$$

$${\text{argmin}}_{\mathbf{X}}\sum _{\ell =1}^L\frac{1}{2}\Vert {\mathbf{y}}_l-{\mathbf{H}}_{\ell }{\mathbf{x}}_{\ell }{\Vert }^2+1_{{\mathbb{R}}^{+}}(\mathbf{X})+h(\mathbf{X}),$$

$$h(\mathbf{X}):=\sum _{\ell =1}^L{\mu }_{\ell }\Vert {\mathbf{W}}_s{\mathbf{x}}_{\ell }{\Vert }_1+\sum _{n=1}^N{\gamma }_n\Vert {\mathbf{W}}_{\lambda }{\mathbf{x}}^n{\Vert }_1.$$

$$\left\{\begin{array}[]{ll}\gamma_{n}&=\gamma\alpha_{n},\quad n=1,\ldots,N\\ \mu_{\ell}&=\mu\beta_{\ell},\quad\ell=1,\ldots,L\\ \end{array}\right.$$

$$\left\{\begin{array}[]{ll}\alpha_{n}&=\frac{1}{\|\boldsymbol{W}_{\lambda}(\boldsymbol{x}^{n})_{0}\|_{2}+\varepsilon}\\ \beta_{\ell}&=\frac{1}{\|\boldsymbol{W}_{s}(\boldsymbol{x}_{\ell})_{0}\|_{2}+\varepsilon}\\ \end{array}\right.$$

$$(u)_{+}:=\left\{\begin{array}[]{ll}u&\quad\text{if}\quad u>0\\ 0&\quad\text{otherwise.}\\ \end{array}\right.$$

$$\text{sat}(u):=\left\{\begin{array}[]{ll}-1&\quad\text{if}\quad u<-1\\ 1&\quad\text{if}\quad u>1\\ u&\quad\text{if}\quad|u|\leq 1.\end{array}\right.$$

$$\tau \left(\frac{\beta }{2}+\sigma \left(\sum _{\ell =1}^L{\mu }_{\ell }^2\Vert {\mathbf{W}}_s{\Vert }^2+\sum _{n=1}^N{\gamma }_n^2\Vert {\mathbf{W}}_{\lambda }{\Vert }^2\right)\right)<1,$$

$$\nabla$$

$$\tilde{\mathbf{X}}$$

$$\mathbf{U}$$

$$\mathbf{V}$$

$$\mathbf{X}$$

$${\mathbf{x}}_{\ell }$$

$${\nabla }_{\ell }$$

$${\mathbf{s}}_{\ell }$$

$${\mathbf{m}}_{\ell }$$

$${\tilde{\mathbf{x}}}_{\ell }$$

$${\mathbf{p}}_{\ell }$$

$${\mathbf{u}}_{\ell }$$

$${\tilde{\mathbf{v}}}^n$$

$${\mathbf{v}}^n$$

$${\mathbf{t}}^n$$

$$\mathsf{R}=\sum _{\ell =1}^L{\mathsf{R}}_{\ell },$$

$${\mathsf{R}}_{\ell }=\mathbb{E}\Vert {\mathbf{H}}_{\ell }({\mathbf{x}}_{\ell }^{\star }-{\mathbf{x}}_{\ell }){\Vert }^2.$$

$$\begin{gathered} \hat{\mathsf{R}}=\sum _{\ell =1}^L{\hat{\mathsf{R}}}_{\ell } \\ . \end{gathered}$$

$${\hat{\mathsf{R}}}_{\ell }=\Vert {\mathbf{y}}_{\ell }-{\mathbf{d}}_{\ell }{\Vert }^2-N{\sigma }_w^2+{\sigma }_w^2\frac{2}{ϵ}⟨{\overline{\mathbf{d}}}_{\ell }-{\mathbf{d}}_{\ell },\text{mn@boldsymbol}{\delta }_{\ell }⟩$$

$${\text{lim}}_{ϵ\to 0}\mathbb{E}\left[\hat{\mathsf{R}}\right]=\mathsf{R}.$$

$${\hat{\mathsf{R}}}_{\ell }^{\prime }=2({\mathbf{H}}_{\ell }{\mathbf{x}}_{\ell }-{\mathbf{y}}_{\ell }{)}^{\top }{\mathbf{H}}_{\ell }{\mathbf{x}}_{\ell }^{\prime }+\frac{2{\sigma }_w^2}{ϵ}⟨{\overline{\mathbf{d}}}_{\ell }^{\prime }-{\mathbf{d}}_{\ell }^{\prime },\text{mn@boldsymbol}{\delta }_{\ell }⟩,$$

$${\text{lim}}_{ϵ\to 0}\mathbb{E}\left[{\hat{\mathsf{R}}}^{\prime }\right]=\mathbb{E}\left[{\mathsf{R}}^{\prime }\right],$$

$$\mathcal{U}(u):=\left\{\begin{array}[]{ll}0&\quad\text{if}\quad u\leq 0,\\ 1&\quad\text{if}\quad u>0,\\ \end{array}\right.$$

$$\Pi(u):=\left\{\begin{array}[]{ll}1&\quad\text{if}\quad-1\leq u\leq 1,\\ 0&\quad\text{if}\quad\text{otherwise},\\ \end{array}\right.$$

$${\mathbf{\theta }}^{\star }={\text{argmin}}_{\mathbf{\theta }}\hat{\mathsf{R}},$$

$${\mathbf{\theta }}^{(i+1)}={\mathbf{\theta }}^{(i)}-\eta \nabla {\hat{\mathsf{R}}}_{0:n_c}({\mathbf{\theta }}^{(i)}),$$

$$\text{mn@boldsymbol}{\theta }^{(q+1)}=\text{mn@boldsymbol}{\theta }^{(q)}-\eta \nabla {\hat{\mathsf{R}}}_{q\Delta :k}(\text{mn@boldsymbol}{\theta }^{(q)}),$$

$${\mathbf{x}}_{\ell }^{\prime }$$

Full text

††pubyear: 2019††pagerange: A parallel & automatically tuned algorithm for multispectral image deconvolution–A parallel & automatically tuned algorithm for multispectral image deconvolution

1 Introduction

The new generation of radio interferometers is characterized by very wide fields of view, high resolution, huge sensitivity and very large bandwidth. The future Square Kilometre Array (SKA) is emblematic of these extreme instruments. In its final phase of construction, SKA will consist of up to one million dipole antennas and two thousands parabolic antennas (Dewdney et al., 2009, 2013). These antennas, comprising a collecting area of approximately one square kilometer, will be located in Australia and South Africa. The SKA radio telescope is expected to provide images of the sky in the giga pixel range with (sub-)arcsec angular resolution, fields of view up to 200 $\text{deg}^{2}$ and more sensitivity than any other present-day instrument with a dynamic range covering up to seven orders of magnitude (see e.g. Table 1 in Dewdney et al. (2013)).

The work presented here addresses the large-scale and high-precision radio-interferometric image reconstruction problem. Radio-interferometric data, the visibilities, are the cross-correlations between voltages from pairs of antennas. In the ideal case, two receivers with baseline $\boldsymbol{b}$ observing in a narrow frequency band $c/\lambda$ measure a complex visibility equal to the Fourier spectrum of the intensity distribution of interest at spatial frequency $\boldsymbol{b}/\lambda$. The partial measurement of the Fourier plane is what makes the image reconstruction strongly ill posed. In practice, the measurement model corresponds to a linear map from the image domain to the visibility space, but includes multiple distortions such as: a) direction-dependent effects (DDE) related to wide-field effects, non-coplanarity of the baselines, the ionospheric perturbations and associated Faraday rotation; b) antenna level effects such as the primary beam patterns, element-based gains and c) computational approximation of operators such as the degridding required to use Fast Fourier Transform algorithms, see for example (Wijnholds et al., 2010; Rau et al., 2009; Bhatnagar et al., 2008). Compensations of the perturbations a) and b) rely on alternating a (self-)calibration step and an image reconstruction step. This paper concentrates on the image reconstruction step.

The imaging task is realized by solving either directly the inverse problem related to the calibrated visibility measurements, or the inverse problem in the image domain after applying the adjoint measurement model to the visibilities. The latter approach relies on the construction of a ‘dirty image’ from calibrated visibilities. The dirty image is related to the sky image by a convolution if the DDEs are properly compensated. This is the approach followed in the classical major/minor loop (Jongerius et al., 2014; Tasse et al., 2013) and in the facet frameworks (Kogan & Greisen, 2009; van Weeren et al., 2016); see (Tasse et al., 2018) for a detailed comparison between these two techniques. The main advantage of the image domain processing is not only the computing efficiency related to the fact that the direct and adjoint operators are simple convolutions, but also the amount of data to process: the dirty image size equals the reconstructed sky image size while the number of raw visibilities can be several order of magnitudes the number of pixels in the sky image.

Interestingly, sparsity was early recognized as a powerful tool for radio-interferometric image reconstruction and has lead to the most populated family of imaging algorithms. Their patriarch is the CLEAN algorithm (Högbom, 1974), which expresses and exploits the sparsity of the sky intensity distribution in the canonical basis. Variants of CLEAN as Cornwell (2008) generalize this approach to extended objects. More recently, various algorithms based on the theory of compressed sensing (Donoho, 2006; Candès & Wakin, 2008) and convex optimization emerged for the reconstruction of radio interferometric images from the calibrated visibilities. These algorithms exploit the fact that the image has a sparse representation in some appropriate transform domain through the use of sparsifying regularizations; see for example (Wiaux et al., 2009; McEwen & Wiaux, 2011; Kartik et al., 2017). A greedy alternative to the minimization of a regularized cost function was proposed by Dabbech et al. (2015). Sparse image representations using wavelet dictionaries has indisputably shown substantial improvement with respect to classical CLEAN-based approaches, as demonstrated by Garsden et al. (2015) using LOFAR data and Dabbech et al. (2018) using VLA data. The price to pay is an increased computational cost. As a result, a lot of effort has gone into designing scalable methods able to cope with large data sets. E.g., the authors of (Onose et al., 2016b, a) achieve parallel computations by splitting the visibilities into multiple blocks and exploiting randomization over data blocks at each iteration.

The SKA will observe the sky in wide sub-bands (covering all together from $50$ MHz to $14$ GHz) and the reconstruction of both spatial and spectral behavior of radio sources is essential to characterize the astrophysical origin of the detected radiation, see e.g. Kraus (1986). As a result, multi frequency reconstruction algorithms are required in order to fully exploit the spectral and spatial resolution of the SKA telescope. Taking the path of CLEAN, the Multi-Frequency Synthesis (MFS) (Conway et al., 1990) and the Multi-Scale Multi-Frequency Synthesis (MS-MFS) (Rau & Cornwell, 2011) deconvolution algorithms were the first wide-band (WB) radio-interferometric reconstruction algorithms. At the time of writing, still relatively few such WB reconstruction algorithms exist (Bajkova & Pushkarev, 2011; Wenger & Magnor, 2014; Junklewitz et al., 2015; Abdulaziz et al., 2016).

In Rau & Cornwell (2011); Bajkova & Pushkarev (2011); Junklewitz et al. (2015) the spectrum of each pixel is modeled as a (generalized) power law, or its Taylor expansion. These ‘semi-parametric’ methods rely on spectral models and thus clearly offer many advantages such as estimation accuracy when the model is appropriate. However, across the broad frequency coverage of current radio facilities, radio sources exhibiting complex spectral shapes (not simple power laws) are expected. For instance, Kellermann (1974) evidences that some sources may exhibit one or more relative minima, breaks and turnovers. For the new generation of low frequency telescopes such as LOFAR, recent studies have also shown that second order broadband spectral models are often insufficient (Scaife & Heald, 2012). In addition to that, not assuming any spectral shape for radio sources in the deconvolution step can be useful for mutli-wavelength polarization studies, as well as for simultaneously identifying spectral lines within wide-band continuum observations. However, few WB reconstruction algorithms consider the general case of reconstructing the spatio-spectral ‘image cube’ without constraining the spectral variation to a power law. In particular, few WB reconstruction algorithms have opted instead for a sparse regularization of the spectra. Among these, Wenger & Magnor (2014) formulate the problem as an inverse problem with a smooth spectral regularization allowing for local deviations; Abdulaziz et al. (2016) constrain the sky image cube to be low rank in order to enforce that only few sources, each with a distinct spectral signature, can exist in the field of view. This approach requires a singular value decomposition of the whole sky image cube at each iteration, so computational costs can become an issue. Finally, WB image reconstruction has also recently been tackled as a joint blind source separation and reconstruction problem by Jiang et al. (2017).

The present article focuses on total intensity wide-band and HI imaging task within the minor loop. Section 2 derives the proposed algorithm which is nonparametric in the sense that it does not rely on a particular assumption of the spectral structure, but rather on a sparse representation of the spectra in a dictionary. The spectral dimension critically blows up the size of the inverse problem and a particular effort is dedicated to the scalability of the algorithm. This results in an algorithm where most computationally intensive steps, such as the spatial wavelets decompositions and the associated adjoint operators, can be computed in parallel for each wavelength. Given the execution time for a large scale reconstruction, the selection of the regularization parameters is particularly critical. Section 3 proposes an automatic tuning strategy based on the Stein Unbiased Risk Estimator. Finally, section 4 presents simulation results which demonstrate the effectiveness of the algorithm.

It is worthy to note that, whereas simulations in section 4 are realized in the deconvolution framework, the proposed algorithm can be used in the same way to reconstruct the sky image cube from dirty images or visibilities, depending on the direct operator.

2 Inverse problem

In this section, we first formulate the multi-frequency image deconvolution task as a regularized and constrained least squares optimization problem. We then propose and discuss a primal-dual algorithm for solving the aforementioned optimization problem.

2.1 Optimization problem

Let $\boldsymbol{x}_{\ell}{\bf{{}^{\star}}}$ and $\boldsymbol{y}_{\ell}$ be $N\times 1$ column vectors collecting the unknown (and sought-after) sky intensity image and the (measured) dirty intensity image respectively. The subscript $\ell$ in the notations $\boldsymbol{x}_{\ell}{\bf{{}^{\star}}}$, $\boldsymbol{y}_{\ell}$ and in all following matrix and vector notations refers to the wavelength band $\lambda_{\ell}$, with $\ell=1,\ldots,L$ and $L$ the number of spectral channels. At each wavelength band, the dirty image $\boldsymbol{y}_{\ell}$ is related to the sky image $\boldsymbol{x}_{\ell}{\bf{{}^{\star}}}$ by a linear mapping and contaminated with noise. This takes the following form:

[TABLE]

where $\boldsymbol{H}_{\ell}$ represents a convolution by the point spread function (PSF), and $\boldsymbol{w}_{\ell}$ is a perturbation vector accounting for noise and modeling error. The perturbation vector $\boldsymbol{w}_{\ell}$ is assumed to be Gaussian, independent and identically distributed (i.i.d) with zero mean. Without loss of generality we will assume that the variance of the noise does not depend on the wavelength and is equal to $\sigma^{2}_{w}$. More sophisticated correlation structures can be taken into account provided that the covariance matrix is known. Given the dirty images and the PSF, Eq. (1) defines an ill-posed deconvolution problem. The problem is ill-posed owing to the partial coverage of the Fourier plane, i.e. due to the fact that $\boldsymbol{H}_{\ell}$ is rank deficient. A meaningful solution can be obtained in a cost minimization framework by incorporating regularizations and constraints reflecting our prior knowledge on the unknown variables $\boldsymbol{x}_{1}{\bf{{}^{\star}}},\ldots,\boldsymbol{x}_{L}{\bf{{}^{\star}}}$, the true sky intensity images. Let $\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{L}$ denote the optimization variables corresponding to $\boldsymbol{x}_{1}{\bf{{}^{\star}}},\ldots,\boldsymbol{x}_{L}{\bf{{}^{\star}}}$, and let $\boldsymbol{X}$ denote the $N\times L$ matrix $\boldsymbol{X}:=[\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{L}]$ obtained by stacking the image vectors next to each other. Let $\boldsymbol{x}^{1},\ldots,\boldsymbol{x}^{N}$ denote the rows of $\boldsymbol{X}$ which represent the pixels spectra. Hence, $\boldsymbol{X}$ can be alternatively obtained by stacking the pixels spectra on top of each other, i.e. $\boldsymbol{X}:=[{\boldsymbol{x}^{1}}^{\top},\ldots,{\boldsymbol{x}^{N}}^{\top}]^{\top}$. With these notations, we propose to minimize the following cost function:

[TABLE]

with

[TABLE]

The first term in the cost function in equation (2) is a data fidelity term penalizing the discrepancy between the measurements and the estimated model. The second term in equation (2) is a positivity constraint enforcing each element in $\boldsymbol{X}$ to be positive, more precisely $\boldsymbol{1}_{\mathbb{R}^{+}}(\boldsymbol{X})=0$ if $\boldsymbol{X}\geqslant 0$ (where the inequality is element-wise) and $\boldsymbol{1}_{\mathbb{R}^{+}}(\boldsymbol{X})=\infty$ otherwise. The third term in equation (2), developed in equation (3), is a spatial-spectral regularization. More precisely, $h(\boldsymbol{X})$ is the sum of a spatial regularization imposed on the images at all the wavelength bands and a spectral regularization imposed on the spectra at all the pixels. The matrices $\boldsymbol{W}_{s}$ and $\boldsymbol{W}_{\lambda}$ represent the spatial and spectral decomposition operators respectively. The scalars $\mu_{\ell}\in\mathbb{R}^{+}$ for $\ell=1,\ldots L$ and $\gamma_{n}\in\mathbb{R}^{+}$ for $n=1,\ldots N$ are the spatial and the spectral regularization weights respectively.

The spatial and spectral regularizations in equation (3) are sparse analysis regularizations (Elad et al., 2007). The sparse analysis approach consists in penalizing the $\ell_{1}$-norm of the correlations between the atoms in a dictionary, i.e. $\boldsymbol{W}_{s}$ and $\boldsymbol{W}_{\lambda}$ in the case of (3), and the images, i.e. $\boldsymbol{x}_{\ell}$ for $\ell=1,\ldots,L$ and the spectra, i.e. $\boldsymbol{x}^{n}$ for $n=1,\ldots,N$. Penalizing the $\ell_{1}$-norm of these decompositions promotes a sparse correlation between the atoms in the dictionary and the variables. Hence, it encourages finding a solution that is not correlated to some atoms of the dictionary and conversely more correlated to others. The dictionaries can be orthonormal transforms or more generally redundant transforms also referred to as over-complete dictionaries. Possible choices for $\boldsymbol{W}_{s}$ are (concatenations of) Daubechies wavelets (Carrillo et al., 2014; Onose et al., 2016b) or the isotropic undecimated wavelet transform (IUWT) (Garsden et al., 2015; Starck & Bobin, 2010; Dabbech et al., 2015) that have proved their efficiency for various image processing tasks in astronomy. One possible choice for $\boldsymbol{W}_{\lambda}$ that will be illustrated in the experiments is the Discrete Cosine Transform (DCT), a Fourier related transform. The DCT spectral regularization assumes that the DCT is compact, an assumption widely used in the sparse reconstruction literature.

Other types of regularizations based on sparse representations, such as synthesis priors or more generally representations based on 3D spatio-spectral dictionaries, can be of course considered. See (Lanusse et al., 2014) for a general discussion on 3D wavelets and e.g. (Paris et al., 2013) where a spatio-spectral dictionary is derived for a specific detection task. The drawback of these regularizations is computational: contrarily to (3), they couple space and frequency, i.e. the lines and columns of $\boldsymbol{X}$, and prevent an efficient implementation of (2), e.g. with parallel steps. This point will be clarified in the next section. Finally, it is essential to note that if $\forall\ell$, $\gamma_{\ell}=0$, the minimization problem in Eq. (2) becomes separable w.r.t. $\boldsymbol{x}_{\ell}$ and reduces to $L$ deconvolution problems that can be solved independently at each wavelength $\lambda_{\ell}$.

The regularization parameters $\mu_{\ell}$ for $\ell=1,\ldots,L$ and $\gamma_{n}$ control the strength of the regularization term w.r.t. the data fidelity term. A large value of $\mu_{\ell}$ tends to over-smooth structures at wavelength $\ell$, while a small value of $\mu_{\ell}$ leads to under-regularization. Similarly, a large value of $\gamma_{n}$ tends to over-smooth spectra at pixel $n$. See, e.g. (Renard et al., 2011) in the case of a single regularization parameter for monochromatic optical interferometry. Although the use of a different regularization parameter per pixel and band could improve the estimation performance compared to the use of one regularization parameter per regularization type, the tuning of $(L+N)$ regularization parameters is in practice very challenging. For this reason, we propose to use two regularization parameters and assign a different (fixed) weight per pixel and per band. By doing this, we reduce the number of regularization parameters from $(L+N)$ to $2$ and we take into account the spatial and spectral variability of the images and spectra respectively. We adopt the following model for the regularization parameters:

[TABLE]

with

[TABLE]

where $(\boldsymbol{x}^{n})_{0}$ and $(\boldsymbol{x}_{\ell})_{0}$ correspond to the spectra and images of the sky image used for initializing the algorithm. The intuition for the weights setting is similar to the one in re-weighted $\ell_{1}$ minimization schemes (Candès et al., 2008) where the strength of the regularization parameter is inversely proportional to the norm of the corresponding regularized norm. The intuition for this setting is that the value of the regularization parameter should be high whenever the decomposition coefficient is close to zero in order to encourage sparsity, and it should be small otherwise. Nevertheless, in re-weighted $\ell_{1}$ minimization schemes, the weights may be updated at each iteration or after a certain number of iterations using the current estimates for $\boldsymbol{x}^{n}$ and $\boldsymbol{x}_{\ell}$. In this work, and in the simulations, we consider fixed weights computed from a coarse initial estimation of $\boldsymbol{X}$.

Finally, note that the proposed optimization problem in equation (2) is convex but does not admit a closed form solution. In order to find the minimizer of (2), we need to resort to an iterative procedure as described in the next section.

2.2 Optimization algorithm

The proposed optimization problem described in equations (2) and (3) is the sum of a smooth and differentiable function (the LS term), and three non-smooth and non-differentiable functions (the indicator function, the sparse analysis spectral regularization, and the sparse analysis spatial regularization). Nevertheless, the non-differentiable functions, the indicator function and the two sparse analysis regularizations, are simple functions in the sense that their proximity operators (Combettes & Pesquet, 2011) have closed form representations. In order to solve the proposed optimization problem, we resort to the primal dual splitting algorithm proposed by Vũ (2011); Condat (2014) also known as the Vũ-Condat algorithm. The Vũ-Condat algorithm is well suited for optimization problems similar to the one proposed in this work i.e. problems that consist of a sum between a smooth function and several proximable functions involving linear operators. Note that other splitting algorithms could have been used such as the Alternating direction Method of Multipliers ADMM (Boyd et al., 2011), Simultaneous Direction of Multipliers SDMM (Setzer et al., 2010; Combettes & Pesquet, 2011) or the primal dual splitting algorithm proposed by Combettes & Pesquet (2012). These algorithms operate by splitting the objective function to be minimized into smaller and simpler functions that are dealt with individually. However, the Vũ-Condat algorithm is more appealing compared to ADMM and SDMM because it does not involve any implicit inversion of the operators $\boldsymbol{I}+\alpha\boldsymbol{W}_{\lambda}^{\top}\boldsymbol{W}_{\lambda}$ and $\boldsymbol{I}+\beta\boldsymbol{W}_{s}^{\top}\boldsymbol{W}_{s}$ which can be very large and cumbersome to invert, see for example the work by Meillier et al. (2018). It can be considered as a generalization of Chambolle-Pock algorithm (Chambolle & Pock, 2010) which includes a differentiable term and multiple analysis regularizations.

Vũ-Condat algorithm is not symmetric with respect to the primal and dual variables. Switching the roles of these variables gives rise to two different versions. A straightforward application of the first Vũ-Condat algorithm (Condat, 2014, Algorithm 3.1) to the optimization problem described in equations (2) and (3) with no relaxation ($\rho_{n}=1$) leads to algorithm 1. These choices are discussed at the end of this subsection. We refer to algorithm 1 as MUFFIN for MUlti Frequency image reconstruction For radio INterferometry. Algorithm 1 requires using the proximity operator of the positivity indicator function, namely the projection on the positive orthant defined as:

[TABLE]

The proximity operator of the $\ell_{1}$-norm is the well-known soft thresholding operator. Algorithm 1 requires using the proximity operator of the Fenchel conjugate of the $\ell_{1}$-norm function which can be easily computed from the soft thresholding operator using Moreau’s identity. This leads to the saturation function defined as:

[TABLE]

Note that both $(\cdot)_{+}$ and $\text{sat}(\cdot)$ used in algorithm 1 are applied element wise when the input is a matrix. The matrices $\boldsymbol{D}_{\mn@boldsymbol{\mu}}$ and $\boldsymbol{D}_{\mn@boldsymbol{\gamma}}$ used in algorithm 1 are diagonal matrices such that $\boldsymbol{D}_{\mn@boldsymbol{\mu}}:=\text{Diag}(\mu_{1},\ldots,\mu_{L})$, and $\boldsymbol{D}_{\mn@boldsymbol{\gamma}}:=\text{Diag}(\gamma_{1},\ldots,\gamma_{N})$, and the parameters $\sigma>0$ and $\tau>0$ are parameters of the algorithm to be fixed appropriately in order to guarantee convergence. More precisely, they must be fixed such that:

[TABLE]

where $\beta$ is the Lipschitz constant of the gradient of the LS term in equation (2).

Hereafter, we briefly describe each step in Algorithm 1. In equation (9), the gradient of the LS term is computed. This consists of applying the linear operator $\boldsymbol{H}_{\ell}^{\dagger}$ on the residue of the model for $\ell=1,\ldots,L$. In equation (10), $\boldsymbol{\tilde{X}}$ is computed. The adjoint operators $\boldsymbol{W}_{s}^{\dagger}$ and $\boldsymbol{W}_{\lambda}^{\dagger}$ are applied on the dual variables $\boldsymbol{U}$ and $\boldsymbol{V}^{\top}$. More precisely, the third term is where the operator $\boldsymbol{W}_{s}^{\dagger}$ is applied on $\boldsymbol{U}$, and the fourth term is where the operator $\boldsymbol{W}_{\lambda}^{\dagger}$ is applied on $\boldsymbol{V}^{\top}$ ($\boldsymbol{V}\boldsymbol{W}_{\lambda}=(\boldsymbol{W}_{\lambda}^{\dagger}\boldsymbol{V}^{\top})^{\top}$). The primal variable $\boldsymbol{X}$ is updated using the gradient and the aforementioned inverse transforms, and a projection on the positive orthant is then applied. The result is then assigned to the intermediate primal variable $\boldsymbol{\tilde{X}}$. In equation (11), the first dual variable $\boldsymbol{U}$ is computed. The linear operator $\boldsymbol{W}_{s}$ is applied on $(2\boldsymbol{\tilde{X}}-\boldsymbol{X})$. The dual variable $\boldsymbol{U}$ is updated using previously computed transform, and the saturation function is then applied. The result is assigned to the dual variable $\boldsymbol{U}$. In equation (12), similarly to the previous step, the second dual variable $\boldsymbol{V}$ is computed. The linear operator $\boldsymbol{W}_{\lambda}$ is applied on $(2\boldsymbol{\tilde{X}}-\boldsymbol{X})^{\top}$. The dual variable $\boldsymbol{V}$ is updated using previously computed transform,, and the saturation function is then applied. The result is assigned to the dual variable $\boldsymbol{V}$. In equation (13), $\boldsymbol{X}$ is updated by simply assigning it to $\boldsymbol{\tilde{X}}$. Finally, note how the proposed iterative algorithm 1 activates separately: the gradient of the LS term, i.e. the linear operator $\boldsymbol{H}_{\ell}^{\dagger}$ for $\ell=1,\ldots,L$, the linear operators $\boldsymbol{W}_{s}$, $\boldsymbol{W}_{\lambda}$, and their adjoint operators $\boldsymbol{W}_{s}^{\dagger}$, and $\boldsymbol{W}_{\lambda}^{\dagger}$. For example, the spatial and spectral operators are applied on images and spectra respectively in separate steps, see for example $\boldsymbol{W}_{s}\boldsymbol{X}$ in eq. (11) and $\boldsymbol{W}_{\lambda}\boldsymbol{X}^{\top}$ in eq. (12). Finally, for a better understanding of the derivation of algorithm 1, the reader can refer to Algorithm 3 (section 4) in (Condat, 2014).

From (11), the size of $\boldsymbol{U}$ is the size of $\boldsymbol{W}_{s}\boldsymbol{X}$ and consequently several times the size of the cube of images $\boldsymbol{X}$ when $\boldsymbol{W}_{s}$ is a redundant transform such as a union of orthogonal wavelets or an IUWT decomposition. Symmetrically to (11), the second version of the Vũ-Condat algorithm (Condat, 2014, Algorithm 3.2) requires, for the update of $\boldsymbol{X}$, the update of $\boldsymbol{U}$ in the current iteration and in the previous iteration. Keeping these two variables in memory can be too cumbersome when $\boldsymbol{X}$ is large and when the spatial decomposition is redundant. For this reason, we privileged the first version of the algorithm (Condat, 2014, Algorithm 3.1). Similarly, in order to reduce memory requirements, the relaxed version of the algorithm is not used ($\rho_{n}=1$).

2.3 Parallel implementation

The parallel implementation consists in parallelizing the steps that are carried on the images at different wavelength bands separately. The advantage of this parallelization is twofold. First the memory load can be distributed among several slave nodes, a slave node only stores the images at certain wavelength bands and not at all wavelength bands. Second, most of the computations and most interestingly the most expensive ones (such as the application of $\boldsymbol{W}_{s}$ and $\boldsymbol{W}_{s}^{\dagger}$) can be also carried out in parallel. More precisely, a group of $L$ slave nodes, where each slave node is assigned a band, handles the gradient computation, the application of the spatial operator $\boldsymbol{W}_{s}$ and its adjoint operator $\boldsymbol{W}_{s}^{\dagger}$.

The parallel implementation of algorithm 1 is described in algorithm 2. The computations in equations (9)-(11), and (13), except for the last term in equation (10) are computed in a parallel manner, for each band. The parallel computations of equations (9)-(11), and (13) in algorithm (1) are described in equations (14)-(20) in algorithm 2. After the slave nodes perform there computations in parallel, each node sends the new value of $\boldsymbol{x}_{\ell}$ and $\tilde{\boldsymbol{x}}_{\ell}$ to the master node. The master node centralizes the computation of the last term in equation (10) and equation (12), since these computations operate on spectra rather than images. It is worthy to note that if $\forall n$ $\gamma_{n}=0$, i.e. the spectral regularization is omitted, computation steps on the master (21)-(23) can be removed: each slave reconstructs the image corresponding to its wavelength independently to the others. From a memory point of view, the master node stores the whole dirty image cube and psf, whereas each slave node as mentioned previously only stores the dirty image and PSF at the corresponding band. Finally, without any loss of generality, a node can be assigned a group of bands, i.e. it can handle a sub-cube of the image rather than just one band. In this case, the communication overhead between the master node and the slave nodes can be reduced compared to the case where each band is assigned to a node.

3 Risk estimation for optimal parameter selection

The regularization parameters values $\mu$ and $\gamma$ defined in (4), control the strength of the spatial and spectral regularizations w.r.t. the strength of the data fidelity term. Setting the regularization parameters values appropriately is crucial to ensure a good performance of the deconvolution algorithm, and to obtain a good estimated image $\boldsymbol{X}$ of the true sky $\boldsymbol{X}^{\star}$. However, judging the estimated image quality and its closeness to the true and sought after image is a nontrivial and complicated task. This is due to the lack of ground truth, i.e. due to the fact that the true image is unknown. In addition to this, a ‘trial and error’ technique consisting in the inspection of multiple reconstructions is extremely time consuming when dealing with large-scale problems and multiple parameters.

The literature on optimal parameter selection, within the frequentist statistical paradigm, provides various quantitative measures for evaluating the quality of an estimated image without prior knowledge of the true image. These measures can be broadly classified as those based on the discrepancy principle (Karl, 2005; Morozov, 1966), the L-curve method (Regińska, 1996; Hansen & O’Leary, 1993), the generalized cross-validation (GCV) (Golub et al., 1979), and the (generalized) Stein Unbiased Risk Estimator (SURE) (Stein, 1981; Eldar, 2009). Among these measures, the SURE is one of the most attractive: unlike the GCV, the L-curve, and the discrepancy principle, the SURE benefits from optimal properties for nonlinear estimators and it is not restricted for linear ones. It provides an unbiased estimate of the risk i.e. the Mean Square Error (MSE), a quantitative measure often used for quality assessment in image processing. The SURE measure was originally applicable in the case of image denoising, i.e. when all the $\boldsymbol{H}_{\ell}$ are equal to the identity matrix in equation (1). Nevertheless, the concepts underlying the SURE have been generalized to the case of image reconstruction. Weighted versions of the SURE allow the unbiased estimation of a weighted MSE in the case where the $\boldsymbol{H}_{\ell}$ are not equal to the identity matrix as it is the case in equation (1). In this work, we consider a weighted version of the SURE known as the Predicted SURE (PSURE). The PSURE and other weighted SURE metrics have been successfully adopted for automatic tuning in image deconvolution, see for example (Ramani et al., 2012; Giryes et al., 2011; Deledalle et al., 2014; Ammanouil et al., 2017) and references therein.

A straightforward strategy for optimal parameter selection based on the SURE criteria or any other metric is by exhaustive search. More precisely, image reconstructions using all combinations of the regularization parameters are tested to achieve the minimum SURE estimate, see for example the work of Eldar (2009); Ramani et al. (2008). Compared to the grid search, a golden section allows to narrow the range of values inside which the extremum exists (Ramani et al., 2012). However, the golden section search is most suited for the case of one parameter value. In the case of multiple regularization parameters it can be either used to successively find the optimal regularization parameters values as in (Ammanouil et al., 2017) or more sophisticated method such as the simplex method can be used. In general, exhaustive search, or bisection strategies for multiple parameters can become rapidly computationally prohibitive.

More recently, Deledalle et al. (2014) proposed to minimize the SURE (or PSURE) using a gradient descent approach. This approach is based on an approximation of the SURE that is weakly differentiable w.r.t. the regularization parameters. This approximation, denoted as SURE-FDMC, relies on finite differences and Monte Carlo approximations (FDMC). The weak gradient of the SURE-FDMC, called the Stein Unbiased GrAdient estimator of the Risk (SUGAR) is used in a gradient descent scheme to find the optimal parameters. Among orther advantages, the gradient descent scheme allows setting multiple parameters more easily than grid search and golden section approaches.

The exhaustive search using SURE and the gradient descent using SUGAR all require running the reconstruction algorithm until convergence various times, with different values for the regularization parameters, which can be prohibitive when the algorithm is computationally complex. Giryes et al. (2011) proposed to update the regularization parameter in parallel to the reconstruction algorithm: at each iteration of the reconstruction algorithm, the regularization parameter is updated by a minimization of the SURE computed a number of iterations ahead using a golden section algorithm. As noticed by Giryes et al. (2011), for a given number of iterations of the reconstruction algorithm and golden section, the computational cost of this greedy approach is larger than the global approach. In this work, we propose to estimate the optimal regularization parameters in a more computationally efficient way. Towards this goal, they are updated in parallel to the reconstruction algorithm using a gradient descent strategy based on SUGAR. More precisely, the regularization parameters are updated every $\Delta$ iterations of the reconstruction algorithm such that the SURE-FDMC is reduced. We refer to this method as ADAptive Parameter Tuning (ADA-PT). To the best of our knowledge, (Garsden et al., 2015; Jiang et al., 2017) are the only works in image reconstruction to propose an ‘online’ self-tuning of the parameters. Their approach relies on the fact that in many sparse reconstruction algorithms based on primal optimization methods, the regularization parameters act at each iteration as thresholds levels on reconstructed images. This is the case in (Garsden et al., 2015) where Beck & Teboulle (2009) FISTA is used with a sparse synthesis prior. The regularization parameters can then be set as a function of a noise level estimation associated to the image update. In contrast with this approach, the application of the ADA-PT approach is not limited to MUFFIN and can be extended to most iterative reconstruction algorithms.

In what follows, we briefly review the essential concepts behind the estimation of the Predicted Stein Unbiased Risk Estimator based on Finite Differences and Monte Carlo simulations (PSURE-FDMC) and the Stein Unbiased GrAdient estimator of the Risk introduced by Deledalle et al. (2014). We derive the iterative steps necessary for the computation of the SUGAR in the case of MUFFIN. We then show how SUGAR is used to automatically and adaptively select the optimal regularization parameters int ADA-PT. We also show that the steps required for self-tuning can be performed in parallel and in a memory efficient way.

3.1 Risk estimation: PSURE-FDMC

The SURE provides an unbiased estimate of the MSE. However, for ill-posed inverse problems, it is generally not possible to estimate the MSE since the measurements may only contain partial information about the true image. The concepts underlying the SURE have been extended to the case of image reconstruction using weighted versions of the SURE measure (Ramani et al., 2012; Giryes et al., 2011): the PSURE. The predicted risk, or Predicted Mean Square Error (PMSE) for model (1) refers to the norm of the discrepancy between the true sky and its estimate in the measurement domain, in other terms it refers to the reconstruction error measured in the $\ell_{2}$-norm sense:

[TABLE]

where

[TABLE]

The PSURE is an unbiased estimator of the PMSE, but unfortunately it is not generally differentiable, e.g. when considering non-smooth regularizers as in (3). For this reason, we will focus from now on on the PSURE-FDMC. The PSURE-FDMC is an asymptotically unbiased estimator for the PMSE in the case of a linear model with additive white Gaussian noise as in (1). Let $\hat{\mathsf{R}}$ be the overall PSURE-FDMC. According to (Deledalle et al., 2014), the PSURE-FDMC is given by:

[TABLE]

where

[TABLE]

where $\mn@boldsymbol{\delta}_{\ell}$ are independent random Gaussian vector direction such that $\mn@boldsymbol{\delta}_{\ell}\sim\mathcal{N}(0,\boldsymbol{I})$, $\boldsymbol{d}_{\ell}=\boldsymbol{H}_{\ell}\boldsymbol{x}_{\ell}$ and $\overline{\boldsymbol{d}}_{\ell}=\boldsymbol{H}_{\ell}\overline{\boldsymbol{x}}_{\ell}$ with $\overline{\boldsymbol{x}}_{\ell}$ the output of the deconvolution algorithm when $\overline{\boldsymbol{y}}_{\ell}:=\boldsymbol{y}_{\ell}+\epsilon\mn@boldsymbol{\delta}_{\ell}$ is used instead of $\boldsymbol{y}_{\ell}$ for $\ell=1,\ldots,L$ i.e. when the observations are slightly perturbed, $\epsilon>0$ is a small constant that should be chosen as small as possible as long as it does not rise to numerical instabilities. Note that $\mn@boldsymbol{\delta}_{\ell}$ does not have to be necessary Gaussian as explained by Avron & Toledo (2011). As mentioned in the beginning of this section, the PSURE-FDMC provides an unbiased estimate of the PMSE:

[TABLE]

The PSURE-FDMC in equation (27) only depends on the measurements, the outputs of the deconvolution algorithm ($\boldsymbol{x}_{\ell}$ and $\overline{\boldsymbol{x}}_{\ell}$, $l\in 1,\ldots,L$), and requires knowing the noise variances $\sigma_{w}^{2}$. Knowledge of the noise variance is generally assumed in interferometric imaging, as in (Onose et al., 2016b; Carrillo et al., 2014), or it can be estimated in the image domain using as in (Garsden et al., 2015) robust statistics such as the median absolute deviation (MAD). Most importantly, note how the PSURE-FDMC in equation (27) does not require the knowledge of $\boldsymbol{x}_{\ell}^{\star}$, unlike its PMSE counterpart in (25). As a result, it measures the quality of the estimation for given values of the regularization parameters, and the regularization parameters that gives the best estimate (i.e. the lowest PSURE-FDMC) are the best one among the tested values.

3.2 Gradient estimation: SUGAR-FDMC

Contrarily to the SURE, the PSURE-FDMC defined in (27) is differentiable in the weak sense w.r.t. the regularization parameters. Noting that $\boldsymbol{x}_{\ell}$ and $\overline{\boldsymbol{x}}_{\ell}$ are the only terms that depend on the regularization parameters, a straightforward derivation of PSURE-FDMC in equation (27) w.r.t. one of the regularization parameters yields the SUGAR-FDMC:

[TABLE]

where $\boldsymbol{d}_{\ell}^{\prime}=\boldsymbol{H}_{\ell}\boldsymbol{x}_{\ell}^{\prime}$ and $\overline{\boldsymbol{d}}_{\ell}^{\prime}=\boldsymbol{H}_{\ell}\overline{\boldsymbol{x}}_{\ell}^{\prime}$, and the notation $(\cdot)^{\prime}$ refers to the weak derivative w.r.t. one of the regularization parameters, in other terms $(\cdot)^{\prime}$ corresponds to either $\frac{\partial(\cdot)}{\partial\mu}$ or $\frac{\partial(\cdot)}{\partial\gamma}$. Interestingly, it was shown by Deledalle et al. (2014) that SUGAR-FDMC defined in (29) as the weak gradient of PSURE-FDMC, is an unbiased estimate of the weak gradient of the PMSE in (25) when $\epsilon$ tends to zero:

[TABLE]

which means that equation (29) gives access to an unbiased estimate of the risk’s gradient without prior knowledge of the true sky.

The computation of SUGAR-FDMC requires $\boldsymbol{x}_{\ell}^{\prime}$ and $\overline{\boldsymbol{x}}_{\ell}^{\prime}$, the weak derivatives of $\boldsymbol{x}_{\ell}$ and $\overline{\boldsymbol{x}}_{\ell}$ respectively, for $\ell=1,\ldots,L$. Given that $\boldsymbol{x}_{\ell}$ and $\overline{\boldsymbol{x}}_{\ell}$ are estimated iteratively, $\boldsymbol{x}_{\ell}^{\prime}$ and $\overline{\boldsymbol{x}}_{\ell}^{\prime}$ are also estimated iteratively. This is done by simply deriving each equation in algorithm 2 w.r.t. the regularization parameter under investigation. The iterative algorithms necessary for estimating $\frac{\partial\boldsymbol{x}_{\ell}}{\partial\gamma}$ and $\frac{\partial\boldsymbol{x}_{\ell}}{\partial\mu}$ are described in algorithms 3 and 4 respectively. The weak derivatives $\frac{\partial\overline{\boldsymbol{x}}_{\ell}}{\partial\gamma}$ and $\frac{\partial\overline{\boldsymbol{x}}_{\ell}}{\partial\mu}$ are estimated by replacing $\boldsymbol{y}_{\ell}$ by $\overline{\boldsymbol{y}}_{\ell}$ in algorithms 3 and 4. Each iteration in algorithm 3 is obtained by deriving the corresponding iteration in algorithm 2 w.r.t. $\mu$. When a variable implicitly depends on $\mu$, which is the case for $\boldsymbol{x}_{\ell}$ and all other intermediate variables, this consists in simply replacing the variable by its derivative, see for example (36). When the variable depends on the product between $\mu$ and another variable (that depends on $\mu$), this requires a simple application of the chain rule, see for example (38). The functions $\mathcal{U}(\cdot)$ and $\Pi(\cdot)$, in equations (40) and (42) respectively, correspond to the weak gradients of the projection on the positive orthant function and the saturation function defined in equations (6) and (7) respectively. They are defined as follows:

[TABLE]

and

[TABLE]

When the input is a matrix, they are applied component-wise. The same concepts are applied for the derivation of algorithm 4.

Running $m$ iterations of algorithms 3 and 4 with algorithm 2 starting from iteration $k$, provides an iterative procedure for estimating the weak derivative of the reconstructed image $m$ steps ahead $k$ w.r.t. $\mu$ and $\gamma$ respectively, and similarly for $\overline{\boldsymbol{x}}_{\ell}$ when $\boldsymbol{y}_{\ell}$ is replaced by $\overline{\boldsymbol{y}}_{\ell}$. Plugging these derivatives in equation (29) provides an asymptotically unbiased estimate of the gradient of the risk $m$ iteration ahead $k$. We will denote the corresponding gradient as $\nabla\hat{\mathsf{R}}_{k:k+m}(\boldsymbol{\theta})$ where $\boldsymbol{\theta}$ represents either $\gamma$, $\mu$ or $(\gamma,\mu)^{\dagger}$.

3.3 Adaptive Parameter selection

Ideally one would like to find the optimal parameters $\mu^{\star}$ and $\gamma^{\star}$ such that the PMSE at the convergence of the MUFFIN algorithm is minimized. An alternative solution for estimating the regularization parameters is to use the PSURE-FDMC instead of the PMSE:

[TABLE]

In order to solve problem (33) avoiding grid search techniques and for multiple patrameters, the authors of (Deledalle et al., 2014) propose to take advantage of the gradient estimate of the PSURE-FDMC defined in the previous section in order to find the optimal parameters using a gradient descent scheme at convergence of the reconstruction algorithm. The gradient descent approach finds the optimal regularization parameters according to the following updates:

[TABLE]

where $\eta$ is a sufficiently small step size, $\nabla\hat{\mathsf{R}}_{0:n_{c}}$ is computed as explained in the previous section and $n_{c}$ corresponds to the number of iterations required to achieve convergence of MUFFIN. Note that $n_{c}$ may depend on $\boldsymbol{\theta}^{(i)}$ but can be set without any additional cost monitoring the successive values of $k\mapsto\|\nabla\hat{\mathsf{R}}_{0:k}(\boldsymbol{\theta}^{(i)})\|$. This gradient descent approach requires at iteration $i$ running until convergence MUFFIN with the regularization parameter value $\boldsymbol{\theta}^{(i)}$ and the derivatives $\nabla\hat{\mathsf{R}}$. Afterwards, the parameters $\boldsymbol{\theta}^{(i)}$ are updated to $\boldsymbol{\theta}^{(i+1)}$, MUFFIN and the derivatives estimation algorithms should be run again and so on.

In order to reduce the computational cost of this algorithm, we propose an adaptive gradient descent approach that simultaneously updates the estimates of $\boldsymbol{x}^{\star}_{\ell}$ and the regularization parameters. This approach, in the spirit of (Giryes et al., 2011), was introduced by Ammanouil et al. (2018) for general iterative algorithms. More precisely, the regularization parameters are updated every $\Delta$ iterations of MUFFIN according to a gradient descent direction such that the PSURE-FDMC at the current iteration is reduced, see Fig. 1. At iteration $k$ of MUFFIN, if $k=(q+1)\Delta$, $\mn@boldsymbol{\theta}^{(q)}$ is updated using

[TABLE]

Hence, rather than running the MUFFIN algorithm 2 until convergence and then updating the regularization parameters using (34), it requires running algorithms 2, 3 and 4 only once. We refer to this method as ADAptive Parameter Tuning (ADA-PT). More precisely, the ADA-PT strategy consists in the iterations described in Algorithm 5. The parameter $\Delta$ sets the trade-off between a greedy algorithm, $\Delta=1$, which will update the parameters at each iteration and the global algorithm (34), $\Delta=n_{c}$. As pointed by Giryes et al. (2011), one of the problems of the greedy strategy ($\Delta=1$) is that it will update the parameters in order to reduce the estimated PMSE of the current iteration but can harm the overall PMSE. Thus, instead of choosing to reduce the estimated PMSE at each iteration, ADA-PT reduces the estimated PMSE every $\Delta>1$ iterations, typically few tens. It is worthy to note that the value of $\Delta$ does not affect the complexity of the algorithm.

Finally, it is worth pointing out that although estimating the risk gradient entails the additional operations described in algorithms 3 and 4, the global memory requirement and computational complexity of ADA-PT’ed MUFFIN are multiplied only by three w.r.t. MUFFIN. Moreover, these additional computations benefit from the same parallel implementation as MUFFIN, described in section 2.3.

4 Tests on simulated data

The reconstruction results reported in this section were obtained with the proposed distributed implementation of MUFFIN. The aim of these simulations is twofold: 1) to demonstrate the efficiency of the distributed implementation; 2) to illustrate the relevance of the ADA-PT strategy for automatically finding the optimal regularization parameters.

4.1 Data generation description

The data sets used in the following are composed by two simulated, multi-wavelength radio sky emissions in which we have introduced complex spatial and spectral features in order to challenge the reconstruction algorithms. The raw data are obtained by convolving these sky emissions with wavelength dependent PSFs and adding noise. We describe in more details the considered sky emissions and the PSFs below. Note that in all Figures the flux scale is arbitrary as we are mostly interested in comparing the relative intensities of reference and reconstructed images/cubes. Point estimation of the Signal to Noise Ratio (SNR) and of the PMSE denoted as WMSE:

[TABLE]

will however be indicated.

4.1.1 Simulated galaxy with HI emission line (‘Galaxy-HI’)

The first data set is a RA-Dec-$\nu$ cube generated from a well known image representing a star-forming region in the M$31$ galaxy (see Fig. 2, top left panel). For each of the $256\times 256$ pixels in this image, we computed a spectrum in $256$ adjacent frequency bands of constant bandwidth ($2\,\text{MHz}$), starting from $950\,\text{MHz}$. Each spectrum follows a first order power-law spectrum model depending on a spectral index. The $256\times 256$ map of spectral indices was constructed following the procedure detailed by Junklewitz et al. (2015): for each pixel, the spectral index is taken as a linear combination of a homogeneous Gaussian field and the reference sky image.

As an additional spectral feature, we inserted in the spectra of the brightest pixels (see legend of Fig. 2) a HI line simulating the localized presence of neutral hydrogen gas. The line is centered at $1.42$ GHz, corresponding to our $235^{\textrm{th}}$ channel. We assume here a simple model where the galaxy has a gas consumption time scale of typically 1 Gyr, a H mass of $10^{9}$ solar masses and thus a star formation rate of one solar mass per year. Assuming further a velocity dispersion of $200$ km/s, we obtain lines whose flux is about twice the flux of the continuum at $1.42$ GHz, a value that may typically fluctuate within one order of magnitude (P. Serra, private communication). Translated into frequency, the velocity dispersion leads to a line width of 1 MHz, with tails that slightly smear out on adjacent channels. In view of these elements, the spectral HI line profile was computed as a continuous Gaussian emission line profile with a maximum corresponding to twice the local continuum and a FWHM of 1 MHz. The profile was further convolved by a boxcar function of width $2\,\text{MHz}$, sampled at the locations of the considered $256$ frequency bands and added to the continuum. The top right panel of Figure 2 shows the simulated galaxy in the channel where the HI line has most energy ($235^{\text{th}}$ wavelength band). Figure 2 bottom panel shows the spectrum of a bright pixel from the Galaxy data set where the HI emission was added. We refer to this simulated sky data cube as ‘Galaxy-HI’ for short.

4.1.2 Simulated galaxy cluster with diffuse halo (‘Halo’)

The second data set is a galaxy cluster radio model based on a refined version of 2D radio simulations considered in our previous works (Ferrari et al., 2015). This model includes various kinds of radio sources that can be found in ‘radio-loud’ clusters, in particular unresolved point-like objects, bright elongated radio galaxies with complex morphologies and a giant low-surface brightness source of synchrotron radiation (or ‘radio halo’), see e.g. (Ferrari et al., 2008; Feretti et al., 2012) and references therein. In addition to these features, this simulation also properly models the number counts of both the cluster radio source and the foregound/background population. The simulated radio halo has been produced following the approach described in Murgia et al. (2004) while discrete radio sources have been modelled according to what described in Loi et al. (2018)

To generate the 3D data cube, we segmented the simulated 2D radio map described above into three contributions: two giant radio galaxies (with the typical head-tail and narrow-angle-tail morphology of FRI in clusters) and a low-surface brightness extended halo. Each of the three contribution was given a specific power law spectrum $\boldsymbol{e}_{k}$ ($k=1,2$ or $3$) that remains spatially constant up to a multiplicative factor $\alpha_{n}(k)$. Formally, spectrum $n$ of the blended data cube is thus generated as:

[TABLE]

The three considered spectra $\{\boldsymbol{e}_{k}\}$, also called endmembers, are shown in the left panel of Figure 3. The three considered segmentation maps (also called abundance maps), which show in each pixel the fraction $\alpha_{n}(k)$ of the $k^{\textrm{th}}$ contribution, are displayed in the three other panels of Figure 3. Note that the halo is normalized such as its power equals approximately 11 % of the power of the final cube.

The purpose of this data set is twofold. First, to evaluate the ability of the proposed method to reconstruct scenes with very high dynamic range and complex structures (the halo in particular), which is one of SKA first objectives. Second, to show that, thanks to MUFFIN’s non parametric approach, the reconstructed cube can be further processed to achieve a specific task such as the deblending of the sources using their spectral properties. For this spectral unmixing case, a two-step procedure (reconstruction and calibration within a major/minor loop followed by the deblending), compared to a joint procedure such as in (Jiang et al., 2017; Abdulaziz et al., 2016), can be advocated since i. it avoids the resolution of global non convex optimization problems and ii. it gives some space to a human expertise during the unmixing task, e.g. for identifying the endmembers when complex objects are superposed to a halo.

The final cube, which will be referred to as ‘Halo’ for short, spans a frequency range from 950 MHz to 1.458 GHz, with 256 spectral channels of equal bandwidth.

4.1.3 Point spread function (PSF)

We simulated the PSFs of SKA1-MID with its 197 dishes using the HI-inator package111https://github.com/SpheMakh/HI-Inator based on the MeqTrees software by Noordam & Smirnov (2010). The simulated PSFs correspond to a Fourier coverage produced by a total observation time of $8$ hours. Figure 4 shows the PSF at the first wavelength band. The PSF in each channel is characterized by a central lobe and ringed side-lobes at larger angular distances. This cube of PSFs was used to convolve channel by channel the Galaxy-HI and the Halo data sets. Gaussian noise was added to the result to produce the dirty data cubes with an SNR of 20 dB for the dirty Galaxy-HI data set and an SNR of 40 dB for the dirty Halo data set.

4.2 Deconvolution examples

4.2.1 Implementation details

We tested the proposed deconvolution approach, namely ADA-PT’ed MUFFIN, on the two data sets described above. In all the simulations, $\textbf{W}_{s}$ and $\textbf{W}_{\lambda}$ were respectively set to the union of the first $8$ Daubechies wavelet bases and to the Discrete Cosine Transform (DCT). The parameter of the deconvolution algorithm $\sigma$ was set to $\sigma=10$ and the corresponding value of $\tau$ was calculated as explained in section 2.2. The parameter $\Delta$ was set to $10$.

We compare the results of the proposed approach with the MS-MFS (Multi-Scale Multi-Frequency) deconvolution algorithm. We used the implementation developed within the SKA Science Data Processor Consortium available online222https://github.com/SKA-ScienceDataProcessor/algorithm-reference-library with the default parameters. In particular, the number of iterations is equal to $10^{3}$, the scales are set to $(0,3,10,30)$ and order 3 polynomials are used to fit the spectra. We tested MS-MFS using the default parameters and using other parameters, the default parameters gave the best results.

This work was granted access to HPC ressources. In particular, the simulations were performed using $80$ compute nodes, 1 master node and 79 slave nodes, with $64$ Go of RAM per node. This setting corresponds to $4$ spectral bands at most per slave node. The algorithm is written in Python, and the Message Passing Interface (MPI) is used to deal with the parallel computing architecture. The computing time to process a $256\times 256\times 256$ cube is: 3 sec./iteration for MUFFIN and 15 sec./iteration for the full ADP-PT’ed code. The code is available at https://github.com/andferrari/MUFFIN.

The most computationally intensive steps of MUFFIN are the $L=256$ spatial wavelets decompositions and recompositions that are performed at each iteration. Nevertheless, these steps are performed in parallel, see Eq. (16, 19) in Alg. 2. On the other hand, MS-MFS requires 12 scale decompositions (4 scales and order 3 polynomials) that are performed at the beginning of the algorithm, i.e. before the first iteration. Memory wise, the amount of memory required by MUFFIN is to store the matrices $\boldsymbol{U}$, $\boldsymbol{X}$, $\tilde{\boldsymbol{X}}$ and $\boldsymbol{V}$. As a result, the global memory footprint of MUFFIN mainly is proportional to the total size of these matrices which is equal to $11LN$. Note that the memory load is spread across the $L$ computing nodes. On the other hand, the MS-MFS memory usage is $12N$. The memory requirements for MUFFIN are doubled for ADA-PT’ed MUFFIN.

4.2.2 Galaxy-HI data set

The spatial and spectral regularization parameters $(\mu,\gamma)$ were initialized to four different couples of values: $(1,10)$, $(1,3)$, $(0,005,3)$ and $(0.005,10)$. The ADA-PT’ed MUFFIN algorithm was ran for $10^{4}$ iterations and Figure 5 shows the evolution of each regularization parameters for the four different initialization settings as a function of the iterations. In all cases, the regularization parameter $\mu$ converges to $0.05$ and the regularization parameter $\gamma$ converges to $7$ (on average). Note that the four initialization settings cover the fours distinct scenari where both parameters values are greater than the optimal ones, less than the optimal ones, and where one is greater and one is less than the optimal value.

Figure 6 shows the evolution of $k\rightarrow\nabla\hat{\mathsf{R}}_{q\Delta:k}(\boldsymbol{\theta}^{(k)})$, i.e. of SUGAR w.r.t. $\mu$ and of SUGAR w.r.t. $\gamma$ as a function of the iteration number for the four considered initializations values. Note how both SUGAR curves converge to zero, which justifies the regularization parameters convergence. Figure 7 shows the variation of SNR and WMSE as a function of the iteration number for the four considered initialization cases. The SNR converges to a value of 24 dB and the WMSE converges to $3.4\times 10^{-4}$. Note that the SNR and WMSE curves are only shown for the first $5000$ iterations as they remain stable after the first $2000$ iterations.

We compare in Figure 9 the proposed approach with Rau & Cornwell (2011) MS-MFS. The third and fourth column in the first row of Figure 9 show that the deconvolved image obtained using MUFFIN and MS-MFS are both very close to the sky image and ‘clean’ when compared to the dirty image shown in the second column. The second row of Figure 9, which highlights the low intensities, reveals that MUFFIN was able to better reconstruct the image while some smooth artifacts are seen in the MS-MFS image. Note that the SNR of the reconstructed image cube is $24.6$ dB for MUFFIN and $18.8$ for MS-MFS. Finally, Figure 8 shows the reconstruction results for the spectrum of the bright pixel in the image characterized by an HI emission and shown in Fig. 2. It can be seen that both MUFFIN and MS-MFS are able to estimate and denoise the continuum part of the spectrum. Turning to the inset panel (zoom on the HI emission line), we see that the line is logically absent in the MS-MFS spectrum, because this algorithm is aimed at reconstructing a power law spectrum. In contrast, MUFFIN is able to preserve the HI emission in the reconstructed spectra.

4.2.3 Halo data set

The proposed ADA-PT strategy was run on the Halo data set, for different initialization settings and for $10^{4}$ iterations. The regularization parameter $\mu$ and $\gamma$ converge respectively to $7\times 10^{-6}$ and to $2.6\times 10^{-2}$. The SNR converges to a value of 31 dB, and the WMSE converges to $1.6\times 10^{-7}$.

Figure 10 compares the performance of the proposed approach with MS-MFS. The third (MUFFIN) and fourth (MS-MFS) columns lead to the same conclusion as for the Galaxy-HI data set. Note that the SNR of the reconstructed image cube is $31.44$ dB for MUFFIN and $14.06$ for MS-MFS.

Finally, Figure 11 shows the source separation result. We chose three spectra belonging to brightest pixels: one from the faint halo region, one from the largest galaxy and one from the smallest galaxy. The spectra of these pixels are shown in the first column of Figure 11. We then solved a positively constrained least squares problem in order to find the contribution of each spectrum for each pixel of the reconstructed data cubes. Columns $2$ to $4$ show the recovered abundance maps (to be compared to the maps shown in Fig.3): the first row corresponds to MUFFIN and the second row to MS-MFS. It can be seen that in the case of MUFFIN, the two galaxies are better separated than in the MS-MFS case. As for the Halo, it is well separated in both cases, but for MS-MFS it presents larger artifacts and differences with respect to the reference.

5 Conclusions

In this work, we have developed a spatial-spectral deconvolution method for large multispectral data named MUFFIN. For this algorithm, we designed a specific parallel implementation and an original procedure for automatically self-tuning the regularization parameters. We demonstrated the efficiency of the proposed deconvolution method on simulated but realistic astronomical multispectral images. The parallel strategy and the SURE-based automatic parameter selection were specifically adapted to the spatial-spectral reconstruction strategy of the MUFFIN algorithm. However, the tools used in this work should provide a foundation for developing analogous parallel implementations and automated parameter selection procedures in the case of different deconvolution problems.

Acknowledgements

This work was supported by the Agence Nationale pour la Recherche, France, (MAGELLAN project, ANR-14-CE23-0004-01). The numerical simulations in this work were performed at the Centre for Intensive Computing ‘Mésocentre sigamm’ hosted by the Observatoire de la Côte dAzur. The authors would like to thank Alain Miniussi for his help on the cluster and P. Serra from Osservatorio Astronomico di Cagliari for his help on the model of the HI line. The authors are very grateful to the referee, Stephen Ord, for his insightful suggestions and helpful comments to improve the quality of the original manuscript.