CLEANing Cygnus A Deep and Fast with R2D2

The R2D2 algorithm (Aghabiglou et al. 2024) is underpinned by a sequence of DNN modules denoted by $\{{{\bf\mathsf{N}}}_{{\widehat{{\boldsymbol{\theta }}}}^{(i)}}\}{}_{1\leqslant i\leqslant I}$, which are described by their learned parameters $\{{\widehat{{\boldsymbol{\theta }}}}^{(i)}\in {{\mathbb{R}}}^{Q}\}{}_{1\leqslant i\leqslant I}$. With the image estimate being initialized as ${{\boldsymbol{x}}}^{(0)}={\bf{0}}\in {{\mathbb{R}}}^{N}$, at any iteration i ∈ {1,...,I}, the network ${{\bf\mathsf{N}}}_{{\widehat{{\boldsymbol{\theta }}}}^{(i)}}$ takes as input both the previous image estimate x ⁽ⁱ⁻¹⁾ and associated residual dirty image r ⁽ⁱ⁻¹⁾. The latter is updated from the dirty image by removing the contribution of the previous image estimate: r ⁽ⁱ⁻¹⁾= x _d − D x ⁽ⁱ⁻¹⁾. The current image estimate is updated from the output of the network as

$$\begin{eqnarray}&&{{\boldsymbol{x}}}^{(i)}={{\boldsymbol{x}}}^{(i-1)}+{{\bf\mathsf{N}}}_{{\widehat{{\boldsymbol{\theta }}}}^{(i)}}({{\boldsymbol{r}}}^{(i-1)},{{\boldsymbol{x}}}^{(i-1)}),\end{eqnarray} \tag{ 2 }$$

where the predicted residual image ${{\bf\mathsf{N}}}_{{\widehat{{\boldsymbol{\theta }}}}^{(i)}}({{\boldsymbol{r}}}^{(i-1)},{{\boldsymbol{x}}}^{(i-1)})$ captures emission from the input residual dirty image r ⁽ⁱ⁻¹⁾, and corrects for estimation errors in x ⁽ⁱ⁻¹⁾, thus progressively enhancing the resolution and dynamic range of the reconstruction. Assuming I DNN modules in the R2D2 sequence, the final reconstruction $\widehat{{\boldsymbol{x}}}$ takes a simple series expression, as the sum of output residual images from all DNN modules. The iteration structure of R2D2 is reminiscent of MP, but its model components at each iteration are identified in a learned manifold of basis functions encoded by a DNN, rather than a handcrafted sparsity basis fixed across iterations.

R2D2 DNN modules are trained in a sequential manner from a data set consisting of K image pairs of the ground-truth images and their associated dirty images $\{{{\boldsymbol{x}}}_{k}^{\star },{{{\boldsymbol{x}}}_{{\rm{d}}}}_{k}\}{}_{1\leqslant k\leqslant K}$. At any iteration i ≥ 1, the network ${{\bf\mathsf{N}}}_{{\widehat{{\boldsymbol{\theta }}}}^{(i)}}$ is trained taking the previous image estimates and associated residual dirty images $\{{{\boldsymbol{x}}}_{k}^{(i-1)},{{\boldsymbol{r}}}_{k}^{(i-1)}\}{}_{1\leqslant k\leqslant K}$ as input, with $\{{{\boldsymbol{x}}}_{k}^{(0)}={\bf{0}}\}{}_{1\leqslant k\leqslant K}$ and $\{{{\boldsymbol{r}}}_{k}^{(0)}={{{\boldsymbol{x}}}_{{\rm{d}}}}_{k}\}{}_{1\leqslant k\leqslant K}$. Its learned parameter vector ${\widehat{{\boldsymbol{\theta }}}}^{(i)}$ is minimizer of the loss function:

$$\begin{eqnarray}&&\mathop{\min }\limits_{{{\boldsymbol{\theta }}}^{(i)}\in {{\mathbb{R}}}^{Q}}\displaystyle \sum _{k=1}^{K}\,\Vert {{\boldsymbol{x}}}_{k}^{\star }-{[{{\boldsymbol{x}}}_{k}^{(i-1)}+{{\bf\mathsf{N}}}_{{{\boldsymbol{\theta }}}^{(i)}}({{\boldsymbol{r}}}_{k}^{(i-1)},{{\boldsymbol{x}}}_{k}^{(i-1)})]}_{+}{\Vert }_{1},\end{eqnarray} \tag{ 3 }$$

where ∥. ∥₁ stands for the ℓ₁-norm, and [.]₊ denotes the projection onto the positive orthant ${{\mathbb{R}}}_{+}^{N}$, ensuring nonnegativity of the image estimate at any point in the iterative sequence. The image estimates and associated residual dirty images $\{{{\boldsymbol{x}}}_{k}^{(i)},{{\boldsymbol{r}}}_{k}^{(i)}\}{}_{1\leqslant k\leqslant K}$ are updated from the outputs of the trained network. Alongside their corresponding ground-truth images, they serve as the training data set for the subsequent DNN module. Training of the DNN series concludes when the reconstruction metrics on the validation data set reach a point of saturation.

The R2D2 algorithm takes three variants. The first, simply referred to as R2D2 henceforth, adopts the well-established U-Net as the architecture of its DNN modules. The second, called R2D2-Net, is an unrolled version of R2D2 itself, trained end-to-end on GPUs. Its architecture consists of J U-Net layers interlaced with J − 1 data-consistency layers. The latter layers leverage a fast and memory-efficient PSF-based approximation of the operator D for the computation of residual images, which is instrumental to alleviate the scalability challenges encountered by unrolled deep-learning approaches. The third is a nested R2D2 architecture, where the U-Net modules of R2D2 are substituted with R2D2-Net. In reference to its nested structure, this variant is dubbed "Russian Dolls" R2D2, in short R3D3. In this landscape, the R2D2-Net variant, as DNN module to R3D3, also represents R3D3's first iteration.

The Cygnus A data considered were acquired at S band (2.052 GHz) with the following observation setting. VLA configurations A and C were combined, with a single pointing centered at the inner core of Cygnus A, given by the coordinates R. A. = 19h59m28.356s (J2000) and $\mathrm{decl}.\,=\,+{40}^{{\rm{o}}}4^{\prime} 2\buildrel{\prime\prime}\over{.} 07$. The total observation duration is about 19 hr conducted over multiple scans, respectively combining 7 and 12 hr with configurations A and C. The data are single-channel, acquired with integration time step of 2 s and channel width of 2 MHz. Careful calibration of the direction independent effects was performed in AIPS (Sebokolodi et al. 2020), following which, the data were averaged over 10 s totaling M = 9.6 × 10⁵ points. Dabbech et al. (2021) showed that these data could also benefit from the calibration of DDEs, likely attributed to pointing errors at the hotspots. However, DDE solutions from the underpinning joint calibration and imaging framework (Repetti et al. 2017) are not considered in this study as none of the benchmark algorithms aside from uSARA and AIRI could benefit from them.

As far as the imaging setting is considered, we target forming an image of size N = 512 × 512, with a pixel size of about 029, corresponding to a superresolution factor of 1.5, from Briggs-weighted data with Briggs parameter set to 0 (weights are generated in WSClean; Offringa & Smirnov 2017). The target dynamic range of the reconstruction is about 1.7 × 10⁵, and is inferred from the data as the ratio between the peak pixel value of the sought image and the image-domain noise level estimated as η_Briggs $\tau /\sqrt{2\parallel {{\boldsymbol{\Phi }}}^{\dagger }{\boldsymbol{\Phi }}{\parallel }_{S}}$ , with η_Briggs > 0 being a correction factor accounting for the Briggs weights (Terris et al. 2022; Wilber et al. 2023a), and ∥. ∥_S denoting the spectral norm of its argument operator.

We have borrowed the first incarnations of the R2D2 variants from Aghabiglou et al. (2024). These were trained in a telescope-specific approach, more precisely for VLA. R2D2 and R3D3 were respectively trained with I = 15 U-Net modules, and I = 8 R2D2-Net modules, all sharing the same architecture composed of J = 6 U-Net layers. The R2D2-Net incarnation considered is exactly the first DNN module of the R3D3 incarnation.

The training data set was composed of 20,000 pairs of ground-truth images and associated dirty images of size N = 512 × 512. Synthetic ground-truth images, endowed with high dynamic ranges in the interval [10³, 5 × 10⁵], were created from a low-dynamic-range database of optical astronomy and medical images via denoising and exponentiation procedures (Terris et al. 2022, 2023). Targeting a robust incarnation, usable across a wide landscape of VLA observation settings, data were created from a variety of Fourier sampling patterns generated using MeqTrees (Noordam & Smirnov 2010), combining configurations A and C. The sampling patterns were generated by sampling uniformly at random: (i) the pointing direction; (ii) the temporal parameters (total observation duration assuming a single scan and a constant time step of 36 s for a combined duration at both configurations ranging between 6 and 13 hr); (iii) the spectral parameters (number of consecutive frequency channels and frequency channel bandwidth) under the assumption of flat spectra of the radio emission. Additionally, random flagging of up to 20% of the Fourier samples was applied. The simulated visibilities of size M ∈ [0.2, 2] × 10⁶ were corrupted with noise levels commensurate with the dynamic ranges of the corresponding ground-truth images. In order to contain the number of varying parameters in the training, the imaging setting was fixed, with a N = 512 × 512 image size, and dirty images obtained: (i) using the same weighting scheme across all data, that is Briggs weighting whose parameter is set to 0; (ii) with a pixel size ensuring a superresolution factor of 1.5.

We note that, while the imaging setting for the pretrained R2D2 variants and the target image reconstruction of Cygnus A are the same, the VLA temporal parameters of the observation settings (total duration, integration time step, number of scans) belong to quite different categories. Our working hypothesis in using these pretrained variants is that the wide variety of observation settings visited at training stage endows the DNN modules with the necessary robustness, when interlaced with the update of the residual dirty image at each iteration, to generalize across different categories of observation settings. The same question in fact formally holds with regards to the target category of images, as the training data set contains no radio image whatsoever. The capability to train on synthetic data and reconstruct radio images was extensively validated in simulation (Aghabiglou et al. 2024), but remains to be demonstrated on real data.

Python and MATLAB implementations of the R2D2 variants that can be executed on flexible hardware are available as part of the BASPLib code library. Extensive validation of image reconstruction in simulation, as well as a detailed analysis of their computational performance at both training and imaging stages are provided in Aghabiglou et al. (2024). In this study, image reconstruction with R2D2 variants was conducted using their Python implementation, executed on a single GPU.

The performance of the three R2D2 variants is studied in comparison with the Hö-CLEAN, CS-CLEAN, MS-CLEAN, uSARA, and AIRI. Parameters of CLEAN variants were adjusted to ensure a good compromise between speed and data fidelity, while achieving a good reconstruction quality. Both uSARA and AIRI are shipped with automated noise-driven heuristics for the choice of their regularization parameter (Terris et al. 2022; Wilber et al. 2023a). However, an adjustment within 1 order of magnitude from their heuristic values appeared to be necessary for best results. For CLEAN variants, we have considered the widely used C++ software WSClean (Offringa & Smirnov 2017; see the Appendix for full commands). Both uSARA and AIRI are implemented in MATLAB as part of the BASPLib code library. All CLEAN variants and uSARA were executed on a single CPU core. AIRI was executed on a single CPU core for its data-fidelity steps and single GPU for its regularization steps (denoising DNNs).

Cygnus A images obtained by the different imaging methods are displayed in ${\mathrm{log}}_{10}$ scale in Figures 1–4. Reconstructions are overlaid with additional panels consisting of (a) the associated residual dirty images displayed in linear scale to visually assess the fidelity to back-projected data, and zooms on selected regions of the radio galaxy, all displayed in ${\mathrm{log}}_{10}$ scale: (b) the inner core, (c) the west hotspots, and (d) the east hotspots. The overall visual inspection of Cygnus A reconstructions shows that R2D2 variants exhibit higher resolution than CLEAN variants, while generally corroborating the achieved depictions by AIRI and uSARA. They provide deep reconstructions, whose pixel values span nearly 5 orders of magnitude, which is in line with the target dynamic range estimate. A close-up inspection indicates that both R3D3 (Figure 3, bottom) and AIRI (Figure 4, bottom) stand out, owing to their high levels of detail and their limited amount of patterns that could be construed as artifacts. R2D2 (Figure 2, bottom) and R2D2-Net (Figure 3, top) seem to lack details in the faint extended emission. uSARA (Figure 4, top) depicts spurious ringing and wavelet-like patterns. As expected, Hö-CLEAN (Figure 1, top) delivers a poor reconstruction with severely limited dynamic range, due to its inherent approximate data model. Both CS-CLEAN (Figure 1, bottom) and MS-CLEAN (Figure 2, top) provide much improved reconstructions, with the former exhibiting grid-like artifacts due to its inadequate sparsity (identity) basis for the complex target radio source.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The inner core. The inner core consists of the point-like active galactic nucleus (AGN) of Cygnus A, from which two jets emanate (panels (b) of all figures). The reconstructions of R2D2 variants, uSARA, and AIRI show a superresolved depiction of the region. In particular, R2D2 variants exhibit continuous emission between the AGN and both jets. uSARA exhibits wavelet-like artifacts around the AGN. CLEAN variants provide unresolved depiction of the source due to the restoring beam.

The lobes. Examination of the west and east lobes of Cygnus A highlights the ability of R2D2 variants to provide a more physical depiction of their filamentary structure than the benchmark algorithms. On the one hand, CLEAN variants deliver a smooth reconstruction. On the other hand, uSARA, and to a much lesser extent AIRI, exhibit ringing artifacts in the west lobe (pointed at with a green arrows in Figure 4). These artifacts are likely induced by pointing errors at the hotspots, resulting in the overfitting of the high-spatial-frequency content of the data by both uSARA and AIRI. Joint DDE calibration and imaging (using either AIRI or uSARA as the imaging module) can drastically reduce (if not remove) these artifacts (see their corresponding reconstructions provided in Dabbech et al. 2024). These findings suggest that R2D2 variants may be less prone to calibration errors than AIRI and uSARA.

Focusing on the faint diffuse emission at the tails of the lobes, both R3D3 and AIRI provide consistent depiction. When flipping between their associated figures (Figures 3–4, bottom), the R3D3 reconstruction appears sharper. One such example is the faint filamentary structure at the tail of the east lobe. However, one cannot ascertain whether this sharpness is physical, particularly in absence of DDE calibration. Both R2D2 (Figure 2, bottom) and R2D2-Net (Figure 3, top) provide a rather smooth structure. uSARA introduces wavelet-like patterns (pointed at with red arrows in Figure 4, top). The faint emission under scrutiny is completely buried in the noise in the Hö-CLEAN reconstruction (Figure 1, top). A blurry and noisy depiction emerges in the reconstructions of both CS-CLEAN (Figure 1, bottom) and MS-CLEAN (Figure 2, top).

The hotspots. As recovered by R2D2 variants, AIRI and uSARA, the hotspots highlight the ability of these algorithms to resolve physical structure beyond instrumental resolution, in contrast with CLEAN variants (see panels (c) and (d) of all figures). Interestingly, where AIRI and uSARA exhibit artificial zero-valued pixels around the hotspots, all R2D2 variants depict continuous emission. This observation suggests the ability of R2D2 variants to achieve a more physical reconstruction.

Image-domain data fidelity. We evaluate data fidelity delivered by the imaging algorithms by scrutinizing their residual dirty images (i.e., back-projected data residual) displayed in panels (a) of all figures, whose standard deviation values are reported in the corresponding captions. AIRI, uSARA, and CS-CLEAN obtain residual dirty images with the lowest standard deviation values, reflecting their high data fidelity. Both MS-CLEAN and R2D2 variants provide comparable values, slightly above the best-performing algorithms. Among R2D2 variants, R3D3 performs better than both R2D2 and R2D2-Net, owing to its underpinning model-informed DNN modules on the one hand, and its iterative structure on the other hand. Hö-CLEAN delivers the lowest fidelity with its standard deviation value being more than 1 order of magnitude higher than the rest. This numerical analysis aligns with the overall visual examination. A closer inspection of the R2D2 variants reveals discernible structure at the pixel positions of the western and eastern hotspots, where the highest emission is concentrated. While similar behavior has been observed in high-dynamic-acquisition regimes in simulation (Aghabiglou et al. 2024), these patterns are possibly amplified by the pointing errors at the hotspots.

It is worth noting that the quest for noise-like data residuals is only really meaningful in combination with the recovery of a model satisfying physical constraints. The fact that CLEAN variants accept negative components eases the reduction of their data residuals. On the contrary, R2D2 variants, uSARA, and AIRI explicitly enforce the positivity of the recovered intensity images.

Computational cost. Computational details of the imaging algorithms are provided in Table 2. R2D2 variants are able to deliver reconstructions in a few seconds only. We note that reported times include a computational overhead for initializing DNN modules at the first iteration, inducing a large difference between the average time of the regularization step between R2D2-Net on one hand, and R2D2 and R3D3 on the other. Both CS-CLEAN and MS-CLEAN require around half a minute, whereas the highly iterative algorithms AIRI and uSARA take about 10 and 20 minutes, respectively. When compared to CS-CLEAN and MS-CLEAN, both R2D2 and R3D3 perform comparable number of passes through the nongridded visibilities. Their DNN modules provide a boost to their speed in comparison with the minor cycles of the advanced CLEAN variants. R2D2-Net is also at least 1 order of magnitude faster than its CLEAN counterpart Hö-CLEAN. However, one must recognize the differences in implementation and hardware between the imaging algorithms, reflected in the varying average time of their data-fidelity step (all calling for the computation of the residual dirty images).

Of note, the training of R2D2 variants required the order of thousands of CPU core hours and GPU hours (Aghabiglou et al. 2024). This is, however, a one-off cost accommodating any number of subsequent reconstructions from VLA data, as exemplified by the fact that the present analysis did not require any training.