# A Decentralized Framework for Radio-interferometric Image Reconstruction

S. Wang, S. Mignot, S. Prunet, L. Di Mascolo, M. Spinelli, A. Ferrari

► To cite this version:

S. Wang, S. Mignot, S. Prunet, L. Di Mascolo, M. Spinelli, et al., A Decentralized Framework for Radio-interferometric Image Reconstruction. The Astronomical Journal, 2025, 169, 10.3847/1538-3881/adc57b. insu-05099859

HAL Id: insu-05099859

https://insu.hal.science/insu-05099859v1

Submitted on 6 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Decentralized Framework for Radio-interferometric Image Reconstruction

S. Wang¹, S. Mignot¹, S. Prunet¹, L. Di Mascolo¹,², M. Spinelli¹, and A. Ferrari¹

¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, 06000 Nice, France; sunrise.wang@oca.eu

² Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD, Groningen, The Netherlands
Received 2025 January 15; revised 2025 March 20; accepted 2025 March 24; published 2025 May 2

## Abstract

The advent of large aperture arrays, such as the ones currently under construction for the SKA project, allows for observing the Universe in the radio-spectrum at unprecedented resolution and sensitivity. To process the enormous amounts of data produced by these telescopes, scalable software pipelines are required. This paper helps address this by proposing a framework that allows for decentralized radio-interferometric image reconstruction, parallelizing by spatial frequency. This is achieved by creating pseudo-full-resolution problems for each node by using the local visibilities together with previous major cycle reconstructed images from the other nodes. We apply the proposed framework to both multiscale CLEAN and sparsity regularized convex reconstruction and compare them to their serial counterparts across four different data sets of varying properties in the context of two visibility partitions. We found that the parallelization framework allows for significantly improved reconstruction times for images of similar quality. This was especially the case for our larger data sets where we were able to achieve close to the optimal 2× speedup.

*Unified Astronomy Thesaurus concepts: Radio interferometry (1346); Astronomy image processing (2306); Computational methods (1965); Distributed computing (1971)*

## 1. Introduction

The two telescopes currently under construction for The SKA project,<fnref n="3" /> SKA-Mid (350 MHz–15.35 GHz) in South Africa and SKA-Low (50–350 MHz) in Australia, will allow us to observe the Universe at these frequencies at unprecedented resolutions and sensitivities. These more detailed observations will result in large amounts of data to be processed. For SKA-Mid, projections estimate up to 2.375 TB/s from the dishes to the beamformer and correlator engines, and 1.125 TB/s from these to the imaging super-computer (G. P. Swart et al. 2022). For SKA-Low, the estimated data transfer is 0.725 TB/s from the antennas to the correlator and 0.29 TB/s from the correlator to the imaging supercomputer (M. G. Labate et al. 2022).

With these data rates, it is imperative for current image reconstruction algorithms to scale. On a coarse level, this involves parallelizing the sample data (i.e., visibilities) processing across multiple nodes of a cluster. This is typically performed along the embarrassingly parallel frequency and time domains, with independent reconstruction of each partition. Approaches that separate the image into facets for direction-dependent calibration, such as T. J. Cornwell & R. A. Perley (1992), R. J. Van Weeren et al. (2016), and C. Tasse et al. (2018), also enable parallelization in the spatial domain. This is more complicated due to the corner turn caused by each facet depending on all visibilities, and vice versa, but there has been recent work by P. Wortmann et al. (2024) that looks to address this.

Another possible axis to parallelize is by spatial frequency. This requires more careful treatment compared to parallelizing

by time and frequency, since image reconstruction methods typically need to process all the spatial frequencies together in order to obtain full-resolution reconstructions, leading to only a modicum of work being done in this respect (A. Onose et al. 2016; L. Pratley et al. 2019; S. Van der Tol et al. 2018; P. Wortmann et al. 2024). This paper addresses this by proposing a general decentralized framework to parallelize radio-interferometric image reconstruction by spatial frequency. This is achieved by filling in the missing information in each node with reconstructed images from the others, allowing each node to simultaneously reconstruct full-resolution images. It then combines these to obtain the final reconstructed image.

We present our framework within the traditional major-minor loop image reconstruction paradigm, and apply it to two different image reconstruction methods, one based on convex optimization with sparsity regularization, such as in Y. Wiaux et al. (2009), R. E. Carrillo et al. (2014), H. Garsden et al. (2015), A. Dabbech et al. (2015), and the other being multiscale CLEAN (T. J. Cornwell 2008). We evaluated our framework by comparing these parallel implementations with their serial counterparts across four different data sets. We found our framework to produce reconstructed images of comparable quality in substantially less time than the serial approaches, particularly in the case of larger data sets.

The remainder of this paper is structured as follows. We provide a brief overview of radio-interferometric imaging as well as a review of the literature in Section 2. We describe our proposed parallelization framework in Section 3. We discuss how it can be applied to two existing reconstruction methods in Section 4. Section 5 details our experimental setup, and Section 6 presents and discusses our results. Finally, we provide conclusions and discuss future avenues of research in Section 7.

<figure><img src="image_2.png" /></figure>

## 2. Radio-interferometric Imaging

Radio-interferometers measure the sky using arrays of antennas (i.e., aperture arrays). Pairs of antennas (i.e., baselines) produce

<footnote n="3">³ https://www.skao.int/</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** A high-level overview of the radio-interferometric pipeline. A reconstructed image $\hat{\tau}^n$ is compared to the measured visibilities $v$ in the de/gridding step, which outputs the residual $\tilde{\tau}^n$ in the spatial domain. This is passed to the estimation algorithm, which generates $\tilde{\tau}^n$, an estimate of $\tilde{\tau}^n$ without the effects of the measurement operator $F^\dagger G G^\dagger F$. The current reconstructed image $\hat{\tau}^n$ is then updated to $\hat{\tau}^{n+1}$ by adding $\tilde{\tau}^n$ to it. Since the estimation algorithm is usually iterative, the overall pipeline often has a major-minor loop structure.</figcaption></figure>

visibilities, the correlated instrumental response of the electrical field for some given time duration and electro-magnetic frequency. These can be related to the sky radiance distribution $I(l, m)$ using the measurement equation (O. M. Smirnov 2011):

$$V(u, v) = C_{uv} \iint_{D_{uv}(l, m)} \times \frac{I(l, m)}{n} e^{-2\pi i(lu + vm + w(n-1))} dl dm, \\ n = \sqrt{1 - l^2 - m^2} \qquad (1)$$

where $C$ represents the direction-independent effects, such as antenna gain, $D$ denotes the direction-dependent effects, such as phase gradients caused by the Earth's ionosphere, $(u, v, w)$ is the difference between antenna coordinates in the frame where $w$ is aligned with the phase center, and $(l, m)$ are the spatial angular coordinates on the celestial sphere, which is also the domain of the integral.

If we ignore the $D$ and $e^{-2\pi i w(n-1)}$ terms, Equation (1) simplifies to a two-dimensional Fourier transform, allowing us to retrieve an image of the sky emission through its inversion. This image contains artifacts due to the partial sampling caused by the nature of the antenna arrays and, especially in the case when imaging larger fields of view, by the omission of the $w$ and $D$ terms. Radio-interferometric imaging aims to correct for these, reconstructing an image that is usable for science.

The general approach of most radio-interferometric imaging algorithms is iterative and is illustrated in Figure 1. The residual $\tilde{\tau}^n$ between the measured visibilities $v$ and the current image estimate $\hat{\tau}^n$, also termed the dirty image when $n=0$, is computed in the de/gridding step. In image space, this can be expressed as

$$\tilde{\tau}^n = F^\dagger G (v - G^\dagger F \hat{\tau}^n), \qquad (2)$$

where $F$ and $F^\dagger$ are the fast Fourier transform (FFT) and its inverse, respectively, and $G$ and $G^\dagger$ are the gridding and degridding operators for resampling visibilities to regular and irregular positions, respectively. There is a myriad of methods for $G$ and $G^\dagger$, most of which also include the correction for $w$, with some examples being the work of T. J. Cornwell et al. (2008), A. R. Offringa et al. (2014), S. Van der Tol et al. (2018), and H. Ye et al. (2022).

The residual image $\tilde{\tau}^n$ is then sent to the estimation algorithm, which aims to remove the effect of the measurement operator $F^\dagger G G^\dagger F$ from the residual in a regularized manner,

producing $\tilde{\tau}^n$. Methods that aim to achieve this include CLEAN and its variants (J. A. Högbom 1974; T. J. Cornwell 2008; U. Rau & T. J. Cornwell 2011) that deconvolve $\tilde{\tau}^n$ in a greedy nonlinear manner much akin to matching pursuit (S. G. Mallat & Z. Zhifeng 1993), and convex optimization methods based on compressive sensing (Y. Wiaux et al. 2009; R. E. Carrillo et al. 2014; A. Dabbech et al. 2015). More recently, approaches that use artificial neural networks have also been used, either in an end-to-end manner (K. Schmidt et al. 2022; L. Connor et al. 2022), in plug-and-play (M. Terris et al. 2022), or unfolded (A. Aghabiglou et al. 2024) frameworks. After estimation, the reconstructed image $\hat{\tau}^n$ is updated to $\hat{\tau}^{n+1}$ with

$$\hat{\tau}^{n+1} = \hat{\tau}^n + \tilde{\tau}^n. \qquad (3)$$

The estimation step is usually iterative, therefore the imaging pipeline usually has a major-minor loop structure, with the loop shown in Figure 1 being the major and the estimation being the minor.

De/gridding is often the bottleneck of the imaging pipeline; for example, C. Tasse et al. (2018) showed that it took up to 94% of the total processing time for a serial implementation. Thus, there is much work that looks to expedite this step using parallelization. This can be done in a coarse or fine-grained manner. Fine-grained approaches parallelize on the local machine at a fine scale, such as per visibility or grid cell. There has been substantial amounts of work done in this area, both for the CPU as in A. H. Barnett et al. (2019) as well as the GPU as in J. W. Romein (2012), B. Merry (2016), and B. Veenboer et al. (2017).

Our method does not deal with fine-grained parallelization, and focuses primarily on coarser scales, aiming to parallelize the imaging pipeline across multiple nodes within a cluster. Typical methods for this are to parallelize according to the frequency and time axes. More recently, there has also been work that looks to distribute the gridding according to facets (N. Monnier et al. 2022) as well as by baseline (S. Van der Tol et al. 2018; L. Pratley et al. 2019) and regions of the uv-plane (A. Onose et al. 2016; L. Pratley et al. 2019). Different strategies can also be combined, such as in C. Gheller et al. (2023), which parallelizes by time and the v-axis, and the swiFTly framework (P. Wortmann et al. 2024), which allows for simultaneous distribution both by facets in the spatial domain and subgrids in the Fourier domain.

Although similar in aim, our work differs from most of the aforementioned methods in that it takes into account the entire imaging process, including the deconvolution step, as opposed to only dealing with distributed de/gridding. To this end, the most similar methods to ours that are mentioned above are those by A. Onose et al. (2016) and L. Pratley et al. (2019). These works employ primal-dual methods for the image reconstruction, and distribute and parallelize by regions of the uv-grid, or also by baseline length in the case of L. Pratley et al. (2019). Contrary to these approaches, which aim to solve directly for the measurement operator $F^\dagger G G^\dagger F$, our work focuses on proposing a general framework that can be applied to a variety of different methods, especially in the traditional major-minor loop paradigm where convolution by the point-spread function (PSF) is used as a surrogate for the measurement operator.

This provides several advantages over previous approaches. First, the general nature of our framework means that it is agnostic to the specific de/gridding and deconvolution algorithms used
within. Thus, our proposed framework can be used in conjunction
with many of the aforementioned methods. For example, for
de/gridding, we currently use the improved w-gridder (H. Ye
et al. 2022), but we could instead use the image-domain gridder
(S. Van der Tol et al. 2018) for a more fine-grained by-baseline
parallelism per node, or use the swiFTly framework (P. Wortmann
et al. 2024) if we are applying our framework to a facet-based
framework to take advantage of its distributed Fourier transform
properties. Similarly, for deconvolution, we study our framework
both with multiscale CLEAN (T. J. Cornwell 2008) and a sparsity
regularized convex optimization method akin to Y. Wiaux et al.
(2009), R. E. Carrillo et al. (2014), and H. Garsden et al. (2015),
but can also look to employ more recent machine learning
approaches (M. Terris et al. 2022; A. Aghabiglou et al. 2024).
Second, presenting our framework in the traditional major-minor
loop paradigm has performance benefits over approaches such as
A. Onose et al. (2016) and L. Pratley et al. (2019) because de/
gridding needs to only be performed every major cycle, rather than
every iteration of the optimization algorithm.

Finally, our framework is also decentralized, as opposed to
the more centralized structure of many of the previous works,
where the gridded visibilities need to be collected in a central
node for deconvolution. Although we do not explicitly
investigate this aspect of our framework in this paper, because
we only study the reduced two-partition case, it may have
performance implications when scaling to more partitions and
larger data sets, because it allows for more local communica-
tion between nodes.

Other works that bear resemblance to ours are those of
X. Cai et al. (2019) and S. Wang et al. (2024). X. Cai et al.
(2019) reconstruct the image progressively by incorporating the
missing information as previously reconstructed images. Their
idea is similar to ours, but it differs in aim because their goal is to
derive an on-line reconstruction method to alleviate memory
costs, whereas we look to improve scalability through
parallelization. To this end, they partition in time, whereas we
partition by spatial frequency. The work in this paper extends the
framework presented in S. Wang et al. (2024) by introducing
decentralized parallelism, since the previous framework only
allowed for the partitioning of visibilities by spatial frequency
but was still serial. This paper also provides multiscale CLEAN
parallelization results, whereas the focus was previously only on
sparsity regularized convex reconstruction.

**3. Decentralized Framework**

We propose a general framework for decentralized coarse-
grained parallelization by spatial frequency of the radio-
interferometric imaging pipeline. Each node reconstructs an
image using only the local visibilities during the first major
cycle. These images feature only partial spatial frequency
information. They are then sent to all other nodes from the
second major cycle onward, allowing all nodes to reconstruct
full-resolution images henceforth. The final reconstructed
image can be obtained by computing the weighted-average of
these full-resolution images. One way to view our framework is
that we perform de/gridding in parallel, while creating pseudo-
full-resolution problems to be solved by each node, using
previously reconstructed images from other nodes as a
surrogate for the missing visibilities. Note that this results in
a one major cycle delay between the residual image of the local
visibilities and the received images from other nodes, which

<figure><img src="image_4.png" /><figcaption>**Figure 2.** An illustration of our decentralized framework in the context of two partitions $V_L$ and $V_H$. We reconstruct two full-resolution images in parallel, achieved by sharing the filtered deconvolved residuals, $\Phi_H^n$ and $\Phi_L^n$, between different nodes after each major cycle, denoted as MC in the figure. The deconvolved residuals are filtered both for normalization purposes, as well as to ensure that only spatial frequencies from their corresponding domains are present. Low- and high-resolution deconvolved images are produced after the first major cycle, denoted as $\bar{\tau}_{V_L}^1$ and $\bar{\tau}_{V_H}^1$ respectively. Afterwards, full-resolution deconvolutions $\bar{\tau}_{V_L}^n$ and $\bar{\tau}_{V_H}^n$, where $n > 1$, are produced. The final reconstructed images of their respective nodes, $\hat{v}_{V_L}$ and $\hat{v}_{V_H}$, are then combined to produce the final reconstructed image $\hat{t}$.</figcaption></figure>

results in full-resolution images only being reconstructed after
two major cycles.

For this paper, we focus on presenting and studying our
framework within a reduced two-partition context, $V_L$ and $V_H$,
which contain the low and high spatial frequency visibilities,
respectively, and are illustrated in Figure 2. Here, after each
major cycle $n$, the nodes corresponding to the partitions $V_L$ and
$V_H$ send their latest filtered deconvolved residuals $\Phi_H^n$ and $\Phi_L^n$
(defined in Section 4), respectively, to the other node. These
images are then used as a surrogate for the missing visibility
information. The deconvolved residuals are filtered both for
normalization purposes and to ensure that only spatial
frequencies from their corresponding domains are present
because the remaining spatial frequency information is available
in the receiving node’s locally stored visibilities.

The partitioning of the visibilities is done according to their
spatial frequency, which we determine using their uv-distance
to zero, and each subset of visibilities is associated to a given
compute node. Figure 3 illustrates an example partition
configuration, where the low-resolution visibilities $V_L$ belong
to the orange domain $\mathcal{L}$, and the high-resolution visibilities $V_H$
belong to the green domain $\mathcal{H}$. We set $\mathcal{L} \cap \mathcal{H} \neq \emptyset$ because
de/gridding can cause frequency spillage after partitioning. We
define the partitioning configurations using variables $\ell$ and $\delta$,
with $\ell$ being the radius of a circle bisecting $\mathcal{L} \cap \mathcal{H}$, and $\delta$ being
the half-width of $\mathcal{L} \cap \mathcal{H}$.

4. Applying our Decentralized Framework to Existing
Reconstruction Methods

The proposed decentralized strategy is general and can be
applied to different major-minor loop reconstruction methods.
We apply our framework to two existing reconstruction
algorithms. One is based on convex optimization with sparsity
regularization (R. E. Carrillo et al. 2014; H. Garsden et al. 2015)
and the other is multiscale CLEAN (T. J. Cornwell 2008).

4.1. Sparsity Regularized Convex Reconstruction

The sparsity regularized convex reconstruction method we
apply within our framework produces a deconvolved residual
$\bar{\tau}^n$ for every major cycle $n$. It achieves this by using FISTA (A. Beck & M. Teboulle 2009) to solve:

$$ 
\begin{align} 
\alpha^n &= \underset{\alpha}{\operatorname{argmin}} \|\bar{\tau}^n - HW\alpha\|_2^2 + \lambda^n \|\alpha\|_1, \nonumber \\
\bar{\tau}^n &= W\alpha^n, \tag{4}
\end{align} 
$$

where $\bar{\tau}^n$ is the current residual image, $H$ the convolution operator by the PSF, $\lambda^n$ the current major cycle regularization parameter, and $W$ a wavelet transform operator. The final reconstructed image after $N$ major cycles is then $\sum_{i=1}^N \bar{\tau}^i$.

We previously studied this method in the context of a multistep reconstruction that partitioned by spatial frequency (S. Wang et al. 2024). To apply our decentralized framework to this method, we augment the objective function in a similar way, with an additional data-fidelity term incorporating the received images from the other node, from the second major cycle onward. To illustrate, the problem for the $V_L$ node becomes:

$$ 
\begin{align} 
\alpha_{V_L}^n &= \underset{\alpha}{\operatorname{argmin}} \|\Gamma_L(\bar{\tau}_L^n - H_L W \alpha)\|_2^2 \nonumber \\
&\quad + \gamma_n \|\rho_H^{n-1} - \Gamma_H W \alpha\|_2^2 + \lambda_{V_L}^n \|\alpha\|_1, \nonumber \\
\bar{\tau}_{V_L}^n &= W \alpha_{V_L}^n, \nonumber \\
\rho_H^{n-1} &= \sum_{j=1}^{n-1} \Phi_H^j - \Gamma_H \sum_{j=1}^{n-1} \bar{\tau}_{V_L}^j, \nonumber \\
\Phi_H^j &= \Gamma_H \bar{\tau}_{V_L}^j, \nonumber \\
\gamma_n &= 0 \text{ if } n=1 \text{ and } \gamma_n=1 \text{ otherwise,} \tag{5}
\end{align} 
$$

where the subscripts $\mathcal{L}$ and $\mathcal{H}$ denote operators and images restricted to the domains of $\mathcal{L}$ and $\mathcal{H}$, respectively, and the subscript $V_L$ refers to variables that are related to the node $V_L$, but not necessarily restricted to its domain (e.g., $\bar{\tau}_L^n$ contains only spatial frequencies from $\mathcal{L}$, whereas $\bar{\tau}_{V_L}^n$ is a full-resolution image produced by the $V_L$ node).

The first data-fidelity term evaluates against the local visibility information $V_L$ with PSF $H_L$ and is akin to the data-fidelity term in Equation (4) and the second data-fidelity term contains the received reconstructed images condensed into the term $\rho_H^{n-1}$, which subtracts the high spatial frequencies of the (already reconstructed) previous image estimate of the $V_L$ node $\Gamma_H \sum_{j=1}^{n-1} \bar{\tau}_{V_L}^j$ from the high spatial frequencies of the previous image estimate of the $V_H$ node, $\sum_{j=1}^{n-1} \Phi_H^j$. This ensures that the previous image estimate of the $V_L$ node gets subtracted from both fidelity terms. We use the filters $\Gamma_L$ and $\Gamma_H$ ensure that the data-fidelity terms evaluate only their respective spatial frequencies, are normalized according to their respective variance, and that the weights in $\mathcal{L} \cap \mathcal{H}$ sum to one.

$\rho_H^{n-1}$ can be viewed as a surrogate of Equation (2) for the high-resolution spatial frequencies, evaluated without the measurement equation $FG^\dagger GF^\dagger$ because there are no locally stored visibilities at these frequencies. The first term is replaced by $\sum_{j=1}^{n-1} \Phi_H^j$ because this is our current estimate of the sky at these spatial frequencies. The second term stays the same as Equation (2), where $\Gamma_H \sum_{j=1}^{n-1} \bar{\tau}_{V_L}^j$ is just the expanded version of the current filtered reconstructed image $\Gamma_H \bar{\tau}_{V_L}^j$.

It should be noted that we do not explicitly compute $\rho_H^{n-1}$, but rather update an aggregate image by adding the received and subtracting the locally deconvolved image every major cycle. This alleviates memory issues, particularly in the case of

<figure><img src="image_5.png" /><figcaption>Figure 3. The full set of visibilities is partitioned according to domains $\mathcal{L}$ (orange), which contains the low spatial frequency visibilities, and $\mathcal{H}$ (green), which contains the high spatial frequency visibilities. These are not mutually exclusive but rather have an overlap region, which is denoted in cyan, and is defined using $\ell$, which is the radius of a circle bisecting $\mathcal{L} \cap \mathcal{H}$, and $\delta$, the half-width of $\mathcal{L} \cap \mathcal{H}$.</figcaption></figure>

imaging large image cubes because we are not required to keep copies of all the reconstructed images. Additionally, we explicitly write out the terms for the computation of $\rho_H^{n-1}$ within the sum in Equation (5) to emphasize that the received images $\Phi_H^j$ are filtered by $\Gamma_H$ before transfer.

The $V_H$ node solves the same problem, with the only difference being that the visibility partition and received image spatial frequencies are swapped for the first and second data-fidelity terms, i.e., we use $\bar{\tau}_H^n$ and $H_H$ for the first data-fidelity term, and the residual image $\rho_L^{n-1} = \sum_{j=1}^{n-1} \Phi_L^j - \Gamma_L \sum_{j=1}^{n-1} \bar{\tau}_{V_L}^j$ for the second, where $\Phi_L^j = \Gamma_L \bar{\tau}_{V_L}^j$.

## 4.2. Multiscale CLEAN

The multiscale CLEAN algorithm iteratively finds the brightest source at the appropriate scale from the residual image, and then updates the residual by subtracting the found source convolved by the PSF. For the $n_{th}$ major cycle, the deconvolved image produced can be denoted as:

$$ \bar{\tau}^n = \text{ms-CLEAN}(\bar{\tau}^n, H, S, K), \tag{6} $$

where $\bar{\tau}^n$ is the residual image, $H$ is the PSF, $S$ the scales, and $K$ the number of CLEAN iterations (i.e., minor cycles).

Applying our decentralized framework to multiscale CLEAN, the problem for the $V_L$ node becomes:

$$ 
\begin{align} 
\bar{\tau}_{V_L}^n &= \text{ms-CLEAN}(\bar{\tau}_{V_L}^n \cup \mathcal{H}_L), H_{(\mathcal{L} \cup \mathcal{H}_L)}^n, \psi^n, K, \nonumber \\
\bar{\tau}_{(\mathcal{L} \cup \mathcal{H}_L)}^n &= \gamma^n \Gamma_L \bar{\tau}_L^n + \mu^n H_H \rho_H^{n-1}, \nonumber \\
H_{(\mathcal{L} \cup \mathcal{H}_L)}^n &= \gamma^n \Gamma_L H_L + \mu^n \Gamma_H H_H, \nonumber \\
\rho_H^{n-1} &= \sum_{j=1}^{n-1} \Phi_H^j - \Gamma_H \sum_{j=1}^{n-1} \bar{\tau}_{V_L}^j, \nonumber \\
\Phi_H^j &= \Gamma_H \bar{\tau}_{V_L}^j. \tag{7}
\end{align} 
$$
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dirty images of the four data sets we use for our experiments. The top two are simulated, whereas the bottom two are real observations performed with ALMA and the Very Large Array, respectively.</figcaption></figure>

Here, multiscale CLEAN deconvolves a dirty image $\tilde{t}_{(\mathcal{L}\mathcal{W}\mathcal{H})}^n$, which incorporates the high-resolution information with $H_H \rho_H^{n-1}$, using the PSF $H_{(\mathcal{L}\cup\mathcal{H})}^n$. $\mathcal{W}$ denotes an asymmetric union operator, since $\tilde{t}_{(\mathcal{L}\mathcal{W}\mathcal{H})}^n$ contains all the spatial frequencies in $\mathcal{H} \cup \mathcal{L}$, but $\tilde{t}_{(\mathcal{L}\mathcal{W}\mathcal{H})}^n \neq \tilde{t}_{(\mathcal{H}\mathcal{W}\mathcal{L})}^n$ since the local visibilities are different for each node. $\gamma^n$ and $\mu^n$ are equal to 1 and 0, respectively, for $n = 1$, and are otherwise weights corresponding to renormalizing the PSF and dirty images globally because these are typically initially normalized based on their respective local visibilities. $\psi^n$ is $S_\mathcal{L} \subseteq S$ for the $n = 1$, which includes only scales present in $\mathcal{L}$, and is otherwise equal to $S$. The other operators and variables are the same as the ones described for Equation (5). Likewise, the $V_\mathcal{H}$ node solves for a very similar problem, just with the various spatial frequency terms swapped, similar to Section 4.1.

## 5. Experimental Setup

We evaluate our framework across two simulated and two real calibrated data sets. This section serves to provide a brief overview of these and how they are partitioned, as well as specify algorithmic parameters and details on the hardware and software.

### 5.1. Data Sets and Partitioning

Our simulated data sets consist of cutouts of regions around Sgr B2 and Sgr C from a Meerkat survey observed at 1.28 GHz (I. Heywood et al. 2022). We simulate observations of these cutouts with pseudo R.A.–decl. coordinates, noise, and per-pixel angular resolutions, and use the telescope configurations of SKA1-Mid AA4 and SKA1-Low AA4, respectively. To create these data sets, we use RASCIL (T. J. Cornwell et al. 2020) to

<table><caption>**Table 1**
Data Set Description</caption><thead><tr><td>Data Set</td><td>No. Vis</td><td>Resolution</td><td>Tel.</td></tr></thead><tbody><tr><td>Sgr B2</td><td>9,361,440</td><td>512 × 512 @ 0.18</td><td>SKA-Mid</td></tr><tr><td>Sgr C</td><td>15,759,360</td><td>512 × 512 @ 0.429</td><td>SKA-Low</td></tr><tr><td>HL Tau</td><td>85,502,624</td><td>1500 × 1500 @ 0.005</td><td>ALMA</td></tr><tr><td>Cygnus A</td><td>82,246,592</td><td>1728 × 1728 @ 0.125</td><td>VLA</td></tr></tbody></table>

**Note.** The resolution is first given in number of pixels in the image along each dimension, and then the angular resolution per-pixel in arc-seconds. The full SKA1 AA4 telescope configurations were used for generating the simulated data.

<table><caption>**Table 2**
Number of Visibilities Per Partition for Test Data Sets</caption><thead><tr><td>Data Set</td><td>ℓ</td><td>V<sub>ℓ</sub></td><td>V<sub>κ</sub></td><td>V<sub>ℓ∩κ</sub></td></tr></thead><tbody><tr><td>Sgr B2</td><td>20</td><td>6.99M(.75)</td><td>2.53M(.27)</td><td>0.16M(.02)</td></tr><tr><td>Sgr C</td><td>25</td><td>9.05M(.57)</td><td>6.81M(.43)</td><td>0.11M(.01)</td></tr><tr><td>HL Tau</td><td>60</td><td>39.86M(.47)</td><td>46.24M(.54)</td><td>0.61M(.01)</td></tr><tr><td>Cygnus A</td><td>40</td><td>41.65M(.51)</td><td>41.83M(.51)</td><td>1.22M(.01)</td></tr></tbody></table>

**Note.** The number of visibilities are given in the millions, rounded off to the nearest 10 k. The fraction of the partition sizes to the total data set sizes are given in parentheses. The values for $ℓ$ are in pixels.

generate the visibility positions with the given observation parameters, and obtain the visibility values by performing degridding on the ground-truth images with the improved w-stacking gridder (H. Ye et al. 2022).

Our real data sets are the HL Tau<fnref n="4" /> band 6 observation from the Atacama Large Millimeter/submillimeter Array (ALMA) long-baselines survey and the first spectral window of the Cygnus A observation from M. L. L. Sebokolodi et al. (2020). Together, our data sets account for a variety of different source types, uv-coverages, fields of view, and numbers of visibilities. The real data sets also account for possible calibration errors and realistic noise properties. We image all data sets using uniform weighting. Table 1 provides a high-level overview of our data sets, with their dirty images shown in Figure 4.

We separate each data set into two partitions, looking to roughly balance the number of visibilities in each because this should result in minimal waiting time between the nodes during reconstruction, with the number of visibilities per partition for each data set shown in Table 2. We found the Sgr B2 data set to be particularly challenging to partition evenly. This is due to the high concentration of low-frequency visibilities present in the SKA1-Mid AA4 configuration. We discuss a possible solution to this issue in Section 6.1. We set $\delta = 1$ (pixel) when partitioning our data sets because we have previously shown in S. Wang et al. (2024) that the size of $\delta$ does not greatly change the final reconstruction quality, and it is best to keep $\mathcal{L} \cap \mathcal{H}$ small for efficiency.

### 5.2. Algorithmic Parameters

For the sparsity regularized convex reconstruction, we use for W a concatenation of the first eight Daubechies wavelets and select $\lambda_n$, $\lambda_{V_{\mathcal{L}_n}}$ and $\lambda_{V_{\mathcal{H}_n}}$ using preliminary experiments.

<footnote n="4">⁴ https://casaguides.nrao.edu/index.php?title=ALMA2014_LBC_SVDATA</footnote>
We found the general form of $\lambda_n = a\|\hat{\tau}^n\|_2 \times \text{pow}(b, n)$ to work well for both the serial and parallel approaches, with *a* and *b* being dependent on the data set. Generally, we found that it was necessary to increase $\lambda$ more aggressively in the $V_L$ node when reconstructing using the parallel approach. We use 50 FISTA iterations as the stopping condition for both the serial and parallel implementations for all data sets.

For multiscale CLEAN, the scales are $S = \{0, 1, 2, 4, 6, 10, 30\}$ for all our data sets because we found that keeping the scales small is necessary to obtain reconstructions of similar resolution to the sparsity regularized convex reconstruction. We set the maximum number of CLEAN iterations to $K=2000$ per major cycle, and set the absolute CLEAN threshold to values between 1e-4 and 2.5e-5 depending on the data set. Unlike the typical approach, we do not convolve our model with the CLEAN PSF when evaluating the reconstructions because this reduces the resolution of the final reconstruction, leading to poorer results.

For purposes of comparison, we terminate all methods after five major-cycle iterations. This includes the case of multiscale CLEAN, where the cleaning can terminate in fewer major iterations. The full-resolution reconstructions from each node are combined using equal weights: $\hat{c} = 0.5\hat{v}_{V_L} + 0.5\hat{v}_{V_R}$ because preliminary experiments using inverse variance weighting produced very similar results. Finally, we do not add the residuals back into the reconstructed images for any of the evaluated methods.

### 5.3. Implementation and Hardware

Our parallelization framework is implemented in Python, using RASCIL (T. J. Cornwell et al. 2020) for de/gridding, multiscale CLEAN, and general quality of life functionality such as reading from measurement sets, and mpi4py (L. Dalcín et al. 2005) for the decentralized computing framework. We implemented the sparsity regularized convex reconstruction method in Julia (J. Bezanson et al. 2017), and integrated it into our Python code. We parallelized the wavelet transforms of this implementation because each individual wavelet in the over-redundant dictionary can be transformed independently.

We use a three-node architecture for our framework, using one head and two worker nodes, with the worker nodes being responsible for the majority of the work and the head node responsible for aggregating statistics. In practice, this architecture can be replaced with a decentralized one with minimal change because our head node is mainly idle and primarily serves to combine the individual reconstructed images and record the elapsed time. Our code and results are all available on our on-line repository, which will be provided upon publication.

We ran our experiments on the Jean Zay supercomputer<fnref n="5" /> using the cpu_p1 partition, which contains nodes each with two Intel Cascade Lake 6248 CPUs and 192 GB of primary memory. Our data sets are loaded from the $WORK directory, which uses the performant parallel Lustre<fnref n="6" /> file system which manages a shared file system across multiple nodes. This should more closely reflect how disks are managed in the computer clusters of large radio interferometers, such as the SKA Science Data Processor.

## 6. Results and Discussion

We evaluate the serial implementations of both multiscale CLEAN and sparsity regularized convex reconstruction, and compare how they perform to their parallelized versions. For this, we primarily look to evaluate the reconstruction accuracy versus wall processing time. We also provide a more detailed breakdown and analysis of the division of computation time in order to provide a clearer picture of our framework’s bottlenecks and possible areas to improve. We also perform an experiment exploring how well our framework scales to larger images. For this, we provide a breakdown of computational times taken by the various subprocesses of our framework when imaging the Cygnus A data set at the same field of view but with a smaller angular pixel resolution, resulting in images of $10 \text{ k} \times 10 \text{ k}$ pixels, and compare these to corresponding processing times obtained from imaging with the original field size of $1728 \times 1728$ pixels.

### 6.1. Reconstruction Speed and Accuracy

We run all methods for five major-cycle iterations for each data set and plot the total processing time against reconstruction quality. Since we have ground truths for our simulated data sets, we evaluate the quality of their reconstructions using signal-to-noise ratio (S/N) in decibels (dB), which we define as $10 \log (\frac{\|i\|_2}{\|i-\hat{c}\|_2})$ where $i$ is the ground truth and $\hat{c}$ is the reconstructed image.

For the real data sets, we instead compare the residuals of the final reconstructed images against a noise reference image that contains just the noise of the measurements. If the image is reconstructed perfectly, the reconstruction will contain only signal and the residual only noise; thus, the distribution between corresponding pixels of the final residual and noise reference should be the same.

To obtain the noise reference image, we use a method based on jackknife resampling which has been previously applied to retrieve and study the noise and systematic properties of single-dish bolometric (A. Weiß et al. 2009; C. Romero et al. 2018, & references therein) and interferometric (L. Di Mascolo et al. 2023; J. Van Marrewijk et al. 2025, & references therein) measurements in both the radio and millimeter domains. For testing the difference between the distributions, we specifically opt to use $\|W(r, n)_5\|_2$, where $r$ is the residual image and $n$ is the noise reference image. $W(x, y)_5$ here is the Wasserstein-1 distance (C. Villani et al. 2009) which produces an image of the same resolution of both $r$ and $n$, containing the minimum distance to move the one-dimensional distribution of the $5 \times 5$ window surrounding each pixel from $r$ to $n$. We opt for this over a statistical test that produces p-values, such as the Mann-Whitney U test (H. B. Mann & D. R. Whitney 1947), because this differentiates between distributions that are far apart. We use the windowed approach because we only have one realization of the residual.

Figure 5 shows the results obtained for both our simulated and real data sets, with each major cycle demarcated with a circle, specifically just after the residual has been computed. The first data point of each plot corresponds to just after the initial dirty image was computed. In the interest of brevity, we refer to multiscale CLEAN and the sparsity regularized convex reconstruction methods as msc and L1, respectively, with their parallel versions having a p- prepended. This is done both for our plots and our discussions of the methods henceforth.

<footnote n="5">⁵ http://www.idris.fr/eng/jean-zay/jean-zay-presentation-eng.html</footnote>

<footnote n="6">⁶ https://www.lustre.org/</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time vs. accuracy for all four of our data sets, with each circle representing a major cycle. The left two simulated data sets use S/N in dB for accuracy, whereas the right two real data sets use the $L^2$ norm of the per-pixel Wasserstein distances between the noisy reference image and the current residual, computed within a 5 × 5 sliding window because we only have one realization of the residual.</figcaption></figure>

We found that for smaller data sets, such as our simulated ones, the parallel methods do not show a large improvement compared to their serial counterparts, whereas they do for our larger data sets. There are two main reasons for this. The first is that the cost of deconvolution is substantial compared to the cost of de/gridding for the smaller data sets, resulting in less overall gain in speed. The second is due to the uneven partitions, such as in the case of our Sgr B2 data set, which resulted in a lot of idling for the $V_H$ node. As mentioned in Section 5.1, we found it difficult to partition this data set equally due to the extremely high density of visibilities in the center of the uv-plane. A possible remedy is to employ a method, such as baseline-dependent averaging (S. J. Wijnholds et al. 2018), which will flatten the density distribution of visibilities in the uv-plane.

Overall, we found that the *p*-L1 and *L1* methods converge to images of similar quality. This is important because it shows that our framework does not introduce significant artifacts for these reconstruction methods, despite the problem being more challenging because the pseudo-full-resolution problems created for each node are less precise than if we were to reconstruct using all the visibilities in a single node. However, we did find *p*-L1 more difficult to regularize, with the $V_L$ and $V_H$ nodes sometimes requiring different regularization parameters. For multiscale CLEAN, we found that for the same number of major cycles, msc always reconstructed a slightly better image than the *p*-msc, which is largely expected due to the less precise nature of the pseudo-full-resolution problems. Despite this, we note that neither msc nor *p*-msc converged after our tested number of major cycles for most of our data sets because we did not explicitly tune the number of minor cycles for each data set; thus, it is not clear whether the final converged images will have similar quality. We did find that in terms of equal-time reconstructions, *p*-msc generally reconstructs images that are of higher quality compared to msc, as illustrated in Figure 5 for all our data sets except Sgr B2. Finally, we do note that for the HL Tau data set where msc and *p*-msc do converge, the final image qualities achieved by the two methods are similar. However, more experiments with a higher number of minor cycles are required to see if this result generalizes.

A drawback to our decentralized framework is that it requires two major-cycle iterations to generate a full-resolution reconstruction. This means that, at least in our demonstrated case of two partitions, there is no large difference in reconstruction speed if the corresponding serial method can reconstruct most of the information in the first few major cycles. This can be seen in some cases, such as in the Sgr C reconstruction, where the *p*-L1 and *L1* methods converge in

roughly the same amount of time, despite the faster major cycles in *p*-L1. This will become less of an issue with larger data sets, more partitions, and more complex reconstructions where more major-cycle iterations are required.

## 6.2. Division of Computation Time

We divide the different subprocesses into six categories. The first three are degridding, gridding, and deconvolution, which refers to their in-memory operations (i.e., not including disk I/O). Next, we have disk read, which is performed every major cycle because we ingest data sets piecemeal one channel at a time during de/gridding to account for data sets that are too large to store in memory. Finally we have idle time and other. Idle time refers to the time that individual nodes spend waiting for the other nodes to synchronize, whereas other encompasses processes such as converting polarizations, memory allocations, computing weights, etc. These processing times are wall times obtained using clock queries within our code.

Although our prototype implementation is largely unoptimized, we believe that this breakdown can nonetheless provide useful insight into the difference between the different methods and data sets because the same implementation is used for all our comparisons. Despite this, some aspects are bound to change with a more optimal implementation, e.g., we envisage the read portion to be drastically reduced.

Figure 6 shows how much time each subprocess takes across the different methods for each data set, with both $V_L$ and $V_H$ nodes being shown for the parallel methods. A result of note is that the reading takes up the lion's share of the de/gridding process, whereas the in-memory de/gridding often has comparable times to deconvolution. This is partially due to us reading all the visibilities from disk during each major cycle, but we noticed after profiling with pidstat (S. Godard 2025) that most of this processing time was not spent directly on disk I/O but rather on overhead conforming the read visibilities to RASCIL's format. This overhead is also the primary culprit contributing to the anomalous result of the parallel methods demonstrating more than the optimal 2× speedup when applied to the Cygnus A data set. A more efficient implementation should allow for computational speeds that are more comparable to current commonly used imaging software, such as WSClean (A. R. Offringa et al. 2014), and also allow us to obtain results for much larger data sets. We also note that the overall conclusions of our results should still hold in this case because both the de/gridding and the read subprocesses still increase with larger data sets.

Finally, we note that although the deconvolution times are less pronounced with our larger data sets, they are by no means negligible, as shown in Figure 7 for the *L1* method. This is
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The breakdown of computation time taken by the various subprocesses for our four methods.</figcaption></figure>

especially the case with our parallel framework, where the cost of de/gridding is reduced per node. Thus, in order to achieve better scaling, it is important to not only look at increasing the number of partitions, but also methods for improving the computational cost of the deconvolution method.

## 6.3. Scaling to Larger Images

Since our framework transfers images between nodes every major cycle, it is important to test how well it scales to imaging larger fields. For this, we image the Cygnus A data set described in Section 5 using *p*-msc and *p*-L1 at 10 k × 10 k pixels with a corresponding decrease in the pixel angular resolution, keeping the field of view the same. This results in the image sizes of 800 MB, as opposed to roughly 24 MB for the original test case of 1728 × 1728 pixels.

Table 3 shows a breakdown of the average wall processing times used by the various subprocesses for both the 10 k × 10 k and 1728 × 1728 configurations for a single major cycle, for both *p*-msc and *p*-L1. Overall, we found deconvolution to be the most expensive subprocess for both methods when increasing the pixel resolution. This was particularly the case with *p*-msc, where we found that RASCIL’s fully serial implementation of the method scales poorly to larger image sizes, leading to its cost dominating all other subprocesses. Additionally, we also found a significant increase in de/gridding’s processing time, although to a lesser extent when compared to deconvolution.

The increases in processing costs of both these subprocesses are roughly in line with the increase in number of pixels

<figure><img src="image_9.png" /><figcaption>**Figure 7.** A percentage breakdown of processing time used by the various subprocesses of the *L1* method across our four data sets.</figcaption></figure>

($\approx$33.5×) and are largely expected because the deconvolution methods operate in image space and the de/gridding algorithm has larger grids to convolve onto. This, in addition to the memory cost explosion, poses a problem if we wish to scale our framework to higher pixel resolutions, e.g., 100 k × 100 k, and necessitates that we investigate integrating our system with more scalable methods. For deconvolution, some examples are methods that take advantage of the sparsity in the domain of
<table><caption>**Table 3**

Cygurus A wall Computation Time Breakdown</caption><thead><tr><th>Alg.</th><th>Node</th><th>Pix. Res.</th><th>Deconv.</th><th>Degrid</th><th>Grid</th><th>Disk I/O</th><th>Transf.</th><th>Other</th></tr></thead><tbody><tr><td rowspan="4">p-msc</td><td>V<sub>L</sub></td><td>1728 × 1728</td><td>354.70 s</td><td>92.27 s</td><td>117.21 s</td><td>342.12 s</td><td>0.01 s</td><td>13.51 s</td></tr><tr><td>V<sub>H</sub></td><td>1728 × 1728</td><td>288.83 s</td><td>109.84 s</td><td>129.62 s</td><td>353.67 s</td><td>0.01 s</td><td>11.88 s</td></tr><tr><td>V<sub>L</sub></td><td>10 k × 10 k</td><td>17778.44 s (×50)</td><td>1450.58 s (×16)</td><td>2519.88 s (×21)</td><td>358.31 s (×1)</td><td>0.45 s (×52)</td><td>21.35 s (×2)</td></tr><tr><td>V<sub>H</sub></td><td>10 k × 10 k</td><td>18014.80 s (×62)</td><td>2159.66 s (×20)</td><td>2533.84 s (×20)</td><td>362.50 s (×1)</td><td>0.67 s (×51)</td><td>21.11 s (×2)</td></tr><tr><td rowspan="4">p-L1</td><td>V<sub>L</sub></td><td>1728 × 1728</td><td>91.66 s</td><td>92.44 s</td><td>111.52 s</td><td>332.31 s</td><td>0.02 s</td><td>11.52 s</td></tr><tr><td>V<sub>H</sub></td><td>1728 × 1728</td><td>89.33 s</td><td>108.17 s</td><td>126.85 s</td><td>359.43 s</td><td>0.02 s</td><td>13.49 s</td></tr><tr><td>V<sub>L</sub></td><td>10 k × 10 k</td><td>3595.75 s (×39)</td><td>1449.08 s (×16)</td><td>2457.80 s (×22)</td><td>365.53 s (×1)</td><td>0.59 s (×30)</td><td>20.16 s (×2)</td></tr><tr><td>V<sub>H</sub></td><td>10 k × 10 k</td><td>3573.73 s (×40)</td><td>2173.30 s (×20)</td><td>2555.44 s (×20)</td><td>363.62 s (×1)</td><td>0.60 s (×25)</td><td>20.22 s (×1)</td></tr></tbody></table>

**Note.** Wall computation times are given per major cycle, averaged across four major cycles (the first major cycle was omitted due to deconvolution times being very different from the rest due to only using a subset of the scales). The other category includes smaller tasks such as subtracting visibilities, adding images, memory allocation, applying weights, and converting polarization frames. The increases in time for the respective nodes between the different field sizes are provided in parentheses, rounded off to the nearest integer.

recovery, such as PolyCLEAN (A. Jarret et al. 2025), and methods that separate the problem into subimages, such as the DDFacet’s SubSpace Deconvolution (C. Tasse et al. 2018), Faceted HyperSARA (P.-A. Thouvenin et al. 2023), or the parallel deconvolution of WSClean (S. J. Wijnholds et al. 2023, & references therein). As for de/gridding, we could look to reduce the grid size by splitting the de/gridding into a sum of smaller grids which we obtain by modulating the visibilities toward the center of the uv-plane, similar to the strategy used in image-domain gridding (S. Van der Tol et al. 2018).

We also found that the inter-node image transfer exhibited a significant increase in processing times, on a similar scale to deconvolution. However, it is currently not the bottleneck because its total processing time is negligible compared to the other subprocesses, taking slightly over half a second for the 10 k × 10 k case. Since this operation is relatively straightforward, we can also predict its processing times when scaling to even larger images. For example, at the transfer rates shown in Table 3, which equates to roughly 1.3 GB/s because the images are 800 MB, a 100 k × 100 k image which has a size of 80 GB will take around 65 s to transfer, which is less than the cost of deconvolution of even the 1728 × 1728 case.

Although not an issue for our studied two-partition case, the transfer costs may become problematic once we extend our framework to more partitions. This is because a naive implementation would require each node to share its deconvolved images with all other nodes every major cycle, resulting in transfer data amounts that grow quadratically with regard to the number of partitions. We envisage several methods that can mitigate this. One is to employ a lossless compression scheme to the images before transfer. This should be achievable because the transferred images are band limited due to the supports of the filters $G_L$ and $G_H$ being finite in the Fourier domain. This also scales well to more partitions as the supports of the visibilities decrease with more partitions, and will be nearly mutually exclusive assuming a relatively small $\delta$ is used when partitioning the data sets. This would lead to near linear scaling, as opposed to quadratic. Another possible method that can mitigate transfer costs is to have more local communication, where nodes only communicate with a subset of the other nodes.

## 7. Conclusion and Future Work

This paper proposes a decentralized radio-interferometric imaging framework that allows parallelizing by spatial frequency

across two partitions. This additional axis of parallelization allows for more flexibility when scaling image reconstruction for future radio interferometers, such as the ones currently under construction for the SKA project. We applied our framework to two existing image reconstruction methods, and evaluated it against the serial implementations across four different data sets with various properties. We found that our framework enables faster reconstruction of images of similar quality compared to the serial implementations, and is particularly advantageous in the case of large data sets, such as the ones that will be produced by the SKAO telescopes, as the relative cost of deconvolution becomes less important.

The main drawback to our proposed method is that it requires at least two major cycles to converge. Thus, the performance advantages that it offers are less obvious if its serial counterpart converges in just a few major cycles. We believe this will not be an issue because the aim is to scale to very large, complex data sets that may require many major cycles to converge, and be separated into many partitions.

There are several avenues in which we are looking to expand our work: (a) investigate our framework using more partitions, on larger data sets, which requires both employing a compression method to the transferred images and defining a more intricate communication scheme in order to optimize communication costs, and a more intricate filter design, and is important to test both the scalability and the decentralized nature of our framework; (b) apply our framework to state-of-the-art image reconstruction methods, both to test how well it generalizes and to reconstruct larger fields in a scalable manner; (c) reimplement our prototype to be more performant, which will allow us to achieve more meaningful comparisons with other software and benchmark resource usage, such as CPU, memory, disk and network I/O, for larger data sets; (d) investigate the use of baseline-dependent averaging (S. J. Wijnholds et al. 2018) or a similar method with our framework, which should reduce the amount of idle time when using telescope configurations that have very high densities in the center of the uv-plane.

## Acknowledgments

As part of the “France 2030” initiative, this work has benefited from a State grant managed by the French national research agency (Agence Nationale de la Recherche) attributed to the Exa-DoST project, and bearing the reference : ANR-22-EXNU-0004. This work was also financed by DARK-ERA
(ANR-20-CE46-0001-01), UCAJEDI (ANR-15-IDEX-01), and
was granted access to the HPC resources of IDRIS under the
allocation AD010415840 made by GENCI.

ORCID iDs

S. Wang https://orcid.org/0000-0002-5038-9531

S. Mignot https://orcid.org/0000-0002-8515-939X

S. Prunet https://orcid.org/0000-0002-1755-4582

L. Di Mascolo https://orcid.org/0000-0003-3586-4485

M. Spinelli https://orcid.org/0000-0003-0148-3254

A. Ferrari https://orcid.org/0000-0002-7297-7625

References

Aghabiglou, A., Chu, C., Dabbech, A., & Wiaux, Y. 2024, ApJS, 273, 3

Barnett, A. H., Magland, J., & Af Klinteberg, L. 2019, SJSC, 41, C479

Beck, A., & Teboulle, M. 2009, SIIMS, 2, 183

Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAMR, 59, 65

Cai, X., Pratley, L., & McEwen, J. D. 2019, NMRAS, 485, 4559

Carrillo, R. E., McEwen, J. D., & Wiaux, Y. 2014, MNRAS, 439, 3591

Connor, L., Bouman, K. L., Ravi, V., & Hallinan, G. 2022, NMRAS, 514, 2614

Cornwell, T. J. 2008, ISTSP, 2, 793

Cornwell, T. J., Golap, K., & Bhatnagar, S. 2008, ISTSP, 2, 647

Cornwell, T. J., & Perley, R. A. 1992, A&A, 261, 353

Cornwell, T. J., Wortmann, P., Nikolic, B., Wang, F., & Stolyarov, V. 2025, Radio Astronomy Simulation, Calibration and Imaging Library (RASCIL), https://developer.skao.int/projects/rascil/en/latest/

Dabbech, A., Ferrari, C., Mary, D., et al. 2015, A&A, 576, A7

Dalcín, L., Paz, R., & Storti, M. 2005, JPDC, 65, 1108

Di Mascolo, L., Saro, A., Mroczkowski, T., et al. 2023, Natur, 615, 809

Garsden, H., Girard, J. N., Starck, J.-L., et al. 2015, A&A, 575, A90

Gheller, C., Taffoni, G., & Goz, D. 2023, RASTI, 2, 91

Godard, S. 2025, System performance tools for the Linux operating system (sysstat), https://github.com/sysstat/sysstat

Heywood, I., Rammala, I., Camilo, F., et al. 2022, ApJ, 925, 165

Högbom, J. A. 1974, A&AS, 15, 417

Jarret, A., Kashani, S., Rué-Queralt, J., et al. 2025, A&A, 693, A225

Labate, M. G., Waterson, M., Alachkar, B., et al. 2022, JATIS, 8, 011024

Mallat, S. G., & Zhifeng, Z. 1993, ITSP, 41, 3397

Mann, H. B., & Whitney, D. R. 1947, Ann. Math. Stat., 18, 50

Merry, B. 2016, A&C, 16, 140

Monnier, N., Guibert, D., Tasse, C., et al. 2022, in Proc. of the 2022 IEEE Workshop on Signal Processing Systems, ed. J. McAllister & M. Pelcat (Piscataway, NJ: IEEE), 1

Offringa, A. R., McKinley, B., Hurley-Walker, N., et al. 2014, MNRAS, 444, 606

Onose, A., Carrillo, R. E., Repetti, A., et al. 2016, MNRAS, 462, 4314

Pratley, L., McEwen, J. D., Avezac, M., et al. 2019, arXiv:1903.04502

Rau, U. & Cornwell, T. J. 2011, A&A, 532, A71

Romein, J. W. 2012, in Proc. of the 26th ACM Int. Conf. on Supercomputing, ed. U. Banerjee & K.A. Gallivan (New York: ACM), 321

Romero, C., McWilliam, M., Macias-Pírez, J.-F., et al. 2018, A&A, 612, A39

Schmidt, K., Geyer, F., Fröse, S., et al. 2022, A&A, 664, A134

Sebokolodi, M. L. L., Perley, R., Eilek, J., et al. 2020, ApJ, 903, 36

Smirnov, O. M. 2011, A&A, 527, A106

Swart, G. P., Dewdney, P., & Cremonini, A. 2022, JATIS, 8, 011021

Tasse, C., Hugo, B., Mirmont, M., et al. 2018, A&A, 611, A87

Terris, M., Dabbech, A., Tang, C., & Wiaux, Y. 2022, MNRAS, 518, 604

Thouvenin, P.-A., Abdulaziz, A., Dabbech, A., Repetti, A., & Wiaux, Y. 2023, MNRAS, 521, 1

Van der Tol, S., Veenboer, B., & Offringa, A. R. 2018, A&A, 616, A27

Van Marrewijk, J., Kaasinen, M., Popping, G., et al. 2025, A&A, 695, A204

Van Weeren, R. J., Williams, W. L., Hardcastle, M. J., et al. 2016, ApJS, 223, 2

Veenboer, B., Petschow, M., & Romein, J. W. 2017, in IEEE Int. Parallel and Distributed Proc. Symp., ed. L. O'Conner (New York: IEEE), 545

Villani, C., et al. 2009, Optimal Transport: Old and New (Berlin: Springer)

Wang, S., Prunet, S., Mignot, S., & Ferrari, A. 2024, A&A, 692, A61

Weiß, A., Kovács, A., Coppin, K., et al. 2009, ApJ, 707, 1201

Wiaux, Y., Jacques, L., Puy, G., Scaife, A. M. M., & Vandergheynst, P. 2009, MNRAS, 395, 1733

Wijnholds, S. J., Willis, A. G., & Salvini, S. 2018, MNRAS, 476, 2029

Wijnholds, S. J., Dijkema, T. J., Groot, H., et al. 2023, in XXXVth General Assembly and Scientific Symp. of the Int. Union of Radio Science (URSI GASS), ed. S. J. Wijnholds & Y. Gupta (New York: IEEE), 1

Wortmann, P., Kent, J., & Nikolic, B. 2024, A&A, 687, A152

Ye, H., Gull, S. F., Tan, S. M., & Nikolic, B. 2022, MNRAS, 510, 4110