Benchmarking

Benchmarking Platform for State-of-the-Art RI Imaging Algorithms

We present a benchmarking platform for the state-of-the-art algorithms for synthesis imaging by interferometry (RI) provided by BASPLib: uSARA, AIRI, and R2D2. We provide access to our simulated comprehensive testbed that is composed of 200 RI inverse problems, and offer detailed guidelines to reproduce the reconstruction results published in our study [1].

RI inverse problem

In radio interferometry, each pair of antennas acquires a noisy complex measurement corresponding to a spatial Fourier component of the target radio sky, at a given time instant. The collection of the sensed spatial Fourier modes accumulated over the total observation period forms the Fourier sampling pattern, describing an incomplete coverage of the 2D Fourier plane. Considering the unknown radio image of interest \({\boldsymbol{x}^{\star}} \in \mathbb{R}_{+}^{N},\) the RI data model reads: \[ {\boldsymbol{y}} = {\boldsymbol{\Phi} {\boldsymbol{x}^{\star}} + \boldsymbol{n}} \] where \({\boldsymbol{y}} \in \mathbb{C}^M\) denotes the data vector, \({\boldsymbol{n}} \in \mathbb{C}^M\) represents a complex Gaussian random noise with a standard deviation \(\tau>0\) and mean 0. The linear operator \({\boldsymbol{\Phi}}\colon \mathbb{R}^N \to \mathbb{C}^M\) is the measurement operator (i.e., describing the acquisition process), consisting in a nonuniform Fourier sampling. The operator \({\boldsymbol{\Phi}}\) can also include a data-weighting scheme to improve the observation's effective resolution.

The radio-interferometric data model can be expressed in the image domain via a normalised back-projection using the adjoint of the measurement operator \(\boldsymbol{\Phi}^\dagger \) such that: \[\boldsymbol{x}_{\textrm{d}} =\kappa \text{Re}\{\boldsymbol{\Phi}^\dagger \boldsymbol{y}\}= \kappa \text{Re}\{\boldsymbol{\Phi}^\dagger \boldsymbol{\Phi}\boldsymbol{x}^{\star}+ \boldsymbol{\Phi}^\dagger \boldsymbol{n}\}\] where \(\boldsymbol{x}_{\textrm{d}} \in \mathbb{R}^N\) stands for the back-projected data, known as the dirty image, and \(\textrm{Re}\{\cdot\}\) stands for the real part of its argument. In this equation, the normalisation factor \(\kappa>0\) is equal to \((\max(\text{Re}\{{\boldsymbol{\Phi}}^{\dagger}{\boldsymbol{\Phi}}\boldsymbol{\delta}\}))^{-1}\), ensuring that the point spread function (PSF), defined as \(\kappa \text{Re}\{\boldsymbol{\Phi}^{\dagger}{\boldsymbol{\Phi}}\boldsymbol{\delta}\} \in \mathbb{R}^N\), has a peak value of 1. The Dirac \(\boldsymbol{\delta} \in \mathbb{R}^N\) is the image with value 1 at its center and 0 elsewhere.

RI testbed description

The testbed is composed of ground-truth images (in .fits format) and associated data files (in .mat format). The ground-truth images are derived from four real radio images, which underwent post-processing to eliminate artefacts such as deconvolution errors and noise, and endow them with high dynamic ranges.

The considered raw radio images are the giant radio galaxies 3C353 and M106, and the radio galaxy clusters A2034 and PSZ2G165.68+44.01, sourced from NRAO Archives and LOFAR surveys [2,3,4]. From each radio image, 50 curated ground-truth images of size 512\(\times\)512 were generated, characterised with large dynamic ranges varying in [10³, 5\(\times\)10⁵]. Simulated radio interferometric data associated with the ground-truth images are obtained from sampling patterns of the Very Large Array (VLA), generated under a wide variety of the observation settings, including total observation time, pointing direction, etc. Sampling patterns combining both VLA configurations A and C were considered for a balanced Fourier coverage. The total number of points in the resulting Fourier sampling patterns ranges from 0.2M to 2M. The data are corrupted with additive Gaussian noise with mean 0 and a standard deviation commensurate of the dynamic range of the corresponding ground truth images.

The simulated data are such that the spatial bandwidth of the Fourier sampling is fixed to 2/3 of the spatial Fourier bandwidth of the ground-truth images. Under these considerations, the testbed comprises 200 unique inverse problems. During imaging, the Briggs weighting scheme was applied to the data (also injected into the measurement operator) to enhance the effective resolution of the reconstructions, where the Briggs parameter was uniformly randomised between \(-1\) and \(1\), with lower values approaching uniform weighting, and higher values approaching natural weighting. The pixel resolution of the dirty images was also randomly chosen to reflect a super-resolution factor during imaging in the range \([1.5, 2.5]\).

Detailed information on the testbed can be found in [1]. The code utilised to generate a wide variety of sampling patterns and corresponding RI data is part of the utilities offered by BASPLib.

RI imaging algorithms

The BASPLib imaging algorithms used for benchmarking are:

uSARA: the unconstrained counterpart of the SARA algorithm [5,6]. It promotes a handcrafted sparsity-based regularisation, and is underpinned by the forward-backward algorithmic structure. The Python implementation of uSARA for small-scale RI imaging is available in uSARA repository.

AIRI: a Plug-and-Play (PnP) algorithm at the intersection of optimisation theory and deep learning. By inserting carefully trained AIRI denoisers into the proximal splitting algorithms, one waives the computational complexity of optimisation algorithms induced by sophisticated image priors. AIRI is the PnP counterpart of uSARA. The Python implementation of AIRI for small-scale RI imaging is available in AIRI repository.

R2D2: a novel paradigm for RI, interpreted as a learned version of the well-known Matching Pursuit algorithm [1,7]. R2D2 reconstruction is formed as a series of residual images, iteratively estimated as outputs of iteration-specific DNNs, each taking the previous iteration’s image estimate and associated back-projected data residual as inputs. We incorporate two distinct DNN architectures as the core components of R2D2:

R2D2\(_{\mathcal{A}_1}\):This configuration employs the widely-used U-Net architecture for its DNNs.

R2D2\(_{\mathcal{A}_2}\):This variant leverages a more advanced architecture, U-WDSR, which enhances performance.

The Python implementation of both R2D2\(_{\mathcal{A}_1}\) and R2D2\(_{\mathcal{A}_2}\) are available in R2D2 repository.

U-Net / U-WDSR: an end-to-end DNN, which also corresponds to the first term of the R2D2 DNN series.

In addition to BASPLib imaging algorithms, the benchmarking results also include the celebrated CLEAN algorithm. CLEAN, recognised as a standard algorithm for RI imaging, is represented here by its multiscale variant implemented in the widely used WSClean software [8].

Evaluation metrics

To assess the performance of the reconstruction algorithms, we use quantitative metrics described below.

Ground truth fidelity (SNR,logSNR)
The Signal-to-Noise Ratio (SNR) metric is defined as: \[\text{SNR}(\boldsymbol{x}^{\star}, \boldsymbol{x}) = 20\log_{10}\left(\frac{\|\boldsymbol{x}^{\star}\|_{2}}{\|\boldsymbol{x}^{\star} - \boldsymbol{x}\|_{2}} \right), \] where \(\boldsymbol{x}^{\star}\) stands for the ground truth, and \(\boldsymbol{x}\) denotes the reconstructed image.

The logSNR metric emphasises faint structures by applying a logarithmic mapping to the images, parametrised by the target dynamic range \(a\): \[ \text{r}\log(\boldsymbol{x}) = x_{\text{max}} \log_{a} \left( \frac{a}{x_{\text{max}}} \cdot \boldsymbol{x} + 1 \right), \] where \(x_{\text{max}}\) indicates the peak pixel value of the image \(\boldsymbol{x}\). Using this mapping, the logSNR metric is defined as: \[ \log\text{SNR}(\boldsymbol{x}^{\star}, \boldsymbol{x}) = \text{SNR}(\text{rlog}(\boldsymbol{x}^{\star}), \text{rlog}(\boldsymbol{x})). \]

Image-domain data fidelity RDR
Image-domain data fidelity metric, known as the residual-to-dirty image ratio (RDR) is assessed using the estimated residual dirty image defined as \(\mathbf{r}\)=\(\boldsymbol{x}_{\textrm{d}} - \kappa \text{Re}\{\boldsymbol{\Phi}^\dagger \boldsymbol{\Phi}\boldsymbol{x}\}\). The RDR metric is defined as: \[ \text{RDR}(\mathbf{r}, {\boldsymbol{x}_{\textrm{d}}}) = \frac{\|\mathbf{r}\|_{2}}{\|\boldsymbol{x}_{\textrm{d}}\|_{2}} \]

Table of Results

The numerical reconstruction results obtained by state-of-the-art imaging algorithms are reported in the table below. The reconstruction quality is assessed using metrics such as SNR, logSNR, and RDR. Computational performance is detailed in terms of the total number of iterations (I) and the total reconstruction time (t_tot). Please note that reported timings depend on the machine used and are provided for reference only. \([\boldsymbol{\Phi}^\dagger \boldsymbol{\Phi}]_{imp.}\) is indicating the measurement operator implementation. The subscript \(\mathcal{T}_2\) in R2D2 varients refer to the generlised VLA-specific training set that is used to train these R2D2 algorithms. All reported values represent mean \(\pm\) standard deviation, calculated over 200 inverse problems.

Algorithm	SNR (dB)	logSNR (dB)	RDR (×10-3)	I	ttot. (s)	\([\boldsymbol{\Phi}^\dagger \boldsymbol{\Phi}]_{imp.}\)
CLEAN	12.0±19.3	9.4±18.9	3.29±28.2	8.4±1.0	106.5±81.6	-
uSARA	28.1±3.4	20.4±3.4	2.15±27.5	1482.2±586.1	368.8±296.8	TorchKbNufft
				1490.7±563.2	216.7±119.3	PyNUFFT
				1483.1±587.0	103.0±39.52	FINUFFT
				1482.6±586.4	88.03±31.98	PSF
AIRI	28.3±3.1	21.1±3.8	2.24±28.0	5000.0±0.0	937.4±801.8	TorchKbNufft
		21.0±3.8			566.5±355.8	PyNUFFT
					157.0±36.92	FINUFFT
					114.2±3.450	PSF
U-Net	17.9±3.0	6.8±3.9	113.3±589.2	1	0.641±0.110	-
U-WDSR	16.0±3.6	6.6±3.8	155.8±815.5	1	0.662±0.031	-
R2D2\({\mathcal{A}_1,\mathcal{T}_2}\)	30.0±3.0	23.4±4.2	4.07±91.6	18.3±5.6	7.243±4.131	TorchKbNufft
					6.932±3.960	PyNUFFT
					3.771±1.224	FINUFFT
					3.342±0.931	PSF
R2D2\({\mathcal{A}_2,\mathcal{T}_2}\)	31.2±2.4	24.6±4.2	2.22±28.1	15.8±5.5	8.831±3.923	TorchKbNufft
					9.023±4.112	PyNUFFT
					5.951±2.199	FINUFFT
					5.649±1.911	PSF

Reconstructions & Residual dirty images

Reconstruction images obtained by the RI algorithms associated with one inverse problem selected from our testbed (ID:80009) are displayed below. In this case, the ground truth is characterised by a dynamic range of 150,000.

Estimated model images:

Residual dirty images:

How to Use the Benchmarking Platform

In what follows, we provide instructions to reproduce the results summarised in the above table.
Setup
Setting up the main directory
Firstly, set up the main directory, which will comprise the benchmark testbed, RI imaging repositories, and pre-trained DNNs (for AIRI and R2D2). Replace /path/to/your/directory with your local path. Subdirectories will be created accordingly.

Downloading testbed
Secondly, download the testbed zip files. These correspond to compressed ground truth files in .fits and data files in .mat format.

Thirdly, extract the compressed files in the dedicated subdirectory.

You can access each inverse problem details via this Excel file.

Setting up imaging algorithms repositories

Python environment

Firstly, create your Python environment using either venv or conda, with Python version >= 3.10.0. In this tutorial, we will be using venv and pip.

Secondly, activate this new environment.

The uSARA, AIRI and R2D2 repository is Pytorch-based for accelerated computation on GPU. To ensure compatibility of your CUDA driver version, follow PyTorch instruction. If your CUDA version is older than what is displayed, please refer to here for the newest PyTorch version compatible for your GPU.

Thirdly, once PyTorch is installed, install the remaining dependencies.

uSARA & AIRI

Clone the repository to the current directory by running the command below.

The --recurse-submodules flag is required so that the required sub-modules RI-measurement-operator will be cloned at the same time.

Pre-trained AIRI denoisers

Pre-trained AIRI denoisers are in .onnx format. From main AIRI directory run commands below to download the denoiser shelf compressed in a .tar.gz file, and extract the DNNs in the appropriate subdirectory.

Two additional collections of pre-trained AIRI denoisers are available [9]. While they are not required for this benchmarking, they can be utilisd if needed.

R2D2

The R2D2 repository encompasses the end-to-end DNNs U-Net, U-WDSR, R2D2 variants: R2D2_{\(\mathcal{A}_1\)}, and R2D2_{\(\mathcal{A}_2\)}. The algorithm of choice is selected from the input parameters of the imager. The currently available R2D2 models are trained specifically for the imaging settings considered in [1]. Deviating from these settings may result in sub-optimal results.
First clone the repository.

Download the pretrained R2D2_{\(\mathcal{A}_1\)} model.

Download the pretrained R2D2_{\(\mathcal{A}_2\)} model.

CLEAN

We use the multi-scale variant of CLEAN in WSClean. Please refer to WSClean official website for full instructions on how to install the software.

Imaging
Instructions to run BASPLib RI imaging algorithms and WSClean are outlined. We recall the imaging settings are as follows: (i) the target image is of size 512x512 pixels, with (ii) a random pixel size corresponding to super-resolution factor between [1.5, 2.5], and (iii) random Briggs weighting with parameter with range between [-1, 1] is applied to the data. All algorithms take as input a configuration file in either JSON (for uSARA and AIRI) or YAML format (for U-Net, U-WDSR and R2D2 variants: R2D2_{\(\mathcal{A}_1\)}, and R2D2_{\(\mathcal{A}_2\)}). Below are example Bash scripts designed to image all inverse problems at once. Users are encouraged to customise these scripts to suit the specifications of their local machines as needed.

uSARA

For uSARA imaging algorithm, all inverse problems were solved with fixed parameters. For best results, specifically the regularisation parameter was set to two times its heuristic value. You can access sample config file for uSARA in uSARA.json.
uSARA imager can be run by calling:

AIRI

The AIRI parameter was fixed to the heuristic value for all RI data, except for those simulated from the radio image of 3C353, where 3 times the heuristic value was used. You can access sample config files for AIRI in AIRI.json and
AIRI imager can be run by calling:

R2D2

You can access sample config file for R2D2 in R2D2.yaml.
The R2D2 imager can be run by calling:

For running R2D2 imager save above bash script into a bash file like R2D2.sh, and run like:

The script will prompt you to enter the R2D2 underpinned DNN architecture name. For example, you can enter U-Net, U-WDSR. The script will then execute the corresponding configuration for the specified algorithm. Later, user need to enter the the number of R2D2 iteration needed. 25 iterations are available. For U-Net (resp. U-WDSR) the number of input channels should be input as 64 (resp. 32). BASPLib provides four options for implementing the RI measurement operator model. Three of these consist in different implementations of the NUFFT including TorchKbNufft [10], FINUFFT [11], and PyNUFFT [12]. The fourth option leverages the PSF, which under the assumption of a narrow field of view, enables approximating the RI mapping operator \(\boldsymbol{\Phi}^\dagger \boldsymbol{\Phi}\) via a convolution with the PSF.

CLEAN

For CLEAN, specifically its multi-scale variant all inverse problems were solved with fixed parameters. WSClean imager can be run by calling:

Note: Additional information on how to convert the data files into a measurement set that is readable by WSClean will be provided soon.

Evaluation metrics calculation

The evaluation metrics can be calculated using this Python file. To run calculate_metrics.py you can use below bash script:

References

[1] Aghabiglou A, Chu CS, Tang C, Dabbech A, Wiaux Y. Towards a robust R2D2 paradigm for radio-interferometric imaging: revisiting DNN training and architecture. arXiv:2503.02554v1

[2] Edler HW, de Gasperin F, Shimwell TW, Hardcastle MJ, Boselli A, Heesen V, McCall H, Bomans DJ, Brüggen M, Bulbul E, Chyży KT. ViCTORIA project: The LOFAR HBA Virgo Cluster Survey. Astronomy & Astrophysics. 2023 Aug 1;676:A24. | DOI:10.1051/0004-6361/202346458

[3] Shimwell TW, Hardcastle MJ, Tasse C, Best PN, Röttgering HJ, Williams WL, Botteon A, Drabent A, Mechev A, Shulevski A, van Weeren RJ. The LOFAR two-metre sky survey-V. Second data release. Astronomy & astrophysics. 2022 Mar 1;659:A1. | DOI:10.1051/0004-6361/202142484

[4] Botteon A, Shimwell TW, Cassano R, Cuciti V, Zhang X, Bruno L, Camillini L, Natale R, Jones A, Gastaldello F, Simionescu A. The Planck clusters in the LOFAR sky-I. LoTSS-DR2: new detections and sample overview. Astronomy & Astrophysics. 2022 Apr 1;660:A78. | DOI:10.1051/0004-6361/202143020

[5]Carrillo RE, McEwen JD, Wiaux Y. Sparsity Averaging Reweighted Analysis (SARA): a novel algorithm for radio-interferometric imaging. Monthly Notices of the Royal Astronomical Society. 2012 Oct 21;426(2):1223-34. | DOI:10.1111/j.1365-2966.2012.21605.x

[6] Onose A, Dabbech A, Wiaux Y. An accelerated splitting algorithm for radio-interferometric imaging: when natural and uniform weighting meet. Monthly Notices of the Royal Astronomical Society. 2017 Jul;469(1):938-49. | DOI:10.1093/mnras/stx755

[7] Dabbech A, Aghabiglou A, San Chu C, Wiaux Y. CLEANing Cygnus A deep and fast with R2D2. The Astrophysical Journal Letters. 2024 May 7;966(2):L34. | ArXiv:2309.03291 | DOI:10.3847/2041-8213/ad41df

[8] Offringa AR, Smirnov O. An optimised algorithm for multiscale wideband deconvolution of radio astronomical images. Monthly Notices of the Royal Astronomical Society. 2017 Oct;471(1):301-16. | DOI:10.1093/mnras/stx1547

[9] Terris M, Tang C, Jackson A, Wiaux Y. Plug-and-play imaging with model uncertainty quantification in radio astronomy. arXiv preprint arXiv:2312.07137. 2023 Dec 12.

[10] Muckley, M. J., Stern, R., Murrell, T., & Knoll, F. TorchKbNufft: A High-Level, Hardware-Agnostic Non-Uniform Fast {Fourier} Transform. ISMRM Workshop on Data Sampling & Image Reconstruction. 2020 May.

[11] Shih, Y.-h., Wright, G., And´en, J., Blaschke, J., & Barnett, A. H: cuFINUFFT: a load-balanced GPU library for general-purpose nonuniform FFTs. 2021 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), 688–697 | ArXiv:2102.08463

[12] Lin, J.-M. Python non-uniform fast Fourier transform (PyNUFFT): An accelerated non-Cartesian MRI package on a heterogeneous platform (CPU/GPU). 2018, Journal of Imaging. 4, 51. | ArXiv:1710.03197 | DOI:10.3390/jimaging4030051

Support

If you need assistance with setting up the BASPLib imaging algorithms or encounter any issues during installation, please don’t hesitate to contact our support team.