Scalable precision wide-field imaging in radio interferometry: I. uSARA validated on ASKAP data

As Part I of a paper series showcasing a new imaging framework, we consider the recently proposed unconstrained Sparsity Averaging Reweighted Analysis (uSARA) optimisation algorithm for wide-field, high-resolution, high-dynamic range, monochromatic intensity imaging. We reconstruct images from real radio-interferometric observations obtained with the Australian Square Kilometre Array Pathfinder (ASKAP) and present these results in comparison to the widely-used, state-of-the-art imager WSClean. Selected fields come from the ASKAP Early Science and Evolutionary Map of the Universe (EMU) Pilot surveys and contain several complex radio sources: the merging cluster system Abell 3391-95, the merging cluster SPT-CL 2023-5535, and many extended, or bent-tail, radio galaxies, including the X-shaped radio galaxy PKS 2014-558 and the ``dancing ghosts'', known collectively as PKS 2130-538. The modern framework behind uSARA utilises parallelisation and automation to solve for the w-effect and efficiently compute the measurement operator, allowing for wide-field reconstruction over the full field-of-view of individual ASKAP beams (up to 3.3 deg each). The precision capability of uSARA produces images with both super-resolution and enhanced sensitivity to diffuse components, surpassing traditional CLEAN algorithms which typically require a compromise between such yields. Our resulting monochromatic uSARA-ASKAP images of the selected data highlight both extended, diffuse emission and compact, filamentary emission at very high resolution (up to 2.2 arcsec), revealing never-before-seen structure. Here we present a validation of our uSARA-ASKAP images by comparing the morphology of reconstructed sources, measurements of diffuse flux, and spectral index maps with those obtained from images made with WSClean.


INTRODUCTION
We are now at the beginning of a booming era for radio astronomy, with a worldwide effort to gear up for the Square Kilometre Array (SKA) -a revolutionary radio telescope with capabilities for sub-arcsecond resolution and ultra-deep sensitivity. Pathfinding radio interferometers -such as the Murchison wide-field Array (MWA; Tingay et al. 2013), the LOw Frequency ARay (LOFAR; van Haarlem et al. 2013), ASKAP (Johnston et al. 2007(Johnston et al. , 2008Hotan et al. 2021), and MeerKAT (Jonas & MeerKAT Team 2016) -are paving the way by affording us the opportunity to expand our capabilities in detection, calibration, and image reconstruction of the unknown radio sky. Ongoing wide-field radio continuum surveys -such as the LO-FAR Two-meter Sky Survey (LoTSS; Shimwell et al. 2017Shimwell et al. , 2019, the LOFAR LBA Sky Survey (LoLSS; de Gasperin et al. 2021), the MeerKAT MIGHTEE survey (Taylor & Jarvis 2017), and the ASKAP EMU survey (Norris et al. 2011) -will be used in conjunction to gather statistics on millions of radio galaxies and thousands of galaxy clusters, leading to new and exciting results on cosmic magnetism, cosmic rays, dark matter, dark energy, and the evolution of ★ E-mail: y.wiaux@hw.ac.uk large-scale structure in the Universe. The quest to convert the sheer quantity of radio data from these surveys into science-ready images has propelled radio astronomers toward developing innovative, stateof-the-art calibration (e.g. Smirnov & Tasse 2015;Williams et al. 2016;van Weeren et al. 2016) and imaging (e.g. Offringa & Smirnov 2017;Tasse et al. 2018;Pratley et al. 2018) techniques throughout the last decade. Through the implementation of these trailblazing techniques, results from LOTSS, the ASKAP-EMU Pilot Survey (Norris et al. 2021), and MIGHTEE are already revealing extraordinary extragalactic radio sources that have not been previously detected, some of which are so complex in their physical properties that they challenge current taxonomies and origination theories (e.g. Botteon et al. 2020;Brüggen et al. 2021;Knowles et al. 2022).
Real intensity structure in radio sources -such as Galactic HI emission, supernova remnants, extended radio galaxies, and galaxy cluster radio halos and relics -often exhibits both compact and diffuse components. A complex radio source might include bright, collimated threads of emission that appear embedded within fainter, dispersed lobes. For instance, intracluster radio sources, generated by large-scale turbulence and shocks, tend to exhibit high-resolution filamentary structures (tracing well-aligned magnetic field lines) as well as low-surface-brightness diffuse structure, which can span Mega-parsec scales (see the radio relic of Abell 2256 as one iconic example e.g. Owen et al. 2014;Rajpurohit et al. 2022). When imaging such complex radio emission, state-of-the-art CLEAN-based algorithms (e.g. Högbom 1974;Clark 1980;Schwab 1984;Wakker & Schwarz 1988;Cornwell 2008;Offringa et al. 2014) are usually implemented with data weighting schemes to adjust the sensitivity to either compact or diffuse emission. By modifying the shape of the synthesised beam and placing weights (such as a -taper) to short-baseline data, the sensitivity to fainter, more extended components of radio emission can be enhanced, albeit with a considerable loss of resolution. Consequently, such imaging methods often fail to accurately reconstruct both compact and diffuse components simultaneously in a single image.
In the last decade, progress in compressed sensing research has led to the development of optimisation-based algorithms to reconstruct true signal from partial, undersampled visibility data. This methodology was first applied and shown to be effective for radio interferometry imaging by Wiaux et al. (2009) and has since led to the development of several designated imaging algorithms (e.g. Li et al. 2011;Carrillo et al. 2012;Garsden et al. 2015;Dabbech et al. 2015). Such innovative approaches -relying on sophisticated sparsity-based image models -have demonstrated their success at capturing both compact and diffuse components in the reconstructed radio images, albeit with an increase in computational expense. One such stateof-the-art optimisation-based class of methods is the "Sparsity Averaging Reweighted Analysis" (SARA) family (Carrillo et al. 2012;Dabbech et al. 2018;Abdulaziz et al. 2019; Thouvenin et al. 2022a, see the Puri-Psi web-page for more details). The monochromatic SARA algorithm, initially proposed by Carrillo et al. (2012), leverages state-of-the-art algorithmic structures from optimisation theory to enforce non-negativity and sparsity in a highly redundant sparsity dictionary of the sought radio image under noise-driven data fidelity constraints (Onose et al. 2016;Onose et al. 2017). Evolved versions of SARA include Faceted HyperSARA for wide-band imaging, shipped with spectral and spatial faceting functionalities to handle large image dimensions (Thouvenin et al. 2022a,b), Polarized SARA for polarisation imaging (Birdi et al. 2018), and a sparsity-based joint calibration and imaging framework Birdi et al. 2020;Dabbech et al. 2021).
In addition to the precision and robustness requirements of RI imaging algorithms, the need for scalability is more critical than ever in light of the extreme data volumes, wide fields-of-view, and broad frequency bandwidths offered by modern radio arrays. In this context, we have recently proposed a parallel and automated framework for wide-field, high-resolution, high-dynamic range monochromatic intensity imaging (Dabbech et al. 2022). The framework encapsulates two imaging algorithms at the interface of optimisation theory and deep learning, recently proposed by Terris et al. (2022): (i) the unconstrained SARA (uSARA) algorithm relying on a handcrafted image model enforced via an optimisation-based denoiser, and (ii) the AI for Regularization in Imaging (AIRI) algorithm relying on an implicitly learnt image model promoted via denoising deep neural networks. The framework offers scalable implementation of both imaging algorithms through two key features: memory-efficient, parallel operators and automated parameter selection.
This article is Part I of a series aiming to showcase and validate the uSARA algorithm on imperfectly calibrated RI observations from the ASKAP radio telescope. In Part II, we expand upon this work to include a validation of the AIRI algorithm. In both articles, we aim to study the imaging performance of uSARA and AIRI in comparison to Multi-scale CLEAN (Cornwell 2008) via the WSClean imager (Offringa et al. 2014), both in terms of reconstruction quality and computational efficiency. For a coherent comparative analysis of the three imaging algorithms -uSARA, AIRI, and WSClean -both Part I and Part II utilise the same RI data from publicly available ASKAP observations collected during Early Science and Pilot surveys and processed through the ASKAPsoft pipeline (Hotan et al. 2021). Targeted fields-of-view -hosting extended, diffuse, and complex radio sources -were carefully selected to test the precision and scalability capabilities of our imaging framework. A comprehensive summary of the considered imaging framework (with interchangeable denoising algorithms), including the wide-field measurement operator model, its distribution through parallelisation, and its implementation using high-performance computing systems (HPC), is provided by Dabbech et al. (2022). A fully detailed analysis of the framework's scalability is the subject of a forthcoming article.
The remainder of this article is structured as follows. In Section 2, we present an overview of the investigated imaging framework from the algorithmic structure underpinning uSARA to the parallelisation and automation functionalities ensuring the computational scalability of the framework to large image and data dimensions in the context of wide-field imaging. In Section 3, we provide details of the scrutinised ASKAP data and the imaging settings of uSARA and the CLEANbased benchmark algorithm. Reconstruction results of our primary targets of interest are exposed in Section 4 and discussed in Section 5. Section 6 documents and discusses the computational cost of the imaging framework. Finally, conclusions are drawn in Section 7.

METHODS
In this section, we present the RI data model in the context of widefield monochromatic intensity imaging. We provide an overview of the uSARA imaging algorithm and its underpinning algorithmic structure, and summarise the scalability features of its encompassing framework (Terris et al. 2022;Dabbech et al. 2022).

Data Model
In the absence of instrumental and atmospheric perturbations, RI visibilities measured at a given observing wavelength are noisy Fourier components of the radio sky, modulated by the so-called -effect, a varying chirp-like phase induced by the non-coplanarity of the radio array. With no loss of generality, the data model can be discretised, such that the measured visibilities ∈ C are modelled from the sought intensity image ∈ R + as follows where ∈ C is a realisation of a zero-mean random Gaussian noise with a standard deviation > 0. The operator ∈ C × is the measurement operator encompassing the Fourier sampling, and the -effect. Due to the large amounts of data, a Direct Fourier transform would be intractable; therefore, the incomplete Fourier sampling is modelled via the non-uniform fast Fourier transform (NUFFT) (Fessler & Sutton 2003;Onose et al. 2016). Furthermore, the -effect is taken into account via a hybrid model combining -stacking (Offringa et al. 2014) and -projection , whereby RI data are grouped by their -coordinates into -stacks. For each data point, the chirp-like phase is decomposed into two terms: a modulation of its associated -stack injected in the measurement operator via image-domain multiplication, and a compact Fourier kernel encoding the resulting -offset modulation, injected via Fourier-domain convolution . As a final consideration, a noise-whitening operation is typically applied to the measured data and injected into the associated measurement operator to ensure constant standard deviation of the noise (see Appendix A of Terris et al. 2022, for more details). On some occasions, the operation is performed in combination with a data-weighting scheme derived from the sampling profile (e.g. Briggs weighting;Briggs 1995) to improve the effective resolution of the observation (due to the highly non-uniform density profile of the RI sampling).
Under these considerations, the measurement operator is decomposed into computationally and memory-efficient blocks as the vertical concatenation of the operators 1≤ ≤ , where for each -stack ∈ {1, . . . , }, the associated measurement operator ∈ C × is given by = G FZ (Dabbech et al. 2022). More specifically, the operator ∈ R × is a diagonal matrix encoding the considered data-weighting scheme. The sparse matrix G ∈ C × is the de-gridding matrix, encompassing convolutions between the NUFFT interpolation kernels and the compact -kernels correcting for the associated -offsets in the Fourier plane. Note that estimates of direction-dependent effects (DDEs) can also be encoded as additional convolutions in the rows of the de-gridding matrix, when available. F ∈ C × is the Discrete Fourier transform and the operator Z ∈ C × encodes the -modulation of the th -stack, the zero-padding operator for a finer grid of the Fourier plane, and the correction for the convolution with the approximate NUFFT interpolation kernels.

uSARA algorithm
Image formation from the noisy and incomplete RI measurements is an ill-posed inverse problem. Here we consider the unconstrained SARA imaging algorithm from optimisation theory (Terris et al. 2022). The algorithm provides an estimate of the radio sky as the minimiser of an objective function posed as the sum of two terms: a data fidelity term , emanating from the nature of the noise, and a regularisation term encoding a prior knowledge of the image to address the ill-posedness of the inverse problem. The minimisation task is of the form where > 0 is the regularisation parameter controlling the balance between the two terms. Given the Gaussian nature of the noise affecting the RI data, is naturally set to ( ; ) = 1/2 − 2 2 , with . 2 denoting the ℓ 2 norm of its argument vector.
The uSARA regularisation function is a multi-term nondifferentiable function composed of a non-convex log-sum function enforcing average sparsity in an overcomplete dictionary ∈ R × , consisting in the normalised concatenation of nine orthogonal bases, and a non-negativity constraint (Carrillo et al. 2012;Terris et al. 2022), which reads where (.) denotes the th coefficient of its argument vector, and (·) † stands for the adjoint of its argument operator. The parameter > 0 prevents the argument of the logarithmic from reaching zero values and can be set to the estimate of the noise level in the sparsity domain (Thouvenin et al. 2022a). The non-negativity constraint is encoded via the indicator function of the real positive orthant, given by R + ( ) = +∞ if ∉ R + and 0 otherwise. As such, the resulting minimisation task is non-convex and is addressed in an iterative manner. More specifically, the problem is approached by solving a sequence of surrogate convex minimisation tasks whereby is replaced by a convex regularisation function of the form ( ) = W † 1 + R + ( ) substituting the log-sum function with the ℓ 1 function, denoted by · 1 , and weighted by the diagonal matrix W ∈ R × . Each of the surrogate weighted minimisation tasks is the form where the convex and non-differentiable function is redefined through the update of its underlying weighting matrix W from the solution of its preceding task (Carrillo et al. 2012;Terris et al. 2022). The convex minimisation task is solved approximately (i.e. for a finite number of iterations > 0) using the forward-backward (FB) iterative scheme Terris et al. 2022), and the overall procedure benefits from convergence guarantees ).
The FB iterative scheme relies on two-step image updates: a 'forward' gradient descent step calling for the gradient of the data fidelity function , given by ∇ ( ) = Re{ † } − Re{ † }, followed by a 'backward' denoising step using the proximal operator of the convex regularisation function (see Terris et al. 2022, for the mathematical details) such that for any ∈ N ( +1) = prox ( ) − ∇ ( ( ) ) .
Let > 0 denote the Lipschitz constant of ∇ given by = Re{ † } , with . denoting the spectral norm of its argument operator. The step-size satisfies the condition 0 < < 2/ to guarantee the convergence of the iterative scheme. Finally, the proximal operator prox , not benefiting from closed-form solutions, is computed sub-iteratively, involving soft-thresholding operations in the sparsity dictionary by .

A scalable and automated imaging framework
To address the scalability requirement to large data and image sizes, our imaging framework provides automated parameter choice and parallel and memory-efficient models of the operators and functions involved (Dabbech et al. 2022), summarised in what follows.
Regularisation parameter selection. The choice of the regularisation parameter , balancing data fidelity and image regularisation, is of paramount importance as it affects the solution of the minimisation task (2). Considering > 0, the estimate of the standard deviation of the image domain noise, Terris et al. (2022) proposed to equate the soft-thresholding parameter , involved in the denoising step, to the estimate of the standard deviation of the sparsity domain noise given by /3 (the factor three emanates from the normalisation of the sparsity dictionary). In the case when a data-weighting scheme is adopted to compensate for the non-uniform density profile of the sampling (e.g. Briggs weighting), additional correlation is induced in the image domain noise. Under this consideration, can be obtained as where the data-weighting operator ∈ R × is a diagonal per block matrix, whose diagonal matrix-blocks are the data-weighting matrices 1≤ ≤ . The correction factor reduces to one otherwise. In our experiments, Briggs weighting was applied to the data in imaging, and the resulting values of were found to be in the interval [0.3, 0.6]. The regularisation parameter can be set around with the step size fixed to = 1.98/ , ensuring the convergence of the FB algorithm. Finally, the parameter involved in (3) is typically set to the estimate of the standard deviation of the sparsity domain noise, = /3.

Denoiser Faceting.
In light of the sub-iterative nature of the denoising operator underpinning uSARA, distribution and parallelisation of the sparsity operator are required, not only to handle the large image dimensions of interest but also to accelerate the computation. For this aim, we have adopted a faceted implementation of the operator (Pruša 2012), enabling image facet denoising. The number of facets is derived from the number of CPU cores of the computing architecture on which the algorithm is to be deployed, and constraints on image facet dimensions from the wavelet transforms underpinning the sparsity operator .
Automated parallelisation of the measurement operator. Three key features are supported in our implementation of the measurement operator to ensure its scalability to large data sizes (Dabbech et al. 2022). Firstly, the choice of the number of -stacks, , defining the decomposition of the measurement operator into the operators is automated via a planning strategy taking into consideration the computational cost derived from the complexity of the application of the measurement operator and the memory constraints of the computing architecture. Secondly, memory-efficient encoding of the resulting operators † , called for in FB, can be achieved through a data dimensionality reduction functionality via visibilitygridding, whereby the de-gridding and gridding operations underpinning † are explicitly encoded via the sparse holographic matrices H = G † 2 G . By doing so, the dimensionality of the measurement operator is effectively driven solely by the image size. The feature is enabled when the memory required to host the degridding matrices exceeds the available resources. Thirdly, further decomposition of each operator † into smaller blocks is enabled via a data-clustering step such that each data block corresponds to the aggregation of radially-neighbouring Fourier modes, identified under memory constraints. The number of CPU cores allocated for the forward step of the FB algorithmic structure is derived from the number of identified data clusters. These computing resources are initially used to compute the de-gridding matrices (G ) or the holographic matrices (H ), underpinning the operator † (only once), and are later used to host them and apply them at each FB iteration.

DATA, IMAGING, AND ANALYSIS
In this section, we provide a full description of the scrutinised data, including pre-processing and calibration steps. We provide the imaging settings of uSARA and WSClean, and outline the computing architecture and resources required to run both algorithms. Finally, we provide procedures for a coherent comparative analysis.

ASKAP Observations
ASKAP consists of 36 12-metre parabolic dish antennas, spanning six kilometres, at the Murchison Radio Observatory in Western Australia. ASKAP's original design includes Phased Array Feeds (PAFs; Hay et al. 2006) at the focus of each antenna, built as dual-polarisation chequerboard grids sensing signals in the frequency range of 700-1800 MHz. Signals received by each of the 188-element-sensor grids are cross-correlated to simultaneously form 36 separate primary beams (or pointings) on the sky (Hotan et al. 2014;McConnell et al. 2016). This PAF technology gives ASKAP an instantaneous field-ofview (FoV) of 30-square-degrees, making it the most rapid surveying radio telescope in the world (DeBoer et al. 2009). ASKAP's EMU Survey will survey the radio continuum over the entire southern sky (up to a northern declination of +30 • ) at a resolution of ∼ 10 arcsec with sensitivities reaching ∼ 10 Jy beam −1 , and is projected to detect more than 70 million galaxies. Science goals of the EMU collaboration include testing fundamental models for dark energy (e.g. Raccanelli et al. 2012), detecting the warm-hot intergalactic medium (e.g. Hodgson et al. 2020), tracing active galactic nuclei (AGN) and star formation up to high-redshifts (e.g. Mancuso et al. 2017), and mapping the radio continuum of the Galactic Plane and Centre (e.g. Riggi et al. 2021).
The ASKAP Early Science Broadband Survey (Project: AS034, PI: Lisa Harvey-Smith) began in 2018 with the aim to test observations using the full ASKAP array while covering scientifically interesting areas of sky, such as the Galaxy And Mass Assembly 23-hour field (GAMA 23). The EMU Pilot Survey (EMU-PS; Norris et al. 2021, Project: AS101, PI: Ray Norris), centred on RA, Dec: 21h00m00s, -55 • 00 00 , was carried out in 2019 covering an area of 270 deg 2 overlapping a portion of sky covered by the first data release of the Dark Energy Survey (DES; Abbott et al. 2018). From these early ASKAP data releases, we have selected three individual Scheduling Block (SB) beam observations covering extended, morphologically complex radio sources which represent robust test cases for image reconstruction with uSARA. ASKAP SBs contain 36 measurement sets each, corresponding to the 36 primary beam pointings. Therefore, to ensure maximum signal-to-noise, we have chosen the single beam measurement sets which were most closely centred on our primary targets of interest (beam 15 of SB8275, beam 12 of SB9351, and beam 35 of SB9442). Beam observations for each SB are carried out with the same specifications over a 10-hour total integration time. The observing band varies slightly between the selected Early Science data (SB8275 -central frequency 1013 MHz) and EMU-PS data (SB9351 and SB9442 -central frequency 943 MHz), yet both have instantaneous bandwidths of 288 MHz with physical channel intervals of 1 MHz. See Table 1 for further details on the selected observations.
The three selected SB-beam measurement sets, containing calibrated visibilities, were generated through the ASKAPsoft pipeline (Hotan et al. 2021). For direction-independent calibration, ASKAPsoft includes bandpass calibration using the standard calibrator PKS B1934-638 -observed for five minutes in each beam before or after the science target. Further calibration includes a cycle of phase-only self-calibration for every 1 MHz of data. Each beam observation in an SB is calibrated and imaged independently, and as a final step, ASKAPsoft stitches images together to form mosaics covering the full 30-square-degree FoV. For our imaging purposes, we used the ASKAP data products after ASKAPsoft processing, with no further flagging or calibration. Prior to imaging, we shifted the pointing centres of our selected beam observations by their beam offsets. Although most RI packages assume a Stokes parameter (intensity) of = ( + )/2 (where and represent the instrumental polarisations), ASKAPsoft uses the IAU (International Astronomical Union) definition = + . We did not apply a factor two correction for the IAU intensity convention in any of our final image reconstructions; therefore, flux densities in our images are halved when compared to the values found in ASKAPsoft mosaic images.  Table 2. Imaging settings of uSARA and WSClean. The notation ( ★ ) refers to the settings of the fields SB8275-15 and SB9351-12 and ( ★★ ) refers to the settings of the field SB9442-35.

Imaging Settings
To perform monochromatic sub-band imaging, we split the selected wide-band data into effective channels, or spectral windows (SPWs), across the full frequency band. More precisely, data from all three fields were binned uniformly into eight spectral windows with a bandwidth of 36 MHz each, under the assumption of nearly flat spectral behaviour within a spectral window. Data sizes per spectral window range from ∼ 0.8 to ∼ 1.2 GB. To further demonstrate the scalability of uSARA, we also imaged the third field (SB9442-35) over the full frequency band (covering 288 MHz) to form a single monochromatic image from ∼ 7.5 GB of data. We chose not to generate full-band monochromatic images for the other two fields (SB8275-15 and SB9351-12) since they host sources known to exhibit very steep spectral behaviour (specifically, the galaxy clusters Abell 3391 and SPT-CL 2023-5535; Brüggen et al. 2021;Hyeong-Han et al. 2020) and because both fields contain bright, large-scale artefacts.
According to the ASKAP Science Observation Guide 1 , recommended image size is set by the full width at half maximum (FWHM) of the primary beam. ASKAP beam pointings have a primary beam FoV approximated by a circular Gaussian with FWHM of 1.09 obs / (Hotan et al. 2021), where the observing wavelength, obs , is about 0.3 meters (at 1000 MHz) and is the dish diameter of a single ASKAP dish (12 m). We calculated the FWHM of a given 1 https://confluence.csiro.au/display/askapsst/ASKAP+ Survey+Science beam pointing to be ∼ 1.56 • at the middle of the frequency band. For the first two selected ASKAP fields (SB8275-15 and SB9351-12), we noted the presence of bright sources lying just outside of the primary beam that created sidelobes and decided to image a FoV covering twice the FWHM. For the third FoV (SB9442-35), we did not find any bright sources directly outside of the primary FWHM and chose to image a FoV covering 1.6 times the FWHM. In both cases, the imaged FoV is well beyond the FWHM of the primary beam. With a maximum baseline of 6268 m, and a native instrumental resolution between 9 and 12 arcsec over the bandwidth, we selected a cell size of 2.2 arcsec pixel −1 , corresponding to a super-resolution factor of 2 at the highest frequency. Under these considerations, the reconstructed images of fields SB8275-15 and SB9351-12 are 5500 × 5500 pixels in size, and those of the field SB9442-35 are 4096 × 4096 pixels. Unlike the output image of the CLEAN-based algorithm -by design restricted to the instrumental resolution through the application of the restoring beam to its estimated non-physical model image -the uSARA image retains super-resolution at the pixel size. On a final note, although the FWHM of the primary beam is used to determine the imaged FoVs, no primary beam correction is applied to the reconstructions in this work.
Systematic imaging was carried out on the three selected ASKAP fields with both uSARA and WSClean using the same imaging settings where applicable. A summary of the imaging settings is listed in Table 2. Our parallelised and automated imaging framework implemented in MATLAB, and the C++ software WSClean, were both run on Cirrus 2 , a UK Tier2 high-performance computing (HPC) service, comprising 280 compute nodes. Each node has 256 GB of memory and 36 CPU cores (with two hyperthreads each). Parameter choice in both algorithms is summarised in the following paragraphs.
WSClean parameters. The state-of-the-art WSClean imager corrects for the -effect via the -stacking approach (Offringa et al. 2014;Offringa & Smirnov 2017). Monochromatic images of each of the selected ASKAP fields were made with Multi-scale CLEAN, using a gain factor of 0.8. Auto-masking was enabled down to 2.5 times the standard deviation of the WSClean residual image, computed at the start of every major deconvolution iteration using the median absolute deviation estimator. The stopping criteria were set to a maximum iteration number of 10 6 and a 'cleaning' threshold equal to the estimated standard deviation of the residual image. Briggs robust weighting, with robust parameter −0.25 was considered in all the imaging experiments. These weights were stored in the IMAGING WEIGHT column of the measurement set tables and were extracted and utilised in uSARA imaging. The number of -stacks considered by WSClean is set automatically based on the theoretical bound derived in Offringa et al. (2014) and the available compute resources (see Sec. 6 for more details). For future reference, we note that CLEAN reconstructions are the so-called restored images, obtained as the sum of the non-physical model image convolved with the associated restoring beam, and the residual image. To support our flux density analysis, we also consider the WSClean model images convolved with the restoring beam, referred to as smoothed model images.
uSARA parameters. With the uSARA algorithm being implemented in MATLAB, data and associated information (including the standard deviation of the data domain noise) were extracted from the measurement set tables as a collection of MAT files using a dedi-   Tables 7-9, for details). In the reconstruction of the sub-band images of all three fields, the operator † was encoded via the underpinning sparse de-gridding matrices (G ). For the full-band monochromatic imaging experiment of the field SB9442-35, the resulting large data size triggered the dimensionality reduction feature. The operator † was therefore encoded via its underpinning holographic matrices (H ), reducing the memory requirements to host it by nearly a factor 5. The regularisation parameter was fixed to the heuristic value proposed in (7) for the full-band image of the field SB9442-35. However, some adjustment was found necessary to achieve a high reconstruction quality in terms of resolution and sensitivity for all sub-band images of the three selected fields, whereby was set from 0.7 to 0.8 times the heuristic. Higher values of the regularisation parameter resulted in somewhat smoother reconstructions and non-recovery of some fainter point sources, whereas lower values led to images highly contaminated by calibration artefacts. The applied adjusting factor may be partly attributed to the imperfection of the DIE calibration and the lack of DDE calibration, affecting the accuracy of the measurement operator. Nonetheless, it is generally consistent with the findings of the theoretical study of uSARA's heuristic in the context of simulated RI data (Terris et al. 2022). Finally, the stopping criteria of uSARA were set to their default values, including a maximum of 10 re-weighted minimisation tasks, and a relative variation between the image iterates of 0.0005.

Quantitative Analysis
Focusing on the target sources of the selected observations, we provide their flux measurements and in-band spectral index maps obtained by the two imaging algorithms. Estimated model images from uSARA are in units of [Jy pixel −1 ], whereas WSClean restored images are, by design, in units of [Jy beam −1 ]. For the sake of comparison, we normalised WSClean restored images by the area of the associated restoring Gaussian beam, denoted by beam and given by where MAJ and MIN are the respective major and minor axes of the restoring beam in pixel units (i.e. normalised by the cell size of 2.2 arcsec pixel −1 ). Diffuse emission of particular interest in our study presents a complex morphology with edges often blended into the background noise, as seen in WSClean maps. As is common practice, we measure the total flux of diffuse structure within manually generated regions which roughly follow ∼ 2 meas contours of the source, where meas is the root-mean-square (rms) noise measured in a nearby region void of sources from the WSClean restored map. Regions were handdrawn in the visualisation software SAOImageDS9 (Joye & Mandel 2003) to closely follow the contours of recovered signal in both the uSARA and WSClean maps, such that the same region was used when measuring flux density of the same diffuse source. Flux density measurements from uSARA images are expected to be lower than those measured from the WSClean restored images, due to the bias introduced from the WSClean residual map. For a more accurate comparison of flux density, we also provide measurements from WSClean smoothed model images. Note that error measurements on flux densities are not reported since the Early Science and Pilot Survey ASKAP observations are not yet validated against standard flux catalogues. All reported flux measurements and statistics were obtained using the SAOImageDS9 software.
Spectral index maps were created to showcase how sources of interest change over their morphology in electron energy distribution and consequently, in their spectral energy distribution. Firstly, the subband maps were smoothed via convolution with a 2D circular Gaussian kernel with axes of 5 arcsec for uSARA images and 20 arcsec for WSClean images. The spectral index maps were then obtained by fitting a first-order polynomial to the function log( ) = − log( ), where is the flux density for a given beam area at a given frequency and > 0 is the spectral index. Only the first six sub-bands are considered when generating these maps since the last two sub-bands were consistently found to recover less diffuse signal for primary targets of interest (possibly attributed to the steepness of their spectra).

RESULTS
In this section, we showcase high-resolution, high-fidelity images of our three selected fields produced by the uSARA algorithm and compare them to images made with multi-scale WSClean. Select images are presented in Figures 1 -6, showing the full imaged fieldsof-view of our chosen ASKAP observations, and include zoomedin views focusing on complex radio emission of interest and their associated optical images and spectral index maps. In Wilber et al. (2022), we provide FITS files of all spectral windows of the three selected fields imaged with both algorithms. For each field, images are also combined into animated GIF files to show how the recovered emission changes over the full frequency band. In what follows, we provide a detailed comparison of the morphology and flux density of specific sources between the uSARA and WSClean images.

First field: SB8275-15
Beam 15 of SB8275 covers a FoV containing the complex merging galaxy cluster system Abell 3391 -Abell 3395. This field has been recently observed with the eROSITA X-ray space telescope where a warm gas bridge has been discovered between the cluster pair as part of a 15 Mpc intergalactic filament (Reiprich et al. 2021). The field also contains multiple bent-tail and Fanaroff-Riley class I and II (FR-I & FR-II; Fanaroff & Riley 1974) radio galaxies, some belonging to the cluster system. A recent paper utilising mosaic images of SB8275 has confirmed more than 20 giant radio galaxies (at various redshifts) in the 30 deg 2 field (Brüggen et al. 2021).
In Figures 1 & 2, we present our images of the full FoV (3.36 • ) of the first sub-band (SPW:1) of SB8275-15, imaged with WSClean and uSARA, respectively. Both figures include zoomed-in views of the FR-I in Abell 3391 (a: top right panels), a FR-II cluster member in the east (b: middle right panels), and radio sources in Abell 3395 (c: bottom panels). The FR-I radio galaxies at the centre of Abell 3391 in the north and Abell 3395 in the south (see Table 1 for source names) are reconstructed with superior resolution by uSARA. This is most evident in the appearance of 'braiding' and gaps in the plasma of the jets, which are not resolved in the WSClean map. The FR-I in Abell 3391 is the brightest source in the field with a peak pixel flux of 20 mJy, as measured in the SPW:1 uSARA image, and 12 mJy as measured in the SPW:1 WSClean image. Calibration errors in this field manifest as strong ring-like artefacts emanating from these bright FR-I radio galaxies. Brüggen et al. (2021) successfully carried out additional direction-dependent calibration to reduce the effect of these large-scale artefacts, however, we note that we performed no such additional calibration prior to imaging these data and therefore the extended radial artefacts remain in our final images. Over the full frequency band, source morphology and the structure of artefacts change per spectral window (see associated GIFs in Wilber et al. 2022).

Abell 3395
The southern cluster Abell 3395 is made up of two sub-clusters, with separate X-ray peaks (Reiprich et al. 2021), indicating that this cluster is undergoing its own merger within the larger merging system. West of the central FR I in Abell 3395, there is a faint, diffuse source with a dim core connected by arms that extend to the north-west and south-west. As reported by Brüggen et al. (2021), the peak intensity of this dim core does not clearly coincide with a host galaxy visible in optical maps, and therefore the source can not be classified as a typical radio galaxy associated with a host AGN. The diffuse source is possibly a so-called radio 'phoenix' (see Kempner et al. 2004, for classification), re-ignited fossil plasma from past AGN activity which is no longer active. In the middle bottom panels of Figures 1 & 2, WSClean and uSARA images are overlaid as contours on an r-band optical image from DES DR1 (Abbott et al. 2018). As seen in these optical overlays, a cluster member galaxy sits just south of the dim radio core, at RA, Dec 06h25m09.95s, -54 • 30'34.4", raising the possible scenario that old AGN emission from this galaxy has drifted or has been disturbed and shifted by hot gas in the cluster environment, leading to the observed faint emission. The structure of this phoenix candidate appears more clearly defined in our uSARA image, while the edges are much more blended into background noise in the WSClean image. Most notably, the north-west and south-west limbs of the phoenix appear to pop out in the uSARA reconstruction, although they appear very faint and undefined in the WSClean map.
We measure the flux density of this candidate phoenix source using identical polygonal regions hand-drawn to closely follow the total recovered signal in each of the uSARA and WSClean sub-band images, such that the same region is used to measure between the two imaging algorithms. In the WSClean map of SPW:1, this polygonal region closely traces the 2 meas contour line, where meas = 2 Jy pixel −1 . Since the morphology of the phoenix source changes dramatically over the frequency band, a polygonal region was drawn for each spectral window. Interestingly, in the uSARA map of SPW:1 we find that the border of the phoenix's recovered signal can be traced by a ∼ 1 Jy pixel −1 contour level (see middle bottom panel of Figure 2), well below the estimated standard deviation of the noise in the image domain (see Eq. 6) calculated as 6 Jy pixel −1 . This finding remains consistent through subsequent spectral windows, indicating that the uSARA algorithm successfully recovers real signal below the estimated noise levels. In the WSClean maps, due to the blending of diffuse emission with background noise, the border of the phoenix is more clearly defined by contour lines at 2 to 3 meas .
For each sub-band, the measured flux densities from both uSARA and WSClean images are listed in Table 3. For the first two spectral windows, uSARA flux measurements of the candidate phoenix are greater; however, for all subsequent spectral windows, the WSClean flux measurements are consistently greater. This is likely due to the fact that the measured flux density of this region in the WSClean map is also integrated over the noise, which may increase for higher spectral windows and is amplified by artefacts emanating from the two bright FR I sources in the field. Indeed, lower flux densities are measured from the WSClean smoothed model images, bringing them closer to uSARA values, particularly at the lower end of the frequency band. For such a faint source, we can see how the flux measurement from WSClean restored images can be easily overestimated when mixed with a noisy background signal.
As apparent in Table 3, the flux of the phoenix source reconstructed by uSARA drops off dramatically as the frequency increases, indicating a steep spectral index ( > 1), as confirmed in Brüggen et al. (2021). However, this fading over the frequency band is less dramatic in our WSClean results, indicating that WSClean may be possibly biased when wide-band imaging is deactivated. Comparing the spectral index maps of the phoenix (shown in the bottom right panels of Figures 1 & 2), the general trend of steepening over the source morphology is similar. In the uSARA map, a steeper spectral index is seen in the phoenix's north-west limb. Interestingly, both spectral index maps show the phoenix hosting a steep core ( ∼ 1.5) surrounded by a halo of flatter emission ( < 1), in contrast to the total intensity. This steep core is more clearly defined in the uSARA spectral index map. The fact that the brightest portion of the core is steeper in its spectral index than the surrounding fainter emission provides further evidence that this source may be an AGN remnant or phoenix. The emission around the potentially dormant core may exhibit a flatter spectral index because it has undergone gentle re-energisation (e.g. de Gasperin et al. 2017) from turbulence or small-scale shocks in the intracluster medium. Likewise, gentle re-energisation may explain why the ultra-steep emission in the north-west and south-west limbs is visible only at the lower end of the band (i.e. subtle re-brightening of old AGN emission).

Second field: SB9351-12
Beam 12 of SB9351 covers a field containing the massive, merging galaxy cluster SPT-CL J2023-5535 (hereafter SPT2023) near the centre of the FoV and the X-shaped radio galaxy PKS 2014-55 on the western edge of the FoV. Since this beam observation lies on the western border of SB9351's full 30-square-degree field, we were unable to choose another beam observation that hosted the Xshaped radio galaxy closer to the pointing centre. Both the cluster and the radio galaxy have been recently studied: HyeongHan et al. (2020) announced the discovery of a radio halo and radio relic in the merging cluster SPT2023 using the same EMU-PS observation we have selected, and Cotton et al. (2020) used MeerKAT total intensity and polarisation observations to investigate the peculiar X-shaped morphology of PKS 2014-55.
In Figures 3 & 4, we present our images of the full FoV (3.36 • ) of the first sub-band (SPW:1) of SB9351-12, imaged with WSClean and uSARA, respectively. Both figures include zoomed-in views of the galaxy cluster SPT2023 (a: top right panels), a field of compact and point sources (b: middle right panels), and the X-shaped radio galaxy (c: bottom panels). The bright quasar RX J2024.3-5723 at the southern edge of the pointing introduces radial ring-type artefacts, which propagate up to 1 deg in the field.
In each of the zoomed-in views, uSARA shows higher resolution and more definition in the reconstruction of both compact and diffuse emission. However, the very faint emission of the radio halo in SPT2023 is not clearly recovered in the uSARA image. It is also apparent that some of the faintest point sources are missing from the uSARA image (see (b) panels of Figures 3 & 4). This loss of the faintest point sources is likely attributed to the choice of the uSARA regularisation parameter. A lower value would enable the recovery  of more of these point sources, but would also increase the amplitude of recovered calibration artefacts.

X-shaped Radio Galaxy
As apparent when comparing panels (c) of Figures 3 & 4, the X-shaped radio galaxy exhibits more clearly defined borders in our uSARA image. In the middle (c) panels of the same figures, WSClean and uSARA emission of the X-shaped radio galaxy are overlaid as contours on an r-band optical image from DES DR1 (Abbott et al.   Wilber et al. (2022) are provided all sub-band images combined into the GIF 'SB9351-12_uSARA', and the spectral index maps of the X-shaped radio galaxy obtained with uSARA and WSClean in the GIF 'SpectralIndexMap_PKS_2014_55', together with a colour blind-friendly version in the GIF 'SpectralIndexMap_PKS_2014_55_colorblind_friendly'. 2018). Again, we find that the border of the recovered uSARA signal of the X-shaped radio galaxy traces a contour level at ∼ 1 Jy pixel −1 , well below the estimated standard deviation of the image noise for this sub-band: = 6 Jy pixel −1 . In contrast, the diffuse edges of the east and west wings of the X-shaped radio galaxy blend into background noise in the WSClean map, such that the border is more clearly defined by 3 meas contour lines, where meas = 1.6 Jy pixel −1 . The total flux density of the radio galaxy is measured by summing the integrated flux density from three separate regions: the east wing, the core, and the west wing. The totalled flux density measurements from hand-drawn polygon regions for the lobes (roughly tracing emission bounded by the 2 meas contour line in WSClean images) and customised ellipse regions 3 for the core are listed in Table 4. The polygonal regions covering PKS 2014-55 were modified per each spectral window to more accurately follow the source morphology over the frequency band; however, identical regions were used to measure between the two imaging algorithms. As recorded in Table 5, the flux density falls off as the frequency increases, indicating a steepening of the spectral index for this source. Except for the first spectral window, we see again that the flux is consistently greater in the WSClean sub-band images, likely due to integration over the noise in the WSClean map. When measuring from WSClean smoothed model images, we find that the WSClean flux densities decrease, bringing them more in line with uSARA measurements.

X-Shaped RG
Spectral index maps constructed from WSClean and uSARA subband images are shown in the bottom right panels of Figures 3 & 4. The general trend of steepening and flattening is consistent between the two maps, with more patches of flatter emission occurring in the lower portion of the east wing. This flattening is indicative of turbulent "hot-spots", coinciding with brightening seen in the intensity maps. Our uSARA spectral index map shows a dramatic steepening on the edges of the wings, but this is likely to be an artificial steepening since the diffuse structure at the edges is not recovered as well at higher frequencies (see associated GIF available in Wilber et al. (2022), demonstrating how the source structure changes with the frequency).

SPT-CL J2023-5535
As shown in panel (a) of Figure 3, the radio halo in SPT2023 is barely recovered by WSClean, appearing as a very faint increase in the noise across the inner regions of the cluster. The SPT2023 radio relic is apparent in the WSClean map as a small, elongated arc at the western side of the cluster (we refer the reader to HyeongHan et al. 3 Since the WSClean image is convolved with the restoring beam, point sources are much more extended than they appear in our uSARA images. Therefore, we use differently sized ellipse regions to measure the flux density of the core of PKS 2014-55. (2020) for a more detailed analysis of these cluster sources). Our uSARA image does not recover diffuse emission resembling a radio halo (see panel (a) of Figure 4), likely due to the choice of regularisation parameter dampening signal below the estimated noise level; nonetheless, the radio relic appears more clearly defined and brighter than in the WSClean image. We note that the classification of the halo and relic from HyeongHan et al. (2020) was made using a full-band (288 MHz bandwidth) multi-frequency-synthesis image and that the signal in our individual sub-band images (36 MHz bandwidth) is much weaker. We report flux density measurements for the recovered relic source in SPT2023 in Table 5. Flux density measurements of the relic from WSClean smoothed model images are more in line with uSARA measurements, except for the last three sub-band where uSARA shows a decrease in flux.

Third field: SB9442-35
Beam 35 of SB9442 is centred on the complex radio source PKS 2130-538. This peculiar source, nicknamed "the dancing ghosts," is shaped by the jets and lobes of two radio galaxies in the Abell 3785 galaxy cluster. With the published catalogue and initial results of EMU-PS, Norris et al. (2021) included an analysis of the morphology of PKS 2130-538. Two AGN hosts contribute to the observed emission: a radio galaxy in the north at the centre of a bright filamentary arch, and a second radio galaxy in the south at the centre of a smaller arch on the eastern lobe. Despite the advantage of ASKAP's resolution -revealing previously unseen structure in PKS 2130-538 - Norris et al. (2021) point out that it is still unclear whether these two radio galaxies are superimposed or actually interacting with each other.
Similarly to the previous two fields, eight sub-band images were reconstructed and used to obtain the flux density measurements, and the first six sub-band images were used to generate spectral index maps (see associated GIF provided in Wilber et al. 2022). We also formed a single full-band monochromatic image, for increased sensitivity, and to demonstrate the scalability of our imaging framework. In Figures  5 & 6, we present our monochromatic images of the full FoV (3.36 • ) of SB9442-35, formed from the full-band data using WSClean and uSARA, respectively. The figures also include zoomed-in views of a region of background sources (a: top right panels), the star-forming galaxy NGC 7090 (c: mid-right panels), and "the dancing ghosts" (b: bottom panels). Unlike the two previous fields, the SB9442-35 images do not exhibit large-amplitude calibration artefacts. In fact, only a few compact radio galaxies, catalogued by the Sydney University Molonglo Sky Survey (SUMSS; Mauch et al. 2003), emanate localised ring-like artefacts, and do not hamper the recovery of our targets of interest.
Overall, there is a clear difference between WSClean and uSARA in terms of resolution. While uSARA recovers structure at higher resolution, diffuse components of the extended radio galaxy in panel (a) and the star-forming galaxy in panel (c) are also fully recovered . This monochromatic image is a uSARA model with a pixel resolution of 2.2 × 2.2 arcsec. Panels are the same as described in Figure 5. Rightmost (b) panel: spectral index map made with the first six sub-band images of uSARA after smoothing with a common circular Gaussian beam of 5 arcsec. In Wilber et al. (2022) are provided all sub-band images combined into the GIF 'SB9442-35_uSARA', and the spectral index maps of "the dancing ghosts" obtained with uSARA and WSClean in the GIF 'SpectralIndexMap_PKS_2130_538', together with a colour blind-friendly version in the GIF 'SpectralIndexMap_PKS_2130_538_colorblind_friendly'. to the sub-band map. As for uSARA images, a clear improvement in both resolution and sensitivity can be observed in the full-band image. The bridges, which consist of the jets from the northern and southern AGN, are more tightly collimated in the full-band image when compared to the sub-band image. Several faint point sources, missing from the sub-band image, also emerge. Interestingly, one can notice that the filamentary structure extending from the eastern lobe is more clearly defined and brighter in the uSARA reconstruction. This structure has a similar appearance to the synchrotron threads that were recently discovered in Ramatsoku et al. (2020); Dabbech et al. (2022). Recovery of this structure is an exciting result, as we may find more evidence of magnetic threads branching from and connecting extended lobes in radio galaxies with ultra-sensitive, super-resolved images.

Dancing Ghosts
The bottom right panels of Figures 5 & 6 show spectral index maps of the dancing ghosts, inferred from the sub-band images. The uSARA spectral index map contains much more detailed information on the spectra of the turbulent lobes of the northern AGN. We also observe that the higher intensity hot-spots exhibit flatter spectral indices, with steepening as the emission traces southward. The second AGN at the centre of the south-east bridge shows a flat core, as expected, and a second set of jets that bend back toward the northeast lobe of the first AGN. In WSClean maps, the second south-east lobe appears to blend back into the emission of the first north-east lobe, however, with uSARA we are able to see a distinct separation of these two portions, indicating that the emission may be somewhat superimposed when seen in projection. The steepness of the spectra in this region does not indicate a sort of re-energising from one lobe pushing onto the other, as it stays fairly consistent in a steep range between 1 < < 2 in both spectral index maps. The north-east thread exhibits a sharp, collimated spectral trend from steep to slightly flat to steep again -in contrast to the turbulent lobes -with an index ranging from 2.9 ≥ ≥ 0.9 and then 0.9 ≤ ≤ 3.2 following the thread from west-to-east in our uSARA map.
In Table 6, flux densities for the Dancing Ghosts are reported. Integrated flux density was measured from identical polygonal regions tracing the full recovered signal in both uSARA and WSClean maps, including the eastern and western lobes, the northern arch, the south-eastern jet, and the north-eastern filament. Unlike other diffuse sources of interest from previous fields, the uSARA flux densities are much more consistent with WSClean flux densities. This presents an interesting case that may be explained by i) the source having overall flatter spectral behaviour, and ii) the lower noise level due to the absence of large-scale calibration artefacts in this field. At higher spectral windows, uSARA flux densities are even slightly greater than WSClean, opposing the trend seen for the fainter, steeper diffuse sources of interest in previous fields. This may indicate that the faintest and steepest sources, with surface brightness near estimated noise levels, are more difficult to recover with uSARA, given our current choice of regularisation parameter. In both uSARA and WSClean full-band images, the borders of the dancing ghosts are clearly defined by a contour level at 3 Jy pixel −1 , which is 3 meas in the WSClean full-band map (where meas = 1 Jy pixel −1 ). The uSARA imaging framework calculated the estimated standard deviation of the noise in the image domain as = 8 Jy pixel −1 for the full-band data, however, both uSARA and WSClean images have captured signal below this value.

DISCUSSION
In this section, we discuss some specific points regarding the experiments performed with our novel automated and parallelised imaging framework, underpinned by the uSARA imaging algorithm (Terris et al. 2022;Dabbech et al. 2022). In comparison to the widely-used WSClean imager, our uSARA-ASKAP images host greater detail at high resolution, resolving never-before-seen structure at scales up to 2.2 arcsec.
Arguably the most exciting feature of our reconstructed uSARA maps is that they successfully capture both compact and diffuse radio emission. Standard CLEAN imaging methods often require a compromise on resolution in order to gain sensitivity to diffuse emission, necessitating multiple maps at various resolutions to accurately recover emission on separate spatial scales. In many scientific publications and published radio surveys, imaging results are often found to be separated into one high-resolution map and another low-resolution map, where point sources and compact emission may be subtracted (usually through a crude technique that can leave holes in the image). Here, we have demonstrated that uSARA reconstructs images with both high resolution and sensitivity to diffuse emission, enabling advanced scientific analyses of physically complex sources. Thanks to uSARA's super-resolution and superior sensitivity to diffuse and extended components, we are also able to generate highly-detailed spectral index maps, which aid in the classification of our targeted sources.
Moreover, we argue that uSARA can catapult the sensitivity and resolution of existing surveys. Unlike WSClean, no residual image is added to the final uSARA reconstruction. It is therefore highly likely that the flux densities of low surface-brightness sources measured from uSARA images are a closer approximation to the true source flux.
Calibration errors. Several of the ASKAPsoft mosaic continuum images from early science and pilot fields show that calibrated data are still affected by radial artefacts that propagate -up to several degrees in some cases -from bright radio sources. For the ASKAP data considered in this work, the largest source of DDEs (manifested as antenna-based image modulations) is most likely attributed to ASKAP's synthesised beams. ASKAP's beamforming hardware applies digitised weights to the 188-elemental receivers and sums them to form a single primary beam at 1 MHz intervals. Currently, ASKAP uses a common algorithm to maximise signal-to-noise ratio to calculate the weights for beam synthesizing (e.g. Jeffs et al. 2008;Ivashina et al. 2011). Holographic observations of ASKAP's beam patterns (see Hotan 2016, for details) show that their sensitivities vary over frequency from antenna to antenna. The complex-valued sensitivity pattern of the PAF beams therefore introduces DDEs which need to be modelled. Furthermore, antenna pointing errors can introduce direction-dependent antenna gains. ASKAPsoft corrects only for DIEs, and imperfections of the calibration undermine the accuracy of the measurement operator model. Consequently, reconstructed images can exhibit imaging artefacts and, more seriously, suffer from severely limited dynamic ranges. Examples of such artefacts can be seen in the field SB8275-15 containing the merging cluster system A3391-95, where large-amplitude ringing artefacts can be seen around the two bright FR-I radio galaxies (at the centres of the Abell clusters) in both uSARA and WSClean reconstructions. In spite of the lack of DDE calibration, overall uSARA exhibits higher reconstruction quality than CLEAN. On a further note, uSARA can be easily plugged as the imaging module into a joint calibration and imaging framework (Dabbech et al. 2021).
In-band spectral index maps. Spectral index maps obtained from sub-band monochromatic imaging with uSARA have shown to be more detailed than WSClean. All sub-band images having a common cell size at least two times beyond the observation's nominal resolution at the highest sub-band, the spectral index maps were inferred using a small blurring kernel, preserving their high level of detail in comparison with WSClean spectral index maps. We have found that the classification of some sources is more clearly defined based on the spectral behaviour exhibited in uSARA maps. Interestingly, some sources have shown steeper spectral indices in uSARA spectral index maps when compared to their WSClean counterparts, which can be the result of the increase of sensitivity brought by uSARA. However, this spectral trend warrants further investigation by moving to wide-band deconvolution for a more precise spectral analysis.

COMPUTATIONAL PERFORMANCE
To assess the computational performance of our imaging framework, we report specific information for all uSARA imaging experiments in Tables 7-9. Details of the measurement operator are listed in these tables, including the number of processed visibilities , the number of -stacks , the memory requirements to host its underlying sparse matrices H/G , and the computational cost in CPU core hour of † pre-computation C † . We also report the total compute time T Image and the computational cost in CPU core hour C Image of the deconvolution in the same tables. For comparison purposes, we report both T Image and C Image of WSClean runs in Table 10. For uS-ARA, the retained number of -stacks in each experiment from the planning step determines the sparsity of the de-gridding/holographic matrices underpinning the measurement operator and consequently the memory requirements to host † . The decomposition of the data and associated measurement operator into smaller blocks, and the resulting number of deployed CPU cores are inferred from the data clustering step (see Section 2.3). From Tables 7-9, one can notice that, in general, a larger number of visibilities will increase the computational cost of the measurement operator's pre-computation. Specific to FB iterations of uSARA, the compute time of the forward step is dominated by the Fourier transforms performed, while that of the backward step is driven by its sub-iterative nature. Both steps, being parallelised, are on par in terms of computing time, with the latter taking about 1.6 times longer on average.
Although the number of -stacks considered in WSClean is significantly more important, the reported computational time and cost in CPU core hour of uSARA is about 20 times higher on average. The superior computational performance of the standard WSClean RI imager is attributed to (i) its fast approximate data fidelity steps, whereby visibility gridding and de-gridding operations are only conducted few times, and (ii) its simplistic image model, in particular given the overall spatial compactness of the radio emission in the selected ASKAP fields, forcing multi-scale CLEAN to operate only on small-scales. However, the simple regularisation approach underpinning WSClean comes at the expense of lower imaging quality.
Even though WSClean is about one order of magnitude faster than our imaging algorithm in its current MATLAB implementation prototype, substantial improvement of uSARA's computational performance is expected when migrated to a production implementation using C++ and parallel libraries.   Table 8. SB9351-12 -uSARA: settings and computational cost of the reconstruction of the sub-band images. See the caption of Table 7 for the complete list of the acronyms.

CONCLUSIONS
This article presents a comprehensive study of a recently proposed monochromatic wide-field imaging framework for RI, validated on GB-sized, imperfectly calibrated data from Early Science and Pilot Survey ASKAP observations. The framework's underlying sparsitybased imaging algorithm, uSARA, yields image reconstructions with both super-resolution and high sensitivity to diffuse radio emission. Relying on a highly parallelised and automated implementation of the operators and functions involved, the imaging framework enables the formation of wide-field radio maps with large image dimensions.
Focusing on science cases characterised by complex morphology exhibiting both compact and faint diffuse emission, we have carried out the validation of the uSARA imaging algorithm through flux density measurements and spectral index maps generated from subband monochromatic images, in comparison with the widely-used WSClean imager. In spite of the large-scale artefacts due to imperfect calibration present in some of the selected fields, the uSARA images Table 9. SB9442-35 -uSARA: settings and computational cost of the reconstruction of sub-band and full-band images. See the caption of Table 7 for the complete list of the acronyms. In the full-band imaging experiment, the operator † is encoded via its underlying holographic matrices. We therefore report the memory occupied by these matrices, denoted by H [GB]. For all sub-band imaging experiments, we report the memory occupied by the de-gridding matrices denoted by G [GB].
show better reconstruction quality overall in comparison to WSClean, both in terms of resolution and sensitivity. Our uSARA-ASKAP images host more detailed structure of the targeted radio emissionmost clearly seen in the super-resolved jets and lobes of radio galaxies in each of the selected FoVs. In addition to high-resolution structure, faint diffuse emission has also been captured by uSARA, revealing more extended emission of intracluster radio sources which appeared blended into background noise in sub-band WSClean maps.
An advantageous result of our super-resolved uSARA-ASKAP images is the ability to generate more detailed spectral index maps. High-resolution structure, resembling turbulent emission in the radio lobes of several imaged radio galaxies, appears to closely trace small changes in the steepness of the observed spectra. Furthermore, each of our primary sources of interest exhibit steeper spectra in the uS-ARA spectral index maps, attributed to the increase in sensitivity and resolution delivered by the algorithm. Nonetheless, our spectral analysis of the target sources remains preliminary and warrants a deeper study using wide-band imaging. Planned upgrades to the uSARA framework -which will incorporate joint DDE calibration ) and wide-band deconvolution (Thouvenin et al. 2022a) -guarantee the robustness of future images, and consequently, more precise spectral information across all frequency channels.
In terms of scalability of the proposed imaging framework, we have demonstrated that its fully automated and parallel measurement operators enable image reconstruction of data up to 7.5 GB in size. Larger data dimensions and fields-of-view necessitate distributing the data and associated measurement operator into more blocks, and decomposing the image into more facets for the parallel application of the sparsity dictionary and its adjoint in the uSARA denoiser. By adding more computational resources, the time to solution can be maintained at a reasonable scale. In its current MATLAB implementation, the computational cost of our imaging framework remains higher than the benchmark imager WSClean. However, its migration to C++ leveraging parallel libraries, can substantially boost its computational efficiency, thus narrowing the computational gap with the state-of-the-art imager.
In the sequel to this series, Part II: "AIRI validated on ASKAP data," we investigate uSARA's sister algorithm AIRI, the second deconvolution algorithm built into our parallel automated imaging framework. AIRI differs from uSARA by exploiting learned Deep Neural Network (DNN) denoisers in lieu of uSARA's proximal operator in the underpinning FB algorithmic structure (5). AIRI was already very recently and briefly demonstrated on MeerKAT data (Dabbech et al. 2022). The algorithm will be validated on the same challenging ASKAP data as here, with the aim to further demonstrate its potential to deliver further imaging precision and faster reconstruction than uSARA.  Table 10. WSClean settings and computational cost of all sub-band and full-band images: for each experiment, the number of -stacks ( ), the computational cost of the deconvolution in CPU core hours (C Image ), and the imaging time in hours (T Image ) are reported. Each experiment used a single node on Cirrus, leveraging its 36 CPU cores. Since each CPU core has two threads, 72 virtual cores were detected, and consequently, the minimum number of -stacks was set automatically to 72.