Pierre-Antoine Thouvenin

Arwa Dabbech

Ming Jiang

Abdullah Abdulaziz

Jean-Philippe Thiran

Adrian Jackson

Yves Wiaux

Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): II. Code and real data proof of concept

Keywords: techniques, image processing, techniques, interferometric

In a companion paper, a faceted wideband imaging technique for radio interferometry, dubbed Faceted HyperSARA, has been introduced and validated on synthetic data. Building on the recent HyperSARA approach, Faceted HyperSARA leverages the splitting functionality inherent to the underlying primal-dual forward-backward algorithm to decompose the image reconstruction over multiple spatio-spectral facets. The approach allows complex regularization to be injected into the imaging process while providing additional parallelization flexibility compared to HyperSARA. The present paper introduces new algorithm functionalities to address real datasets, implemented as part of a fully fledged MATLAB imaging library made available on Github. A large scale proof-of-concept is proposed to validate Faceted HyperSARA in a new data and parameter scale regime, compared to the state-of-the-art. The reconstruction of a 15 GB wideband image of Cyg A from 7.4 GB of VLA data is considered, utilizing 1440 CPU cores on a HPC system for about 9 hours. The conducted experiments illustrate the reconstruction performance of the proposed approach on real data, exploiting new functionalities to leverage known directiondependent effects (DDEs), for an accurate model of the measurement operator, and an effective noise level accounting for imperfect calibration. They also demonstrate that, when combined with a further dimensionality reduction functionality, Faceted HyperSARA enables the recovery of a 3.6 GB image of Cyg A from the same data using only 91 CPU cores for 39 hours. In this setting, the proposed approach is shown to provide a superior reconstruction quality compared to the state-of-the-art wideband CLEAN-based algorithm of the WSClean software.

INTRODUCTION

The Faceted HyperSARA approach [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF]) is a novel wideband imaging technique for radio interferometry, enabling a fully distributed reconstruction over the wideband data and image via data blocking and spatial and spectral faceting. The approach thus paves the way towards scalability to the large data and image sizes expected from modern radio telescopes. This paper presents a new MATLAB library implementing the Faceted HyperSARA approach for real radio interferometric (RI) data with a focus on (i) its structure for an enhanced user experience and (ii) the added functionalities for a fast and accurate RI measurement operator in the presence of unknown noise levels due to calibration errors. More specifically, models of the RI measurement operator can encapsulate the 𝑤-term, a direction-dependent effect resulting from the non-coplanarity of the radio array, via 𝑤-projection (Cornwell 2008;[START_REF] Wolz | [END_REF]Dabbech et al. 2017). Available direction-dependent effect solutions of other origins (e.g., atmospheric or instrumental) obtained from a ★ E-mail: y.wiaux@hw.ac.uk calibration pre-processing step (Repetti et al. 2017;[START_REF] Dabbech | [END_REF], can also be integrated into the measurement operator. Furthermore, unknown effective noise levels induced by the imperfect calibration are estimated using an adaptive strategy inspired from Dabbech et al. (2018). Finally, large data volumes can be accommodated via a visibility gridding-based data dimensionality reduction scheme (Kartik et al. 2017).

Two proof-of-concept experiments are proposed to assess the algorithm's performance on real wideband RI data acquired with the Very Large Array (VLA). The first experiment is conducted to evaluate the reconstruction performance of Faceted HyperSARA in terms of imaging quality through the reconstruction of a 15 GB image cube of Cyg A from 7.4 GB visibilities across 480 channels, spanning the frequency range 3.979-8.019 GHz . Utilizing the same data, a second experiment is performed with Faceted HyperSARA, combined with a joint image and a visibility gridding-based data dimensionality reduction technique to accommodate larger data volumes. More precisely, given the moderate spectral variation of Cyg A in the frequency range of interest, we target a 16-fold reduction of the spectral resolution. The visibility gridding offers a 7-fold reduction in the data volume.

The article is organized as follows. In Section 2, the faceted prior model and objective function underlying the Faceted HyperSARA approach are summarized. Newly added functionalities for a more robust and scalable measurement operator are detailed in Section 3. Section 4 introduces a fully fledged MATLAB library, addressing the RI imaging problem with algorithms from the SARA family. These are the monochromatic SARA approach (Carrillo et al. 2012;Onose et al. 2017) and the wideband approaches HyperSARA (Abdulaziz et al. 2016(Abdulaziz et al. , 2019b) ) and Faceted HyperSARA. To study the performance of the latter approach, a first experiment is conducted utilizing VLA data in Section 6, assessing its imaging quality and computing time in comparison with the SARA approach (Onose et al. 2017). Section 7 illustrates the potential of combining Faceted HyperSARA and a dimensionality reduction technique to handle larger data and image sizes, utilizing the same data as in Section 6. The approach is assessed in comparison with the state-of-the-art wideband CLEANbased method of the WSClean software (Offringa & Smirnov 2017). Finally, conclusions and perspectives are reported in Section 8.

FACETED HYPERSARA: A SUMMARY

The problem and notations used to formulate the RI imaging problem are recalled in this section (notation summary provided in Table 1). A short summary of the Faceted HyperSARA approach, introduced in the companion paper [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], Sections 2 and 3), is also provided.

Problem formulation

Wideband RI imaging refers to the reconstruction process of unknown images of the sky over 𝐿 frequency channels from noisy data -referred to as visibilities -acquired by an array of antennas. Assuming the image to be reconstructed covers a small field of view on the celestial sphere with perfect antenna calibration, each visibility point, acquired by a pair of antennas at a given frequency, corresponds to a spatial Fourier component of the sky surface brightness. The associated spatial frequency (in units of the observation's wavelength) corresponds to the projection of the corresponding baseline onto the plane perpendicular to the line of sight [START_REF] Thompson | Interferometry and Synthesis in Radio Astronomy[END_REF]). These provide an incomplete coverage of the 2D spatial Fourier plane (also called 𝑢𝑣-coverage) of the image. The measurement process can be modeled for each frequency channel 𝑙 ∈ {1, . . . , 𝐿} as (Abdulaziz et al. 2016(Abdulaziz et al. , 2019b)

y 𝑙 = 𝚽 𝑙 x 𝑙 + n 𝑙 , with 𝚽 𝑙 = 𝚯 𝑙 G 𝑙 FZ, (1)
where y 𝑙 ∈ C 𝑀 𝑙 is the vector of the whitened visibilities (a process known as natural weighting in RI imaging) acquired at the frequency channel indexed by 𝑙 ∈ {1, . . . , 𝐿}, x 𝑙 ∈ R 𝑁 + denotes the underlying unknown image, and n 𝑙 ∈ C 𝑀 𝑙 represents the measurement noise modeled as a realization of a standard complex independent and identically distributed random variable. The measurement operator 𝚽 𝑙 represents a non-uniform fast Fourier transform (NUFFT) (Fessler & Sutton 2003), composed of a zero-padding and scaling operator Z ∈ R 𝐾 × 𝑁 (with 𝐾 the number of points after zero-padding), the 2D Discrete Fourier Transform represented by the matrix F ∈ C 𝐾 ×𝐾 , and an interpolation matrix G 𝑙 ∈ C 𝑀 𝑙 ×𝐾 whose rows contain a compact-support interpolation kernel centered at the corresponding 𝑢𝑣-point. Note that available models of direction-dependent effects can be incorporated in the matrix G as convolution with their compact-support Fourier representation (Repetti et al. 2017;[START_REF] Dabbech | [END_REF]. The diagonal noise whitening matrix 𝚯 𝑙 contains the inverse of the noise standard deviation associated with the raw measurements.

To efficiently address the large number of visibilities in the reconstruction algorithm, the visibility vectors y 𝑙 can be further decomposed into 𝐵 blocks, y 𝑙 = (y 𝑙,𝑏 ) 1≤𝑏≤ 𝐵 (Onose et al. 2016;Abdulaziz et al. 2019b), where y 𝑙,𝑏 ∈ C 𝑀 𝑙,𝑏 is composed of 𝑀 𝑙,𝑏 visibilities. The corresponding decomposition of the measurement operator is hereafter denoted by 𝚽 𝑙 = (𝚽 𝑙,𝑏 ) 1≤𝑏≤ 𝐵 .

Faceted HyperSARA

Optimization problem

From the perspective of optimization theory, reconstructing the original wideband image X = (x 𝑙 ) 1≤𝑙 ≤ 𝐿 from incomplete noisy Fourier measurements can be formulated as the resolution of the following constrained problem minimize

X=(x 𝑙 ) 1≤𝑙 ≤ 𝐿 ∈R 𝑁 ×𝐿 + 𝐿 ∑ 𝑙=1 𝐵 ∑ 𝑏=1 𝜄 B (y 𝑙,𝑏 , 𝜀 𝑙,𝑏 ) 𝚽 𝑙,𝑏 x 𝑙 + 𝑟 (X), (2)
where the indices (𝑙, 𝑏) ∈ {1, . . √ 𝑀 𝑙 , reflecting the statistics of the noise in the corresponding data block, assuming a Chi-squared distribution (see Carrillo et al. 2012;Onose et al. 2016, for further details). The function 𝑟 denotes a regularization term. Finally, a non-negativity constraint is imposed on the wideband image, as it corresponds to an intensity image.

In the minimized objective function, the indicator functions 𝜄 B (y 𝑙,𝑏 , 𝜀 𝑙,𝑏 ) act as data-fidelity terms, ensuring the consistency of the modeled visibilities with the measurements (Carrillo et al. 2012), whereas the function 𝑟 encodes a prior model for the unknown image cube to address the ill-posedness of the inverse problem described in (1). Note that the prior function involves regularization parameters, which are set from the noise level approximated through the heuristics proposed in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF]. Both the HyperSARA and Faceted HyperSARA approaches rely on a low-rankness and average joint-sparsity prior (see [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], Section 2)), 1 Let C be a non-empty, closed, convex subset of C 𝑁 , then 𝜄 C denotes the indicator function of C, defined by 𝜄 C (z) = +∞ if z ∈ C and 0 otherwise.

MNRAS 000, 1-15 (2022) whereas the SARA approach relies on a sparsity prior. The Faceted HyperSARA approach is summarized in the next paragraph.

Faceted HyperSARA approach

Building on the original HyperSARA approach, Faceted Hyper-SARA relies on the two distinct levels of faceting described below.

Spectral faceting. The wideband image cube can be reconstructed as a collection of 𝐶 image sub-cubes, each recovered from a subset of the frequency channels and leveraging an independent model prior. In [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], channel-interleaved spectral sub-cubes have been considered to ensure that each imaging problem can exploit the spectral information over the extent of the observed frequency band. In this setting, a collection of 𝐶 problems are solved, where the datafidelity terms are frequency channel-based and the regularization terms acts on the image sub-cubes separately. Solving the overall objective function involved in (2) reduces to solving the sum of independent objective functions associated with the 𝐶 sub-cubes as follows minimize

X 𝑐 ∈R 𝑁 ×𝐿𝑐 + 𝐿 𝑐 ∑ 𝑙=1 𝐵 ∑ 𝑏=1 𝜄 B (y 𝑐,𝑙,𝑏 , 𝜀 𝑐,𝑙,𝑏 ) 𝚽 𝑐,𝑙,𝑏 x 𝑐,𝑙 + 𝑟 𝑐 (X 𝑐 ), (3)
where 𝑐 ∈ {1, . . . , 𝐶} and 𝐿 𝑐 denotes the number of channels of 𝑐-th sub-cube. For each channel 𝑙 ∈ {1, . . . , 𝐿 𝑐 } and data-block 𝑏, y 𝑐,𝑙,𝑏 ∈ C 𝑀 𝑐,𝑙,𝑏 denotes the vector of 𝑀 𝑐,𝑙,𝑏 visibilities, 𝚽 𝑐,𝑙,𝑏 the associated measurement operator, and 𝜀 𝑐,𝑙,𝑏 the ℓ 2 ball radius. Finally, the matrix

X 𝑐 = (x 𝑐,𝑙 ) 1≤𝑙 ≤ 𝐿 𝑐 ∈ R 𝑁 ×𝐿 𝑐 is the 𝑐-th sub-cube of the full image cube X = (X 𝑐 ) 1≤𝑐≤𝐶 , with x 𝑐,𝑙 ∈ R 𝑁 the 𝑙-th image of the sub-cube X 𝑐 .
The regularization term 𝑟 𝑐 : R 𝑁 ×𝐿 𝑐 → ] -∞, +∞] acts on the 𝑐-th image sub-cube. In addition to nonnegativity, it enforces low-rankness and average joint-sparsity of the image sub-cube, encapsulated in non-convex log-sum functions (see Thouvenin et al. 2021, Section 2.3, for full details).

Spatial faceting.

The regularization term 𝑟 𝑐 from (3) can be further decomposed into a sum of 𝑄 spatial facet-specific terms, i.e., exclusively affecting contiguous pixels. On the one hand, the wavelet dictionaries involved in the average joint-sparsity prior from Hyper-SARA can be exactly decomposed over overlapping facets (Pruša 2012, Chapter 4). The associated regularization term can thus be reformulated as a sum of several per-facet terms. On the other hand, a faceted low-rankness prior is considered to provide a more tractable alternative to the global low-rankness prior from HyperSARA involving the singular value decomposition (SVD) of the full image cube.

On a final note, the HyperSARA approach can be seen as a particular case of the Faceted HyperSARA approach with (𝐶, 𝑄) = (1, 1).

Algorithm

To efficiently solve (3) and address its non-convex regularization terms, a reweighting iterative scheme (Candès et al. 2009) relying on the primal-dual forward-backward (PDFB) algorithmic structure (Condat 2013;[START_REF] Vũ | [END_REF]Pesquet & Repetti 2015) is considered (see Thouvenin et al. 2021, Section 3). More precisely, the outer reweighting algorithm allows the resolution of the non-convex problem via the resolution of a sequence of convex sub-problems (see [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], Section 3.1, Algorithm 1). These subproblems are solved efficiently using PDFB as inner algorithmic structure, where full parallelization of all block-specific data-fidelity terms and facet-specific regularization terms is achieved (see Thouvenin et al. 2021, Section 3.2, Algorithm 2).

Compared to HyperSARA, Faceted HyperSARA thus offers more parallelization flexibility in that:

(i) spectral faceting can be used to independently reconstruct the wideband image in parallel over non-overlapping groups of channels, which lead to 𝐶 independent optimization problems of the form (3), formulated for a single wideband sub-cube X 𝑐 . This feature can be used to exploit more computational resources for the outer reweighting algorithm;

(ii) spatial faceting allows a larger number of workers to be involved in the inner PDFB algorithm underlying the approach via splitting. Furthermore, per-facet regularization terms act on a smaller part of the wideband image, and can be assigned to workers on a perfacet basis. Further details about the communications between the workers involved in the PDFB algorithm are provided in Thouvenin et al. (2021, Section 3.3).

TOWARDS AN ACCURATE AND SCALABLE MEASUREMENT OPERATOR

In this section, we describe new functionalities compared to [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF]. These are aimed at (i) enhancing the accuracy of the RI measurement operator by incorporating available models of the DDEs and correcting for the 𝑤-term in the context of wide-field imaging, (ii) addressing the limitations of the imperfect calibration and its impact on the noise level estimate and (iii) improving the scalability potential of the measurement operator via data dimensionality reduction. These features are crucial to achieve high precision imaging from large data volumes.

Accounting for estimates of the direction-dependent effects

Based on a compact-support model in the spatial Fourier domain, DDE estimates can be accounted for in the measurement operator when building the de-gridding matrix G (Repetti et al. 2017). More precisely, each row of G, associated with a data point, results from a linear convolution between the compact Fourier representation of the DDEs of the associated antenna pair and a non-uniform Fourier interpolation kernel [START_REF] Dabbech | [END_REF]. The particular case of the direction independent effects (DIE), where the antenna gains are modeled via complex values, is addressed in the imaging by applying the available DIE solutions directly to the data.

Correction of the 𝑤-term via 𝑤-projection

The non-coplanarity of the radio array results in variations in the projections of the baselines onto the direction of observation, defining their 𝑤-coordinates. The so called 𝑤-term -a phase modulation in the image domain -becomes non-negligible, especially for wide-field imaging (Cornwell et al. 2005). To take this geometric DDE into account in the measurement operator, we consider the 𝑤-projection approach (Cornwell 2008). Each visibility is modeled from the Fourier transform of the radio map using a linear convolution with the spatial Fourier transform of the 𝑤-term. Building on [START_REF] Wolz | [END_REF]; Dabbech et al. (2017), we consider compact 𝑤-terms with adaptive resolution based on their estimated spatial bandwidth. The resulting 𝑤-kernels are further sparsified via adaptive sparsification procedures. The kernels are incorporated into the de-gridding matrix G via a linear convolution with the non-uniform Fourier interpolation kernel, resulting in accurate and fast measurement operators.

Effective noise level estimation

Given the implication of the noise level in the definition of the minimization problem through (i) the ℓ 2 bounds (𝜀 𝑙,𝑏 ) (1≤𝑙 ≤ 𝐿,1≤𝑏≤ 𝐵) involved in the data-fidelity constraints in (2), and (ii) the regularization parameters resulting from the heuristics proposed in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], it is of paramount importance to provide reliable noise estimates to the intensity imaging algorithm, to ensure its effectiveness in terms of reconstruction quality. In the case of perfect calibration, the standard deviation of the naturally-weighted data is assumed to be equal to one. Yet, in practice, calibration errors tend to compromise the accuracy the measurement operator, resulting into an effective noise level higher than the theoretical instrumental noise, especially in the context of the highly sensitive modern telescopes. In this case, estimates of the effective noise can be obtained by (i) relying on model visibilities obtained from a calibration preprocessing step, where the noise level per each channel is derived as the standard deviation of the associated residual visibilities or (ii) activating a new functionality for adaptive estimation of the effective noise level within a few PDFB iterations inspired by Dabbech et al. (2018); Abdulaziz et al. (2019b). Starting from inactive image priors, the strategy consists in progressively updating the noise level estimates and associated parameters (i.e. ℓ 2 bounds and regularization parameters), throughout the iterative process. Once the ℓ 2 bounds are satisfied, the noise-related parameters are fixed, thus defining the final problem to be solved, and the intensity imaging algorithm is run to convergence.

Data dimensionality reduction

To handle large data volumes, the RI data and, consequently, the associated measurement operator, can be embedded into a lowerdimensional sub-space, especially when the number of measurements in each channel exceeds the spatial dimension of the target model cube. Kartik et al. (2017) proposed various data dimensionality reduction schemes. In this paper, we consider a visibility griddingbased scheme due to its simplicity. An overview of the technique is given in the following lines in the context of Faceted HyperSARA.

Consider 𝐶 independent channel-interleaved sub-problems, and the linear embedding operators 𝚪 𝑐,𝑙,𝑏 ∈ C 𝑀 ′ 𝑐,𝑙,𝑏 × 𝑀 𝑐,𝑙,𝑏 . The forward model (1) after data gridding reads

y ′ 𝑐,𝑙,𝑏 = 𝚽 ′ 𝑐,𝑙,𝑏 x 𝑙 + n ′ 𝑐,𝑙,𝑏 , (4)
where 𝚽 ′ 𝑐,𝑙,𝑏 = 𝚪 𝑐,𝑙,𝑏 𝚽 𝑐,𝑙,𝑏 , and

y ′ 𝑐,𝑙,𝑏 = 𝚪 𝑐,𝑙,𝑏 y 𝑐,𝑙,𝑏 ∈ C 𝑀 ′ 𝑐,𝑙,𝑏
and n ′ 𝑐,𝑙,𝑏 = 𝚪 𝑐,𝑙,𝑏 n 𝑐,𝑙,𝑏 ∈ C 𝑀 ′ 𝑐,𝑙,𝑏 are the gridded data and noise vectors, respectively. In this formulation, the resulting measurement operator 𝚽 ′ 𝑐,𝑙,𝑏 ∈ C 𝑀 ′ 𝑐,𝑙,𝑏 × 𝑁 becomes blind to the original data dimension 𝑀 𝑐,𝑙,𝑏 , and the reduced data dimension 𝑀 ′ 𝑐,𝑙,𝑏 can be significantly smaller than both the original data and image dimensions (𝑀 𝑐,𝑙,𝑏 , 𝑁). The gridding operator reads 𝚪 𝑐,𝑙,𝑏 = 𝚲 𝑙,𝑏,𝑐 G † 𝑐,𝑙,𝑏 𝚯 † 𝑐,𝑙,𝑏 , with the diagonal matrix 𝚲 𝑐,𝑙,𝑏 ∈ R 𝑀 ′ 𝑐,𝑙,𝑏 × 𝑀 ′ 𝑐,𝑙,𝑏 acting as a noise-whitening operator. The resulting measurement operator 𝚽 ′ 𝑐,𝑙,𝑏 explicitly reads

𝚽 ′ 𝑐,𝑙,𝑏 = 𝚲 𝑐,𝑙,𝑏 H 𝑐,𝑙,𝑏 FZ, (5)
that is encapsulating the sparse holographic matrix given by

H 𝑐,𝑙,𝑏 = G † 𝑐,𝑙,𝑏 𝚯 † 𝑐,𝑙,𝑏 𝚯 𝑐,𝑙,𝑏 G 𝑐,𝑙,𝑏 ∈ C 𝑀 ′ 𝑐,𝑙,𝑏 × 𝑀 ′
𝑐,𝑙,𝑏 . The diagonal matrix 𝚲 𝑐,𝑙,𝑏 acts effectively on the nonzero rows of the matrix H 𝑐,𝑙,𝑏 , with its entries corresponding to the inverse of the square root of the nonzero diagonal elements of the matrix H 𝑐,𝑙,𝑏 and zero otherwise. On the one hand, this form of weighting ensures that the embedded noise is approximately an independent and identically distributed Gaussian noise. On the other hand, it operates as a uniform weighting, as the diagonal elements of the holographic matrix encapsulate the density of the sampling of the gridded visibilities (Kartik et al. 2017).

MATLAB LIBRARY

In this section, we describe our proposed MATLAB library, that is focused on RI wideband intensity imaging. The library offers a collection of utility functions and scripts from data extraction from an RI measurement set (MS Table ) to the reconstruction of a wideband intensity image over the field of view and frequency range of interest. The library is composed of three modules associated with (i) the formation of the measurement operator, including the new functionalities presented in Section 3, (ii) a faceted implementation of the sparsity dictionary involved in the average sparsity prior of the SARA family, and (iii) the imaging pipeline which enables interfacing with an RI measurement set and wideband intensity imaging via any of the three following approaches: SARA, HyperSARA, and Faceted HyperSARA. The Faceted HyperSARA library (Thouvenin et al. 2022a) is accessible from the Puri-Psi webpage [START_REF] Wiaux | Puri-Psi library: Parallel Proximal Scalable Imaging[END_REF].

Library structure

The proposed library is composed of a core module, dubbed Faceted-HyperSARA, and two standalone sub-modules which are RI-measurement-operator and SARA-dictionary. Fig. 1(a) provides an illustration of the structure of the library. Details of the three modules are provided in what follows.

RI-measurement-operator. The module covers all the functionalities associated with the formation of the measurement operator. These include the non-uniform Fourier transform, 𝑤-correction via 𝑤-projection, incorporation of available DDE solutions, and data dimensionality reduction via a visibility gridding-based scheme. The module is available at https://basp-group.github.io/RI-m easurement-operator. SARA-dictionary. The library provides faceted and non-faceted implementations of the SARA dictionary and comes with a detailed documentation, selected tests, and code examples. The module is available at https://basp-group.github.io/SARA-diction ary.

Faceted-HyperSARA. The core module, illustrated in Fig. 1(b), enables the following procedures. (ii) Imaging. The library provides the implementation of the solvers (i.e. algorithmic structures) underpinning the SARA, Hyper-SARA and Faceted HyperSARA approaches. Compared to [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], the solvers can seamlessly be applied to the inverse problem resulting from data dimensionality reduction via visibility gridding. The user interacts with the solvers solely through the dedicated imaging script main_input_imaging.m. The script requires a few user-defined parameters, described in Section 4.2 and calls the imaging function imaging.m which uses reliable default algorithmic parameters whenever left unspecified by the user. Finally, the library systematically uses the heuristic regularization parameters described in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], Section 2) for the SARA, HyperSARA, and Faceted HyperSARA approaches, respectively.

Associated documentation is available online (Thouvenin et al. 2022a).

Input parameters

The imaging pipeline depends on a number of parameters, which can be classified into the three following groups, depending on the level of interaction expected from the user.

User-defined parameters. In order to run an image reconstruction task, a few parameters need to be specified in the main_input_imaging.m from the Faceted-HyperSARA module. These are specifically related to (i) the input data (e.g. the name of target source, as used in the data extraction step, and the indices of the frequency channels to image), (ii) the output image cube (e.g. the image size and the pixel size), (iii) the definition of the optimization problem (2) (e.g., the selected approach among SARA, HyperSARA, and Faceted HyperSARA, and the scale factors to be applied to the heuristic values of the associated regularization parameters if different from 1), and (iv) algorithmic parameters (e.g. the number of facets with the associated overlap for Faceted HyperSARA). A full list of these parameters is provided in the online documentation (Thouvenin et al. 2022a).

Default internal parameters.

Several parameters, internal to the proposed imaging approach, are set to reliable default values in the auxiliary file default_parameters.json. These include the kernel type and size taken for the NUFFT, and all the algorithmic parameters necessary to configure the reweighting procedure and inner PDFB algorithm (see [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], Algorithms 1 and 3)). A complete list of these parameters with the associated description is provided in the associated documentation (Thouvenin et al. 2022a).

Users are invited to rely on the default parameter configuration, unless specifically aware of their impact on the image reconstruction task. The default values can however be overridden by creating a new .json file, passed to the imaging script. In this case, users are strongly encouraged to start from the default configuration file and edit the values to be modified.

Optional parameters. The library can take advantage of the results of a preliminary calibration step, where the pre-estimated DDE solutions are incorporated into the measurement operator or, alternatively, the DIE solutions are applied directly to the data (see Section 3.1).

Residual visibilities can be leveraged to obtain estimates of the effective noise level as described in Section 3.3. Finally, an estimate of the image can be used to initialize the intensity imaging algorithm, for acceleration purposes.

DATA DESCRIPTION AND EVALUATION METRICS

To validate the proposed approach and MATLAB library at scale, a large RI wideband dataset is considered. Prior to the wideband intensity imaging, a pre-processing step, consisting in the monochromatic joint calibration and imaging is performed to enhance the accuracy of the measurement operator model, thus enabling higher precision imaging. This section provides the details of the utilized dataset and the pre-processing step. Evaluation criteria considered to assess the compared approaches are presented.

Data description

The measurements for configurations A and C, respectively), processed as double precision complex numbers. The data were flagged in CASA (McMullin et al. 2007) in two steps. Firstly, known spectral regions affected by satellite transmission were removed. Secondly, sporadic data affected by radio frequency interference were removed using a simple thresholding procedure. The data have then been selfcalibrated on a spectral window-basis, by alternating between an imaging step, which consists in the formation of a single image from the full spectral window, and a calibration step.

Data pre-processing: monochromatic joint calibration and imaging

To improve the accuracy of the modeled measurement operator, a joint calibration and imaging pre-processing step was conducted for each frequency channel, separately [START_REF] Dabbech | [END_REF]). The approach, originally proposed by Repetti et al. (2017), consists in the alternate estimation of the unknown DDEs and the image of interest, where a spatio-temporal smoothness prior is adopted for the DDEs, and an average sparsity prior in addition to non-negativity is considered for the image (Repetti & Wiaux 2017;[START_REF] Thouvenin | [END_REF].

Ideally, a wideband joint calibration and imaging approach would have been the most appropriate. Although preliminary work in this direction has been recently proposed [START_REF] Dabbech | International BASP Frontiers workshop[END_REF]), a fully parallelizable algorithm remains to be developed in this context. Following the discussion in Section 2.1, the calibration solutions can be easily integrated into the forward model (1). Furthermore, the resulting estimated model visibilities can be leveraged to estimate the noise statistics, and consequently, the heuristic regularization parameters introduced in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], Section 2), as well as estimates of the ℓ 2 bounds defining the data-fidelity terms in (2), taken as the ℓ 2 norms of the residual visibilities.

Evaluation metrics

Imaging quality is evaluated by visually inspecting the reconstructed images2 , in comparison with SARA (Onose et al. 2017) and the state-of-the-art wideband CLEAN-based method implemented in WSClean -a C++ imager for RI (Offringa & Smirnov 2017). Note that HyperSARA (Abdulaziz et al. 2019b) is not considered as a benchmark approach due to its prohibitive runtime.

For both Faceted HyperSARA and SARA, the analysis is conducted using the estimated model cube X and the naturally-weighted residual image cube R whose columns, indexed by 𝑙 ∈ {1, . . . , 𝐿}, are given by r 𝑙 = 𝜂 𝑙 𝚽 † 𝑙 (y 𝑙 -𝚽 𝑙 x 𝑙 ), with y 𝑙 the naturally-weighted RI measurements. Let 𝜹 ∈ R 𝑁 denote an image with value 1 at the phase center and zero elsewhere. The normalization factor 𝜂 𝑙 is such that the per-channel point spread function (PSF), given by 𝜂 𝑙 𝚽 † 𝑙 𝚽 𝑙 𝜹, has a peak value equal to 1. For WSClean, optimal results are obtained with Briggs weighting (Briggs 1995). In addition to the Briggs-weighted residual cube R = ( r 𝑙 ) 1≤𝑙 ≤ 𝐿 , the restored image cube T = (t 𝑙 ) 1≤𝑙 ≤ 𝐿 is considered. Its columns are defined as t 𝑙 = x 𝑙 * c 𝑙 + r 𝑙 , for each channel index 𝑙 with x 𝑙 being the estimated model image (consisting of the CLEAN components) and c 𝑙 is the CLEAN restoring beam (typically a Gaussian fitted to the primary lobe of the associated PSF). We also consider in our analysis the image cube defined as the difference between the restored cube and the residual cube, whose columns are the model images convolved with the associated restoring beams. As a quantitative measure of data fidelity, the standard deviation of the recovered residual images is reported.

All the methods compared in this paper were run on compute nodes of Cirrus, one of the UK's Tier2 HPC services3 . Cirrus is an SGI ICE XA system composed of 280 compute nodes, each with two 2.1 GHz, 18-core, Intel Xeon E5-2695 (Broadwell) series processors. The compute nodes have 256 GB of memory shared between the two processors. The system has a single Infiniband FDR network connecting nodes with a bandwidth of 54.5 GB/s.

Computing cost is assessed in terms of the following metrics (i) runtime per iteration (run pips ). For SARA and Faceted Hy-perSARA, this metric is reported per PDFB iteration per sub-cube;

(ii) (active) CPU time (cpu pips ), omitting the communication time between workers. For SARA and Faceted HyperSARA, the metric is reported per PDFB iteration per sub-cube;

(iii) reconstruction runtime (maximum over all subcubes), assuming all the subcubes are reconstructed in parallel (run). This metric includes the communication time required between the facet and data workers. Note that the subcubes are reconstructed independently in parallel, as no communication is required between these imaging tasks;

(iv) total (active) CPU time used to reconstruct the full image cube (cpu).

Per iteration metrics are reported in terms of their average and associated standard deviation over all iterations of the algorithm. On a final note, the pre-processing step is not included in the computing time reported for the different approaches.

EXPERIMENT 1: WIDEBAND IMAGING AT FULL SPECTRAL RESOLUTION

In this section, we study the imaging performance and the computational cost of Faceted HyperSARA via the reconstruction of a 15 GB image cube of Cyg A from 7.4 GB of VLA data. Leveraging both the spectral and spatial faceting functionalities, the algorithm is mapped on 1440 CPU cores on the high performance computing system Cirrus. The performance of the proposed approach is evaluated in comparison with the monochromatic imaging approach SARA (Onose et al. 2017). In this experiment, DDE estimates obtained from the pre-processing step described in Section 5.2 are considered. Since these estimates cannot be directly transferred to the WSClean software, the wideband CLEAN-based method is not considered for this experiment.

Experimental setting

Imaging is conducted over the field-of-view (FoV) Ω 1 = 2.56 ′ × 1.536 ′ at the spatial resolution given by the pixel size 𝛿𝑥 = 0.06 ′′ (in both directions), leading to 𝑁 = 1536×2560. Full spectral resolution is considered, i.e., 𝐿 = 480 channels. Hence, the wideband imaged cube is of size 1536 × 2560 × 480. The selected pixel size is such that the spatial bandwidth of the recovered signal is up to 1.75 times the instrumental resolution at the highest channel 𝜈 𝐿 = 8.019 GHz, and 3.53 times the instrumental resolution at the lowest channel 𝜈 1 = 3.979 GHz. Considering the same FoV and spatial resolution, the monochromatic joint calibration and imaging pre-processing step was performed separately on each frequency channel, providing antenna based and time dependent DDE estimates of spatial dimension 5 × 5 in the spatial Fourier domain [START_REF] Dabbech | [END_REF]. For all (𝑐, 𝑙) ∈ {1, . . . , 𝐶} × {1, . . . , 𝐿 𝑐 }, these estimates were incorporated in the de-gridding matrices G 𝑐,𝑙 (see related discussion in Section 2.1). From now on, 𝚽 𝑐,𝑙 refer to the resulting measurement operators.

In both imaging algorithms, 𝐵 = 2 data blocks per channel are considered. The blocks are associated with the VLA configurations A and C, each having different noise statistics. For each data block indexed by (𝑐, 𝑙, 𝑏) ∈ {1, . . . , 𝐶} × {1, . . . , 𝐿 𝑐 } × {1, . . . , 𝐵}, the noise level 𝜎 𝑐,𝑙,𝑏 is estimated as the standard deviation of the residual visibilities y 𝑐,𝑙,𝑏 -𝚽 𝑐,𝑙,𝑏 x (0) 𝑐,𝑙 , where x (0) 𝑐,𝑙 is the model image obtained from the pre-processing step. The ℓ 2 bounds defining the data-fidelity constraints in (3) are set to 𝜀 𝑐,𝑙,𝑏 = ∥y 𝑐,𝑙,𝑏 -𝚽 𝑐,𝑙,𝑏 x (0) 𝑐,𝑙 ∥ 2 . More specifically to Faceted HyperSARA, 𝐶 = 16 channel-interleaved sub-problems are considered with 𝐿 𝑐 = 30 for any 𝑐 ∈ {1, . . . , 𝐶}, using 𝑄 = 10 × 6 facets along the spatial dimension, with an overlap of about 50%. The resulting number of spatio-spectral facets is 𝑄 ×𝐶 = 960. Given the estimates of the noise level, the regularization parameters of both algorithms are set from the heuristics proposed in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF] as follows. In Faceted HyperSARA, the lowrankness regularization parameter is fixed to the associated proposed heuristic, and the joint sparsity regularization parameter is set to 3 times the associated heuristic for optimal dynamic range. The regularization parameter of SARA is set to 3 times the proposed heuristic for the same reason. This deviation from the heuristics proposed in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF] can be explained by the need to compensate for residual calibration errors affecting the accuracy of the measurement operator model. For both imaging algorithms, a maximum of 5 reweighting steps is considered. To reduce the runtime to convergence, one could initialize the algorithms with the wideband model cube X (0) = (X (0) 𝑐 ) 1≤𝑐≤𝐶 obtained from the joint DDE calibration and imaging step. However, in this experiment, to avoid any bias in the reported computational cost, both imaging algorithms have been initialized with a zero wideband cube.

Results and discussion

Imaging quality

To assess the imaging quality of Faceted HyperSARA in comparison with SARA, we examine the estimated images of the selected frequency channels 𝜈 102 = 4.787 GHz and 𝜈 459 = 7.851 GHz, displayed in Figures 2 and3 respectively. The channels correspond to the respective channel indices 7 and 29 of sub-problems 6 and 11. The model images, focusing on Cyg A, are displayed in log 10 scale, and overlaid with the corresponding residual images. Zooms on selected regions of the radio galaxy have been included to enhance the visual comparison.

Although both algorithms yield very comparable reconstructions for high-intensity pixels (≥ 1 mJy, see the hotspot regions), Faceted HyperSARA (top rows of Figures 23) supersedes SARA (bottom rows of Figures 23) at low intensity by capturing a more physicallyconsistent structure. For instance, the annulus at the high end of the western lobe is consistently recovered by Faceted HyperSARA at both frequency channels, while it is fully dominated by artefacts in the form of stripes in the SARA model image at the low frequency channel. This behaviour is also observed in other regions of the radio galaxy (e.g. see both zooms of Figure 3). Moreover, the improved dynamic range of Faceted HyperSARA can be observed in the recovery of faint emission (≤ 0.01 mJy), such as additional segments of the jets and the faint extended emission at the ends of both lobes (see both zooms of Figure 2). These emissions were previously detected in Cyg A image at S band formed from highly sensitive monochromatic VLA observations [START_REF] Dabbech | [END_REF]. The superior imaging quality of Faceted HyperSARA is the result of the correlation imposed across channels, making it more robust to calibration errors than its monochromatic counterpart. Note that calibration artefacts observed in SARA images can be reduced by enforcing more sparsity through a higher regularization parameter, yet at the expense of a much smoother reconstruction. Importantly, the spatial tessellation underlying Faceted HyperSARA image model does not introduce any discernible artefact in the estimated images.

Full FoV residual images obtained with the two algorithms are Full image cubes are made available online in FITS and GIF formats (Thouvenin et al. 2022b). Their visual examination indicates a generally smooth spectral behaviour of Faceted HyperSARA wideband image, in line with expected power law decays associated with synchrotron emission. In contrast, SARA wideband image exhibits prominent channel-varying fluctuations mostly in the form of stripe-like artefacts, confirming the above analysis. Although a spectral index map can be directly inferred from the reconstructed images, such analysis is beyond the scope of this work. Nonetheless, we provide a spectral analysis based on selected emissions from different regions in the radio galaxy. For Faceted HyperSARA, the spectra at different pixel positions are generally consistent with a continuum emission following a typical power-law decay, with the exception of a noticeable spectral discontinuity observed at channel 𝜈 257 = 5.979 GHz. Since the data were acquired separately in the frequency ranges [𝜈 1 , 𝜈 256 ] = [3.979, 6.019] GHz and [𝜈 257 , 𝜈 𝐿 ] = [5.979, 8.019] GHz, this discontinuity is most likely due to differences in the calibration errors and noise statistics between the two frequency ranges. A small-scale spectral discontinuity demarcating the spectral windows can be further noticed, which most likely results from calibration errors induced by the initial calibration procedure conducted on a spectral window basis (see Section 5.1). As for SARA, additional irregularities of the spectra at different pixel positions are observed, particularly important at pixel positions with intensity values below 1 mJy, reflecting the channel-varying fluctuations due to the inherent calibration errors and the lack of correlation across the channels.

Computing cost

Faceted HyperSARA and SARA were run on multiple nodes of Cirrus. For Faceted HyperSARA, each sub-problem (3) (composed of 30 channels) was solved using 91 CPU cores, distributed as follows; 1 master CPU core, 30 CPU cores to process the 2 × 30 data-fidelity terms (four data-fidelity blocks per core), and 60 CPU cores to handle the 60 spatio-spectral facets. As for SARA, the method was applied in the parallel configuration proposed by Onose et al. (2016) with 𝐶 = 𝐿 and 𝑄 = 1 ) in terms of: number of CPU cores, average number of PDFB iterations (across the 𝐶 independent problems), runtime per iteration per sub-cube and total runtime (run pips in seconds with the associated standard deviation, run in hours), CPU time per iteration per subcube and total CPU time (cpu pips in seconds with the associated standard deviation, cpu in hours). The total memory allocated for the computation of G 𝑐,𝑙,𝑏 and the memory needed to store them are reported.

precisely, each channel was imaged using 12 CPU cores, nine of which were exploited by the average sparsity terms (associated with the nine bases of the SARA dictionary). The computing time and resources required by Faceted HyperSARA and SARA are summarized in Table 2. Faceted HyperSARA required nearly 50% more resources than SARA in terms of cpu hours in this case. This could be explained by the larger number of workers used for the Faceted HyperSARA regularization and the higher number of iterations required by the algorithm to reach convergence, as already observed in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF].

EXPERIMENT 2: WIDEBAND IMAGING WITH DIMENSIONALITY REDUCTION

In the previous section, spectral faceting has been leveraged to divide the imaging problem into smaller, independent sub-problems.

Faceting allowed a wideband model cube of Cyg A to be formed at the full spectral resolution of the data. The image cube estimate exhibits a relatively flat spectrum across the frequency range of interest. This observation suggests that a much coarser spectral resolution could have been adopted a priori without violating the Nyquist limit. Data points from consecutive channels can therefore be combined to form an image cube with a lower spectral dimension. In this setting, the amount of data per effective channel increases significantly, leading to large memory requirements to reconstruct the image cube. Memory costs required to store the de-gridding matrix G involved in the measurement operator (up to 2.5 TB for the 7.4 GB data utilized in Section 6) are particularly large. To tackle this challenge, we leverage the data dimensionality reduction procedure relying on a visibility gridding-based scheme described in Section 3.4. The resulting joint image and data dimensionality reduction can significantly reduce both memory requirements and computing cost to form RI image cubes. In this section, we study the performance of Faceted Hyper-SARA in comparison with the CLEAN-based wideband imaging algorithm in WSClean.

Experimental setting

Given the reconstruction results reported in Section 6.2 at the spectral resolution 𝛿𝜈 = 8 MHz, we consider herein a coarser spectral sampling at the resolution 𝛿𝜈 = 128 MHz. In other words, the 16 consecutive channels composing each spectral window are combined to form a single effective channel image. In this setting, a joint DIE calibration and imaging pre-processing step was applied to each of the effective channels separately, with the aim to improve the accuracy of the measurement operator while ensuring a fair comparison with the state-of-the-art CLEAN-based wideband approach in WSClean. As in the previous experiment, the FoV adopted in the pre-processing step is Ω 1 = 2.56 ′ × 1.536 ′ , with the spatial resolution given by the pixel size 𝛿𝑥 = 0.06 ′′ (in both directions). Antenna gain solutions, which are now reduced to complex values, were applied directly to the data. The corrected data were then easily stored in the measurement sets and thus leveraged by both Faceted HyperSARA and the wideband imaging approach in WSClean. Note that the 29 th spectral window (spanning the frequency range [7.767, 7.895] GHz) was removed due to large calibration errors. The target wideband cube is therefore composed of 𝐿 = 29 effective channels. Despite the absence of DDE correction, this pre-processing step results in a good image reconstruction quality as will be shown in the following paragraphs.

For wideband imaging, we consider the FoV Ω 2 = 2.1709 ′ × 3.6181 ′ , with a spatial resolution given by the pixel size 𝛿𝑥 = 0.0424 ′′ (in both directions). The spatial bandwidth of the recovered signal is up to 2.4906 times the instrumental resolution at the highest effective channel 𝜈 29 = 7.959 GHz, and 4.9245 times the instrumental resolution at the lowest effective channel 𝜈 1 = 4.039 GHz. This wideband imaging setting is different from the pre-processing step where a coarser resolution and smaller FoV are considered to reduce its computational cost. The target wideband image cube is of size 𝑁 × 𝐿 = 3072 × 5120 × 29. Given the significant reduction in the spectral dimension, only one spectral facet is considered in Faceted HyperSARA (i.e. 𝐶 = 1, and 𝐿 𝐶 = 𝐿 = 29). The index 𝑐 will thus be omitted in the following paragraphs for simplicity.

Initially composed of about 4.5 × 10 8 points, the data are reduced down to 6.5 × 10 7 points after gridding, that is about 7-fold reduction in data volume, from 7.2 GB down to 1.04 GB. Since DIE solutions are adopted in this experiment, the memory required to store the de-gridding matrices (G 𝑙,𝑏 ) (1,1) ≤ (𝑙,𝑏) ≤ ( 𝐿,𝐵) is about 530 GB, a decrease of a factor ≈ 4.55 in comparison with the experiment conducted in Section 6. The data dimensionality reduction yields a decrease in the memory requirements of the resulting measurement operator, with the underlying holographic matrices (𝚲 𝑙,𝑏 H 𝑙,𝑏 ) (1,1) ≤ (𝑙,𝑏) ≤ ( 𝐿,𝐵) occupying about 108 GB, a 4.8-fold reduction of the memory needed to store the de-gridding matrices.

For each effective channel, the data associated with the two VLA configurations and the different channels forming a spectral window are aggregated into 𝐵 = 1 block after gridding. The noise level after gridding and the associated data constraints are set from the model visibilities ȳ𝑙,𝑏 , obtained from the preprocessing step. More precisely, for each data block indexed by (𝑙, 𝑏) ∈ {1, . . . , 𝐿} × {1, . . . , 𝐵}, the associated ℓ 2 bound is estimated as 𝜀 ′ 𝑙,𝑏 = ∥y ′ 𝑙,𝑏 -𝚲 𝑙,𝑏 G † 𝑙,𝑏 𝚯 † 𝑙,𝑏 ȳ𝑙,𝑏 ∥ 2 . Faceted HyperSARA was applied with 𝑄 = 10 × 6, i.e. 𝑄 × 𝐶 = 60 spatio-spectral facets with an overlap of about 50%. Using the estimates of the noise level, the regularization parameters are fixed using the heuristics proposed in [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF]. The lowrankness regularization parameter is set to the proposed heuristic, whereas the joint sparsity regularization parameter is set to 0.1 times its associated heuristic for optimal resolution. Once again, the choice of the regularization parameter associated with the average joint sparsity prior deviates from the proposed heuristic, which can be explained by the inaccuracy of the measurement operator in the presence of residual calibration errors. In the re-weighting procedure, a maximum of 5 reweighting steps was considered. Finally, Faceted HyperSARA was initialized with a zero wideband cube, since the monochromatic joint DIE calibration and imaging pre-processing step was conducted at a coarser resolution and over a smaller FoV compared to the wideband imaging step.

For WSClean, Briggs weighting with the robust parameter set to -0.25 was considered to enhance the overall resolution. The wideband deconvolution was performed using a spectral polynomial fit of order 3 and 10 6 maximum number of iterations. The command used in WSClean is reported in Appendix A.

Results and discussion

Imaging quality

Reconstructed images at two selected channels 𝜈 10 = 5.191 GHz and 𝜈 29 = 7.959 GHz are reported in Figures 4 and5 respectively. From top to bottom, are shown the estimated model image of Faceted HyperSARA, the restored image of the CLEAN-based approach in WSClean and the associated model image convolved with a restoring beam and obtained by subtracting the residual image from the restored image. All the images are displayed in log 10 scale. WS-Clean images are normalized by the flux of their associated restoring beam, so that all images are in the same unit Jy/pixel. In addition, the negative pixels in WSClean images are set to zero for visualization purposes. All images are displayed with zooms on selected regions of the radio galaxy to better assess their dynamic range and resolution. Finally, full FoV residual images, displayed in linear scale, are overlaying the respective estimated model and restored images of Faceted HyperSARA and WSClean.

Looking at the estimated model images of Faceted HyperSARA (top rows of Figures 45) in comparison with the restored images of WSClean (middle rows of Figures 45), one can clearly see the superior performance of the proposed imaging algorithm, in terms of both resolution (see the zoomed region from the eastern hotspot, displayed on the bottom right of each panel of Figure 4) and dynamic range, particularly noticeable at the tails of the lobes (see the zoomed regions of the tail of the western lobe and the jet emission, displayed on the top right of each panel of Figure 5). By design, the dynamic range of the WSClean restored image is severely limited by the addition of the residual image, a standard post-processing step for CLEAN-based approaches. Indeed, when the latter is removed, the model image convolved with a restoring beam (bottom rows of Figures 45) reveals the captured faint and extended emission. Nonetheless, the image exhibits a large number of artefacts in the form of multi-resolution 2D Gaussian-like components, not to mention its physically-inconsistent negative components. 4 Calibration errors can be observed in the reconstructions of both algorithms. Firstly, these artefacts take the form of structures emanating from the hotspots in the model images of Faceted HyperSARA, and ringing artefacts around the hotspots 4 Further tuning WSClean by increasing the auto-masking threshold has shown to mitigate the number of the artefacts in the model image convolved with the restoring beam, at the expense of resolution, particularly noticeable at the jets and the ends of both lobes. This however has no impact on the quality of the restored image, dominated by the calibration artefacts left in the residual image.

in WSClean images. Similar structures were observed in reconstructions of Cyg A at similar frequencies [START_REF] Dabbech | [END_REF]) and were attributed to pointing errors. Secondly, structures around the core of the radio galaxy are observed in the low frequency-channel image recovered by Faceted HyperSARA (top row of Figure 4) and WSClean reconstructed images at both channels (bottom rows of Figures 45).

Residual images displayed in linear scale, obtained with both algorithms, are displayed over the associated images (top and middle rows of each figure). Note that the images associated with both methods are not obtained using the same measurement operator. Faceted HyperSARA adopts a naturally weighted de-gridding operator. In contrast, the CLEAN-based approach adopts a Briggs-weighted degridding measurement operator, thus emphasizing the high spatial frequency content. The resulting residual images are therefore not comparable in terms of patterns and structure. On the one hand, ringing structures centered at the hotspots are noticed in Faceted Hy-perSARA residual images for the two selected channels. On the other hand, the residual images of WSClean present high signal at the pixel positions of the hotspots. One can also notice ringing structures of somewhat smaller scale than those observed in Faceted HyperSARA residual images due to the different data weighting schemes adopted in both algorithms. The standard deviations of both residual images, computed after normalization with the flux of the associated restoring beams and reported in the captions of Figures 45, indicate the higher fidelity to data of Faceted HyperSARA.

Examining the reconstructed image cubes, available online in both FITs and GIF formats (Thouvenin et al. 2022b), confirms the analysis reported above. As a single effective channel per spectral window is considered, the different emissions in Cyg A recovered by our proposed approach exhibit smooth spectra in general, confirming a typical power-law decay. One spectral discontinuity is observed the effective channel 𝜈 17 = 6.039 GHz, that is at the intersection between the two effective frequency ranges [ 𝜈 1 , 𝜈 16 ] = [4, 039, 5.959] GHz and [ 𝜈 17 , 𝜈 29 ] = [6.039, 7.959] GHz. As already reported in Section 6.2.1, this behaviour is most likely induced by the different calibration errors affecting the independently acquired observations in these two frequency ranges. On a final note, the examination of the spectra at different pixel positions in WSClean restored cube is hampered by the channel-varying errors present in the associated residual cube.

Computing cost

Faceted HyperSARA exploited 90 CPU cores for about 32 hours as follows; 1 master CPU core, 60 CPU cores handling the spatial facets and 29 CPU cores handling the data-fidelity terms, and 2903 cpu hours in total. Further details on the computational and memory cost are reported in Table 3. WSClean was run with 29 CPU cores for about 17 hours, yielding a total computing cost of about 493 cpu hours (assuming all the cores are active the whole time). Though the numbers are sensitive to the resources' allocation, the total cpu time of the C++ RI imager remains significantly lower than that required by the MATLAB implementation of Faceted HyperSARA.

CONCLUSION AND FUTURE WORK

In this paper, new functionality of Faceted HyperSARA, a wideband intensity imaging approach for RI proposed in the companion paper [START_REF] Thouvenin | Parallel faceted imaging in radio interferometry via proximal splitting (Faceted HyperSARA): I. Algorithm and simulations[END_REF], has been introduced for enhanced robustness and scalability when accommodating large scale RI data. On the one hand, available direction-dependent models, including 3. Experiment 2: Computing cost of Cyg A imaging at the spectral resolution 𝛿𝜈 = 128 MHz and spatial resolution 𝛿 𝑥 = 0.0424 ′′ . Results are reported for Faceted HyperSARA (𝐶 = 1, 𝑄 = 60) in terms of number of CPU cores, number of PDFB iterations, runtime per iteration and total runtime (run pips in seconds with the associated standard deviation, run in hours), CPU time per iteration and total CPU time (cpu pips in seconds with the associated standard deviation, cpu in hours). The total memory allocated for the computation of the holographic matrices and the the memory needed to store them are reported.

the 𝑤-terms, are supported in the measurement operator to enhance its accuracy. In addition, an adaptive strategy for on the fly estimation of the effective noise level can be activated, in case of unknown noise levels due to the imperfect calibration, ensuring the robustness of the imaging algorithm in terms of reconstruction quality. On the other hand, significant reduction of the dimension of the data is achieved via a visibility gridding-based scheme.

A fully documented MATLAB library (available on Github ht tps://basp-group.github.io/Faceted-HyperSARA/ (Thouvenin et al. 2022a)) was introduced as part of a wider effort towards a larger adoption of the algorithms of the SARA family by the community.

The enhanced algorithm and code are assessed on two real data experiments in comparison with SARA (Carrillo et al. 2012;Onose et al. 2017) and the wideband CLEAN-based method implemented in WSClean (Offringa & Smirnov 2017). Firstly, leveraging the power of a large scale High Performance Computing system, Faceted Hyper-SARA was validated on the reconstruction of a 15 GB image cube of Cyg A from 7.4 GB of VLA data. The obtained results are a practical proof of concept of the parallelization flexibility of Faceted Hyper-SARA, which is also shown to provide a significant improvement in the imaging quality with respect to SARA. Secondly, the image formation process of Cyg A was revisited considering an important reduction of the image dimension via a reduction in the spectral resolution, and of the data dimension via a visibility gridding-based scheme. In this setting, significant reductions in memory requirements and computing cost are observed. More interestingly, Faceted HyperSARA has shown to be effective in recovering a higher quality wideband image, compared with the state-of-the-art CLEAN-based approach of the WSClean software.

Having addressed the computational bottlenecks raised by both the volume of the data and the size of the image, future work could contemplate the use of a faceted Fourier transform to improve the data and image locality in the proposed algorithm [START_REF] Wortmann | Streaming distributed Fourier transform for large-scale interferometry imaging[END_REF]. Other perspectives include the extension of the approach to polarization imaging and the development of a production C++ version of Faceted HyperSARA, building on the existing C++ version of HyperSARA (see the Puri-Psi webpage [START_REF] Wiaux | Puri-Psi library: Parallel Proximal Scalable Imaging[END_REF]). Ultimately, beyond image estimation, the faceted approach should also be integrated into a joint calibration and imaging approach (Repetti et al. 2017) and combined with the recently proposed methods for uncertainty quantification by convex optimization (e.g. [START_REF] Repetti | 2018 IEEE EUSIPCO[END_REF][START_REF] Repetti | [END_REF]Abdulaziz et al. 2019a).

  (i) Data extraction from a measurement set. A dedicated Python script enables the extraction of the RI data from one or multiple MS Tables into a collection of .mat. The script relies on Pyxis, standing for Python Extensions for Interferometry Scripting, which is part of the MeqTrees software package (Noordam, J. E. & Smirnov, O. M. 2010). The extracted data files are then systematically fed as input to our MATLAB imaging pipeline;

Figure 1 .

 1 Figure1. Proposed library structure and illustration of the Faceted HyperSARA imaging pipeline. The core module (Faceted-Hypersara) appears in red, whereas auxiliary submodules (RI-measurement-operator and SARA-dictionary) appear in blue and yellow, respectively.

Figure 2 .

 2 Figure 2. Experiment 1: Imaging results of Cyg A at the spectral resolution 𝛿𝜈 = 8 MHz and spatial resolution 𝛿 𝑥 = 0.06 ′′ . Reconstructions of channel 𝜈 102 = 4.787 GHz. Top and bottom rows: estimated model images of the proposed Faceted HyperSARA with (𝑄 = 60, 𝐶 = 16) and SARA, respectively, in units of Jy/pixel and displayed in log 10 scale. They are overlaid with zooms on several regions of the radio galaxy, highlighted in dashed circles to enhance the visual comparison between the two algorithms, and the associated residual images displayed in linear scale. The respective standard deviations of Faceted HyperSARA and SARA, after a normalization by the flux of the restoring beam, are 0.021 mJy/pixel and 0.020 mJy/pixel. Full image cubes are available online (Thouvenin et al. 2022b).

Figure 3 .

 3 Figure 3. Experiment 1: Imaging results of Cyg A at the spectral resolution 𝛿𝜈 = 8 MHz and spatial resolution 𝛿 𝑥 = 0.06 ′′ . Reconstructions of channel 𝜈 459 = 7.851 GHz. Top and bottom rows: estimated model images of the proposed Faceted HyperSARA with (𝑄 = 60, 𝐶 = 16) and SARA, respectively, in units of Jy/pixel and displayed in log 10 scale. They are overlaid with zooms on several regions of the radio galaxy, highlighted in dashed circles to enhance the visual comparison between the two algorithms, and the associated residual images displayed in linear scale. The respective standard deviations of Faceted HyperSARA and SARA, after a normalization by the flux of the restoring beam, are 0.026 mJy/pixel and 0.027 mJy/pixel. Full image cubes are available online (Thouvenin et al. 2022b).

Figure 4 .Figure 5 .

 45 Figure 4. Experiment Imaging results of Cyg A at the spectral resolution 𝛿𝜈 = 128 MHz and spatial resolution 𝛿 𝑥 = 0.0424 ′′ . Reconstructions of channel 𝜈 10 = 5.191 GHz. From top to bottom: estimated model image of Faceted HyperSARA, the restored image of the wideband CLEAN in WSClean, and the associated model image convolved with a restoring beam, obtained after the subtraction of the residual image, all in units of Jy/pixel and displayed in log 10 scale (for visualisation purposes, WSClean images are normalized with the flux of the associated restoring beam). The images are overlaid with zooms on several regions of the radio galaxy, highlighted in dashed circles. Residual images of Faceted HyperSARA and WSClean, displayed in linear scale, are overlaying the respective model image (top row) and restored image (middle row). Their respective standard deviations, after a normalization by the flux of the associated restoring beams, are 0.006 mJy/pixel and 0.011 mJy/pixel. Full image cubes are available online (Thouvenin et al. 2022b). MNRAS 000, 1-15 (2022)

Table 1 .

 1 Problem dimensions and mathematical notations

	Problem dimensions
	𝑁	image size
	𝐾	size of the 2D oversampled Fourier space
	𝐵	number of data blocks per channel
	𝐿	number of spectral channels
	𝑀	number of visibilities
	𝐶	number of spectral sub-cubes
	𝑄	number of spatial facets
	Mathematical notations
	Z †	adjoint (conjugate transpose) of the matrix Z
	𝜄 𝐶	indicator function of the set 𝐶

  real data considered in this paper are part of wideband VLA observations of the celebrated radio galaxy Cyg A, acquired in the frequency range 2-18 GHz over the years 2015-2016 (courtesy of R. Perley). A number of 𝐿 = 480 channels in the C band were used. The selected frequencies span the range [𝜈 1 , 𝜈 480 ] = [3.979, 8.019] GHz, characterized by a channel-width 𝛿𝜈 = 8 MHz and a total bandwidth of 4.04 GHz. The observations' phase center is located at the position RA = 19h 59mn 28.356s (𝐽2000) and DEC = +40 • 44 ′ 2.07 ′′ . The dataset was acquired at four instances, corresponding to the frequency ranges [𝜈 1 , 𝜈 256 ] = [3.979, 6.019] GHz and [𝜈 257 , 𝜈 480 ] = [5.979, 8.019] GHz, and VLA configurations A and C. The wideband data cube consists of 30 spectral windows, each composed of 16 channels spanning 128 MHz. The dataset contains approximately 10 6 complex visibilities per channel (about 8 × 10 5 and 2 × 10 5

Table 2 .

 2 . More Experiment 1: Computing cost of Cyg A imaging at the spectral resolution 𝛿𝜈 = 8 MHz and spatial resolution 𝛿 𝑥 = 0.06 ′′ . Results are reported for Faceted HyperSARA (with 𝐶 = 16, 𝑄 = 60) and SARA (i.e.

		Faceted HyperSARA SARA
	CPU cores (𝑁 𝑐 )	91	12
	Total CPU cores (=𝑁 𝑐 × 𝐶)	1456	5760
	average PDFB iterations	4598	2513
	run pips (s)	20.05 (±0.81)	6.25 (±0.56)
	run (h)	34.01	8.04
	cpu pips (s)	1824.55 (±73.71)	75.02 (±6.74)
	cpu (h)	37324	25141
	Memory allocated to compute the G 𝑐,𝑙,𝑏 matrices	3.36 TB	3.36 TB
	Memory occupied by the G 𝑐,𝑙,𝑏 matrices	2.5 TB	2.5 TB

Table

		Faceted HyperSARA
	CPU cores	90
	PDFB iterations	2090
	run pips (s)	55.57 (±3.70)
	run (h)	32.26
	cpu pips (s)	5001.3 (±333)
	cpu (h)	2903
	Memory allocated to compute 𝚲 𝑙,𝑏 H 𝑙,𝑏 matrices	203 GB
	Memory occupied by 𝚲 𝑙,𝑏 H 𝑙,𝑏 matrices	108 GB

). MNRAS 000,

1-15 (2022)

© 2022 The Authors

MNRAS 000,1-15 (2022)

All figures in this article are displayed using the color scheme cubehelix proposed inGreen (2011).

https://epsrc.ukri.org/research/facilities/hpc/tier2/ MNRAS 000,1-15 (2022)

ACKNOWLEDGEMENTS

The first two authors contributed equally to this work. The authors warmly thank R.A. Perley (NRAO, USA) for providing the VLA observations of Cyg A. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. We also thank H.L. Bester (SARAO, South Africa) for useful discussions on data dimensionality reduction and A. Wilber (Heriot-Watt University, UK) for helpful feedback on the MATLAB library. This work was supported by EPSRC, grants EP/M011089/1, EP/M008843/1, EP/M019306/1, EP/T028270/1, EP/T028351/1, and ST/W000970/1 the Swiss-South Africa Joint Research Program (IZLSZ2_170863/1), and used the Cirrus UK National Tier-2 HPC Service at EPCC (http://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC (EP/P020267/1).

DATA AVAILABILITY

The images underlying this article are available from the Heriot-Watt University research portal, with the digital object identifier 10.17861/f127fd8d-9394-4d74-abf0-7af35bd7192e. MAT-LAB codes associated with this article are made available online on the Faceted HyperSARA page https://github.com/basp-g roup/Faceted-HyperSARA of the Puri-Psi library Wiaux (2021).

APPENDIX A: WSCLEAN COMMAND USED IN EXPERIMENT 2

The command used in WSClean to obtain the results of Experiment 2, displayed in Figures 4 and5, is provided below.

wsclean -join-channels -channels-out 29 \ -spws 0, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,\ 15,16,17,18,19,20,21,22,23,24,25,26,27,29 \ -fit-spectral-pol 3 -reorder -scale 0.0424arcsec \ -weight briggs -0.25 -data-column DATA \ -pol I -size 5120 3072 -niter 1000000 \ -nwlayers 1 -field 2 -mgain 0.8 -gain 0.08 \ -auto-mask 3 -auto-threshold 1 \ -no-small-inversion -multiscale \ -name CYGA-EXP2 CYGA.MS This paper has been typeset from a T E X/L A T E X file prepared by the author.


Patrick Ballard

FORMULATION AND WELL-POSEDNESS OF UNILATERAL MULTIBODY DYNAMICS

The classical theory of rigid bodies systems undergoing perfect bilateral constraints has received firm mathematical foundations for a long time. We extend it to the case where there are perfect unilateral constraints also. The formulation follows the line initiated by Schatzman (1978) and Moreau (1983). We give a parametrization-free formulation and try to identify the most general form of impact constitutive equation which is compatible with well-posedness. Then, well-posedness is proved under the assumption of analyticity of the data, since it is known, thanks to Bressan (1960) and Schatzman (1978), that uniqueness does not hold in the C 00 setting.

1.

Background: multibody systems with perfect bilateral constraints

Consider a finite collection of rigid bodies undergoing perfect bilateral constraints. It is classically associated with a manifold Q ( called the configuration manifold of the system) of finite dimension d {called the number of degrees-of-freedom of the system).

A motion of the system is a curve on Q, that is a mapping q(t) from a real interval I into Q. When a motion is smooth enough, there exists a time-derivative q(t) in tangent space T q ( t ) Q also called a (generalized) velocity of the system. We shall also denote it by (q(t), q(t)). An arbi trary element v ( also denoted by ( q, v), q being the basepoint of v) of the tangent bundle TQ is called a {virtual) state of the system.

Given the mass distribution in each of the solids, one classically ob tains the kinetic energy K which appears to define a positive definite quadratic form on each tangent space T q Q of the configuration man ifold, endowing it with a Riemannian structure. The induced scalar product and norm will be denoted by (•, • ) q and 11 • ll q , whereas the du ality product between tangent space T q Q and cotangent space T;Q will be denoted by(•, •) q -With this notation, we have:

( 1 2 K q,v) = 2 llvll q •
The modelling of forces makes use of duality. Internal and external forces acting on the system in the configuration q are represented by a linear form f E T;Q. Given an arbitrary virtual velocity (q, v) of the system, the real number(/, v} q is the virtual power of internal and external forces for the virtual velocity (q, v). The cotangent vector f is usually allowed to depend on the current state (q, v) and also on time t.

Hence, we are given a mapping/ : TQ x JR + ➔ T*Q satisfying:

\/(q,v)E TQ, VtEJR + , II,if( q,v ; t ) ) = q, ( 1
)
( IIQ is the natural projection on T*Q), called the force mapping.

Since such a structure appears systematically in the modelling of the dynamics of rigid bodies systems with perfect bilateral constraints, we are led to make the following definition.

Definition 1 A simple discrete mechanical system is a pair (Q, /) whe re:

■ Q is a finite-dimensional Riemannian manifold called the configu ration manifold.

■ f: TQ x JR + ➔ T*Q is a mapping satisfying (1), called the force mapping.

Now, consider an arbitrary smooth motion q ( t) of the system. The power of inertial forces is by definition the time derivative of the kinetic energy:

: t K ( q, q) = : t � llq( t )ll! ct > = (�q( t ) , q( t )) = / b �q( t ) , q( t )) , q(t) \ q(t)
where D/dt denotes the covariant derivative along the curve q(t) associ ated with the Levi-Civita connection, and b the isomorphism from TQ onto T*Q canonically associated with the Riemannian metric. For any virtual velocity ( q, v) of the system, the real number (b Dq / dt, v) q is nat urally called the virtual power of inertial forces in the virtual velocity (q, v). Given any local parametrization of the system ( local chart), it admits the well-known representation: b�q = (:

t : i K ( q, q )-� i K ( q, q )) dq i .
The fundamental principle of classical dynamics asserts that the vir tual power of inertial forces should equal the virtual power of external and internal forces, for any virtual velocity. We obtain immediately the equation of motion:

'tit, b�cj(t) = f(q(t), cj(t);t).
Next, we are given an initial state ( qo, vo) E TQ. Then, the evolution problem associated with the dynamics of rigid bodies system with perfect bilateral constraints is the Cauchy problem:

Problem I. Find T > 0 and q E 0 2 {[0, T[ ; Q) such that: • (q(0), cj(0)) = (qo, vo), • 't/tE[0, T[, b�cj(t) = f (q(t ), cj(t), t).
Well-posedness is ensured by the following.

Theorem 2 (Cauchy) The Riemannian configuration manifold is as sumed to be of class 0 2 and the mapping f:

TQ x JR + -+ T*Q is of class 0 1 .
Then, there exists a unique maximal solution for problem I.

More precisely, theorem 2 states that there exists Tm > 0 {T m E JR + U { +oo}) and Qm E 0 2 ([0, T m [, Q) being a solution of problem I such that any other solution of problem I is a restriction of Qm• Of course, we expect that Tm = +oo, in which case the dynamics is said to be eternal. This situation cannot be taken for granted, in general. In the usual cases encountered in mechanics, eternal dynamics is ensured by the following general sufficient condition.

Theorem 3 The configuration manifold Q is assumed to be a complete Riemannian manifold. The effort mapping f is supposed to admit the following estimate: ll f(q, v ; t) II; � l(t) (1 +d (q, q o ) + llvll q ), for all ( q, v) E TQ and almost all t E [to, +oo[, where d( •, •) is the Riemannian distance and l(t), a (necessarily nonnegative) function of L}oc ( JR ; JR) .

Then, the dynamics is eternal: Tm = +oo.

Formulation of unilateral multibody dynamics

The consideration of elementary examples shows that the dynamics of rigid bodies systems can lead to predictions of the motion in which some bodies of the system ov erlap in the real world. Of course, this cannot be allowed. Hence, very often, one has to add the statement of non-penetration conditions to a simple discrete mechanical system. This is a simple occurrence of a unilateral constraint. Here, we shall recall briefly the general formulation of the resulting evolution problem. This is due essentially to Schatzman {1978) and Moreau {1983).

Consider a simple discrete mechanical system according to defini tion 1. A unilateral constraint is a restriction on the admissible of the system which is expressed by means of a finite number n of smooth real-valued functions 'Pi on the configuration manifold Q, so that the set of all admissible configurations A is given by: A= {q E Q; 'v'i E {1,2, ... ,n}, '-Pi(q) � 0}.

(

) 2
The set of all active constraints in the admissible configuration q E A is defined by: J(q ) = {iE{l,2, ... ,n}; 'Pi(q) =O} .

The following hypothesis is usual in this framework. The functions 'Pi are assumed to be functionally independent in the sense that, for all q EA, the dipi(q) (i E J(q )) are linearly independent in T*Q.

The cone of admissible right velocities in the configuration q is defined by: V(q) = {v E T q Q; 'v'i E J(q ), { dipi(q), v} q � O}, and we denote by N*(q) its polar cone for the duality (T q Q, T;Q):

the polar cone of V(q) for the euclidean structure of T q Q being N(q) = #(N*(q)) (# = b-1 ).

The equation of motion has to be corrected with some reaction efforts R:

D • + L q -f( •+ . ) R
"Ttq,q , t + .

(3)

We require the unilateral constraint to be perfect, that is, the reaction effort R has to take values in the normal cone N*(q). It is well known that the dynamics of such a system involves some time-discontinuities of velocity called impacts. Therefore, we cannot expect that the equation of motion (3) should be satisfied in a classical sense, but rather in a distributional sense. Actually, it is usual to require that R should be a vector-valued mesasure rather than a general distribution. Hence, we define the class of motions MM A(I; Q) (motions with measure acceler ation) to be the set of all absolutely continuous motions q(t) from a real interval I to Q admitting a right velocity q + (t) at every instant t of I and such that the function q + (t) has locally bounded variation over I. Bounded variation is classically only for functions taking values in a normed vector space. However, for any absolutely continuous curve q(t) on a lliemannian manifold, parallel translation along q(t) classically provides intrinsic identification of the tangent spaces at different points of the curve and so, the definitions can easily be carried over to this case (for a precise mathematical setting, see [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF]). Any mo tion q E MMA(I; Q) admits a left and right velocity, qand q+, in the classical sense at any instant. Moreover, any motion q E MMA(I; Q) is intrinsically associated with the covariant Stieltjes measure Dq + of its right velocity q + . The equation of motion takes the form:

where dt denotes the Lebesgue measure. It remains to express the con dition that R should take values in N*(q) when Risa measure.

Convention. We shall write: RE -N*(q(t)) to mean: there exist n nonpositive real measures Ai such that: n R = L Ai dc,oi(q(t)) and i=l 'v'iE{l,2, ... ,n } , SuppAiC{t; c,oi(q(t))=0} . ( 4) With this convention, the final form of the equation of motion is: R = �Dq + -f(q(t),q + (t);t) dt E -N *(q(t))

(5)

It is classical that this equation of motion has to be completed with some constitutive impact equation describing the outcome of any impact occurrence. We shall assume an impact constitutive equation of general form:

(6) To ensure compatibility with the equation of motion ( 5), the mapping F should satisfy:

'v'q EA, 'v'v-E -V(q), F (q, v-) E V(q), F ( q, v-) -v-E -N(q). ( 7 )
Moreover, we add the assumption that the kinetic energy of the system cannot increase during an impact:

'v' q EA, 'v' v-E -V(q), (8)
Now, we formulate the evolution problem associated with the dy namics of rigid bodies systems with perfect bilateral and unilateral con straints. The initial condition is assumed to be compatible with the realization of the constraint: v 0 V(qo).

Problem II. Find T > 0 and q E MMA([0, T[; Q) such that:

• (q(0), q + (0)) = (qo, vo),

• 'v't E [ 0, T[ , q(t) EA, • R � �D q + -f (q(t),q + (t) ; t) dt E -N* (q(t)), • 'v'tE)0,T[ , q + (t)=.r(q(t),q-(t)).
The equation of motion is understood in the sense of convention ( 4), and the impact constitutive equation is supposed to fulfill requirements ( 7) and ( 8).

Well-posedness of unilateral multibody dynamics

To study existence and uniqueness of solution for problem II, we need to state some regularity hypotheses on the data. It was pointed out by [START_REF] Bressan | Incompatibilita dei teoremi di esistenza e di unicita del moto per un tipo molto comune e regolare di sistemi meccanici[END_REF] and [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF], through a striking example, that under the strong hypotheses that the data are all of class C 00 there may exist non-unique solutions for problem II. It was [START_REF] Percivale | Uniqueness in the elastic bounce problem, I[END_REF] and [START_REF] Schatzman | Uniqueness and continuous dependence on data for one dimen sional impact problems[END_REF] who were the first to notice that this trouble could be overcome by assuming that the data are analytic. But their analysis was limited to one-degree-of-freedom systems. The general case is treated in [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF], [START_REF] Ballard | Formulation and well-posedness of the dynamics of rigid bodies systems with perfect unilate ral constraints[END_REF] and is now briefly recalled.

From now on, we assume the following Regularity hypothesis. The Riemannian configuration manifold Q, the functions 'P i and the mapping f: TQ x JR + ➔ T * Q are analytic.

Proposition 4 Let qo EA and vo E V(q 0 ). Then, there exist Ta > 0, an analytic curve Qa: [ 0, Ta[ ➔ Q and n analytic functions Aa i : [ to, Ta[ ➔ IR such that:

• (qa(0),qf(0)) = (qo,vo),

• Vt E [0, Ta[, � � Qa(t) = f(qa(t), Qa(t); t) + t Aai(t) dcpi{ Q a{t)), i=l • �� :, [ 1 0, 2 T�[_ , .
Aa i{t ) :'.S 0, 'Pi(qa(t)) :'.S 0, Aa i{t ) 'Pi(Qa(t)) = 0 .

V Z 1 1

1 n, Moreover, the solution of this evolution problem is unique the sense that any other analytic solution (T, q, .\ 1 , ... , A n) is either a restriction or an analytic extension of ( Ta, Q a, Aa 1 , ... , Aan)- An early proof of proposition 4 can be found in L6tstedt {1 982).

Corollary 5 There exists an analytic solution (T a , q 0 ) for problem II.

Naturally, the analytic solution furnished by corollary 5 will cease to exist at the first instant of impact. This is the reason why we have considered the wider class MMA which contains motions which are not differentiable in the classical sense. Considering motions in MMA will allow us to extend the solution beyond the first instant of impact. But, we mu st ensure that admitting the wider class of solutions MMA will not introduce parasitic solutions. This is the aim of the following theorem.

Theorem 8

The configurotion manifold Q is assumed to be a complete Riemannian manifold and the mapping f is supposed to admit the fol lowing estimate: ll f ( q ,v ; t )II; � l (t ) (1 + d( q,qo) + llvll q ), for all ( q , v ) E TQ and almost all t E [O, +oo[, where d( •, •) is the Riemannian distance and l (t ) , a (necessarily nonnegative) function of L} oc: ( JR; R).

Then, the dynamics is eternal, that is, the maximal solution for prob lem II is defined on [O, +oo[.

Theorem 6

 6 Let ( Ta , Q a) be the solution for problem II furnished by corollary 5, and (T, q) be an arbitrary solution for problem II. Then, there exists a real number To {0 < To :'.S min{T a , T}) such that: Ql[O,To[ = Qal[O,To[ • In other terms, there is local uniqueness for problem II. Local uniqueness is the difficult part in the study of well-posedness of problem IL Corollary 7 There exists a unique maximal solution for problem II. It was noticed above that the analytical solution for problem II fur nished by corollary 5 fails to exist at the first instant of impact. To overcome this fact, we have proved that local uniqueness still holds in the wider class of motion MMA which allows impacts. However, the maximal solution for problem II may stop to exist at finite time for un physical reasons. In other terms, we still do not know if the class MMA is wide enough. Actually, it is wide enough as shown by the following theorem which should be brought aside theorem 3.


Tom Archibald

Richard T W Arthur

Giovanni Ferraro

Jeremy Gray

Douglas Jesseph

A question of fundamental methodology: Reply to Katz and his co-authors

In the last ten years or so, Mikhail G. Katz has published with various co-authors more than 60 papers in the history and philosophy of infinitesimal mathematics. 1 These papers are part of a sustained programme to rewrite the history of the calculus from a point of view where non-standard analysis is assumed to give the correct framework for historical interpretation. As described in Bair et al. (2013a, p. 887) the aim is "a re-evaluation of the history and philosophy of mathematics, analysing the shortcomings of received views, and shedding new light on the deleterious effect of the latter on the philosophy, the practice, and the applications of mathematics". Here the "received view", which has had such a "deleterious" effect, is presented as that of most historians of mathematics in the domain under consideration. Scholars who have not adopted the new non-standard framework are vilified as being "against" infinitesimals. 2 Needless to say, this is a gross oversimplification. Whether infinitesimal methods are advocated by the historical actors has not been disputed by any of the historians of mathematics these authors criticize, whose concern is to articulate how these methods were understood by their proponents.

What is at stake here is a question of fundamental methodology in the history of mathematics. For the structure of the papers by Katz and co-authors is always the same: it consists mainly of intemperate criticism of the supposed "received views" based on chosen extracts from the said historians of mathematics and a list of chosen extracts from the primary sources, most of the time decontextualized and presented as testifying that the nonstandard interpretation is the only one able to provide a faithful account of these sources. Many examples could be given, involving different historians of mathematics, and different pieces of the history of analysis. But what is relevant is not who these scholars are or which aspects of this history are concerned, nor whether the theses criticized are in fact well argued for, plausible, or questionable. The delicate point is rather with the form of these criticisms. They do not consist in reasoned historical argument, but rather in piecemeal confutation of isolated quotations of their opponents taken out of context, whose sole aim is to show them as enemies of the group's understanding of infinitesimals. The aim is not to contribute to historical clarification, nor even to rectify possible shortcomings or errors, but to put various scholars in the pillory-with accompanying abusive epithets-for not enthusiastically recognising the pervasive presence of non-standard analysis in the history of mathematics.

For a long time, we have refrained from sending "replies" underlining the various shortcomings of these papers, because we thought that the debate should mainly consist in rational arguments and studies of implications of newly discovered documents. But, over the years, it became clearer and clearer that our interlocutors do not care much about rational discussion and scientific dialogue from different perspectives, but seek rather to disparage Preprint. Please quote the published version: Tom Archibald, Richard T. W. Arthur, Giovanni Ferraro, Jeremy Gray, Douglas Jesseph, Jesper Lützen, Marco Panza, David Rabouin and Gert Schubring, "A Question of Fundamental Methodology: Reply to Mikhail Katz and His Coauthors", The Mathematical Intelligencer, 2022, online first: https://doi.org/10.1007/s00283-022-10217-7 their alleged enemies, meanwhile changing their views in an inconsistent way in order to accommodate the various objections they receive. Our main concern here is to draw attention to this strategy, which we consider unacceptable in a properly rational debate. The latest example of that approach is provided by the paper recently published in the Mathematical Intelligencer, "Two-Track Depictions of Leibniz's Fictions" [START_REF] Katz | Two-Track Depictions of Leibniz's Fictions[END_REF]. Rather than attempting to counter it point by point, which would involve us argument and detailed historical contextualization (and also probably encourage more of the same scatter-gun responses from Katz and his collaborators), we decided to focus instead only on its main line of argument concerning the so-called "A" and "B" tracks, a central theme of Katz and co-authors since 2012. We will limit ourselves to demonstrating that these designations of theirs, and related theses such as the contradictory nature of infinitesimals in Leibniz, have undergone a radical change of meaning between 2012 and 2021, finally attaining a full-blown inconsistency.

In the first series of papers published by M. Katz and al., "A" and "B" were meant to describe "separate methodologies in Leibniz", one based on exhaustion arguments à la Archimedes, and the second relying "more directly on infinitesimals" à la Bernoulli (2013b, p. 575). The two methodologies were presented as corresponding to two pictures of the continuum: "• an A-continuum (for Archimedes), and

• a B-continuum (for Johann Bernoulli, following Leibniz).
The former is a 'thin' continuum, exemplified by what are called today the real numbers; the latter is a 'thick' continuum incorporating infinitesimals." (2013b, p. 598) As can be seen from this quotation, it was not a matter of conflicting interpretations between scholars, but of methodologies existing in Leibniz. The main claim was that Leibniz shared with Bernoulli a certain view of the continuum as consisting of infinitesimal numbers in addition to ordinary (or "assignable") numbers. We may note in passing that this already involves anachronism at odds with a properly historical approach. For Leibniz did not conceive of numbers as constituting a continuum, nor did he allow infinite sets (infinite wholes, in his terminology). In the case he discussed with Bernoulli in 1698-1699, the question was rather about whether there exists an infinitieth term in an infinite series, which in the case of a decreasing series would stand for an infinitesimal quantity. Bernoulli insisted that there would indeed be such an infinitieth (although not necessarily last) term, thus entailing the existence of an infinitesimal. According to Katz et al. this entails a conception of infinite series as consisting of an infinite sequence of standard numbers followed by an infinitesimal part. Bernoulli was not very clear about what happens after the infinitieth term, but Katz et al. (2013b, p. 599) explained how this can be conceived by giving a construction using nonstandard analysis. In that paper, as can be seen in the summary above, Leibniz and Bernoulli were presented as being in strong agreement. Infinitesimals were conceived as numbers, and the view of numerical sequences described by nonstandard analysis was presented as compatible with the philosophers' views. The assertion that infinitesimals were numbers was even presented as a feature of crucial importance: "the addition or concatenation of infinitesimals (of the same dimension) is no more difficult to conceive of than adding or concatenating finite magnitudes. This is especially important, because it allows one to represent infinitesimals by means of numbers and so apply arithmetic operations to them. This is the fundamental difference between the infinitary methods of Archimedes (and later Cavalieri) and the infinitary methods of Leibniz and his followers." (2013b, p. 575; our emphasis).

Many objections have been raised to this general picture, the most obvious one being that it testifies to a non-historical and decontextualized methodology. Indeed, we possess the documents in which Leibniz discusses these questions with Bernoulli and everybody can see, by just reading them, that the two mathematicians strongly disagreed on all the issues presented above. There is no way that one can claim that Bernoulli defended a certain picture of the continuum "following Leibniz". Accordingly, there is no B-methodology stricto sensu in Leibniz. Leibniz's main argument is that it is not possible to treat infinitesimals as existing entities because this amounts to the introduction of an infinite number, which he takes to be a contradictory notion. 3After receiving this objection, Katz and al. completely changed their position, but without acknowledging this change, and as if it did not ruin their previous argument. In the paper published in your journal, which is supposed to give a survey of a long-lasting debate, "A" and "B" are no longer presented as a pair of methodologies in Leibniz, but as positions endorsed by commentators to understand the term "fiction" in Leibniz. "B" does not refer to Bernoulli any longer, but to the views endorsed by a generic interpreter named "Bob" (alias Katz and his co-authors). And here is the funny part: Bob now insists that Leibniz disagreed with Bernoulli and that this disagreement is based on the fact that an infinite number is a contradictory notion (2021, p. 1)! In order to "save" his position, Bob, who cannot be Bernoulli anymore, introduces a new distinction: one should not confuse, he says, infinite (and infinitesimal) number and infinite (and infinitesimal) magnitude -with infinite magnitude exemplified by the idea of a linea infinita terminata. A crucial point is that the second pair are supposed to be non-contradictory objects. According to Bob, this is now supposed to be the "proper" meaning of fiction in Leibniz.

We will not enter in the discussion of this new claim here, 4 but just content ourselves with pointing to another paper published by Bob, alias Katz and Sherry, in Studia Leibnitiana (2012). The aim of the paper was to explain the "proper meaning" of the term fiction in Leibniz in the direct line of the argument developed in the paper published in your journal. But here is how it goes: "Infinitesimals and imaginary roots are not, therefore, objects of mathematical science, and their fictional status rests ultimately upon this absence of objectivity. These Preprint. Please quote the published version: Tom Archibald, Richard T. W. Arthur, Giovanni Ferraro, Jeremy Gray, Douglas Jesseph, Jesper Lützen, Marco Panza, David Rabouin and Gert Schubring, "A Question of Fundamental Methodology: Reply to Mikhail Katz and His Coauthors", The Mathematical Intelligencer, 2022, online first: https://doi.org/10.1007/s00283-022-10217-7 concepts, because they contain a contradiction, prevent us from imagining objects in accordance with their definitions" (2012, p. 179; our emphasis) At the time, this went along with the following claim: "because they involve contradictions, Leibniz rejects infinitesimals, along with negatives and imaginaries as objects of mathematical science" (2012, p. 190). Moreover, the conclusion of the paper made it clear that the contradiction related to infinitesimals was for Bob of the same type as the one based on infinite cardinalities, in concordance with the fact that at this time Bob thought that infinitesimals were truly represented by numbers. Nonstandard techniques of the kind invented by Robinson were presented as a way to extend the realm of fictions related to infinite cardinalities by dealing with them in a consistent logical theory (how the two claims were compatible is not clear, but that is not the issue here). 5 But if Leibnizian infinitesimals have nothing to do with infinite cardinalities, if they are non-contradictory objects, as claimed by Bob2021, what he is criticizing is not the Ainterpretation, but the B-interpretation in the sense of Bob2012. We simply do not understand what the "B" position is here, as we do not understand whether or not Leibniz was supposed to have shared the views of Bernoulli (as claimed by Bob2013 but refuted by Bob2021). What we suspect is that "A" and "B" are mere artefacts invented to hide inconsistent positions and we recommend that Bob 2021 engage in a polemic with his own self from 2012-2013 directly, rather than confronting historians of mathematics with six papers of this sort a year.

Finally, we note that the interest of Katz et al. in infinitesimals is not confined to reaching an understanding of the ideas of Leibniz and Johann Bernoulli, but extends to suggesting a tradition that includes Euler and reaches to Cauchy, before being shut down by Weierstrass, upon which the "leading historians" are supposed to have anachronistically relied. This extension rests on the confused views about Leibniz and Bernoulli described above, but it aims at promoting an interpretation of Cauchy's analysis as relying on the kinds of intuitions that can be made rigorous in a modern non-standard infinitesimal-based theory of analysis.

Here the A-track is defined as "the Weierstrassian approach (in the context of an Archimedean continuum)", and the B-track is "the approach with indivisibles and/or infinitesimals (in the context of what we will refer to as a Bernoullian continuum)". In order to protect themselves from what would be an absurd claim that Cauchy possessed the rigorous machinery of modern non-standard analysis, the authors supply the concept of hidden lemmas. These are provable results in non-standard analysis that have intuitive, if unrigorous, meanings; the suggestion is that Cauchy used these in his work in the B-track. This is to wrench Cauchy's ideas out of their historical context and misinterpret what he meant by 5 "Infinite sets have the property that they can be shown to be equinumerous with a proper part of themselves. This contradicts the principle that the whole is greater than the part, a principle to which Leibniz was firmly committed. Infinite aggregates would, therefore, be fictions from a Leibnizian perspective. Robinson's creation of a consistent theory of infinitesimals consists, roughly, in defining the hyper-reals by set theoretic means and showing that any first order statement that holds for a real function holds for the hyper-real natural extension of that function. (…). To be sure, the demonstration of this result requires resources unavailable to Leibniz. But for Robinson those resources occupy the same role as infinitesimals and imaginary roots did for Leibniz." (2012, pp. 191-192).

  Preprint. Please quote the published version: Tom Archibald, Richard T. W. Arthur, Giovanni Ferraro, Jeremy Gray, Douglas Jesseph, Jesper Lützen, Marco Panza, David Rabouin and Gert Schubring, "A Question of Fundamental Methodology: Reply to Mikhail Katz and His Coauthors", The Mathematical Intelligencer, 2022, online first: https://doi.org/10.1007/s00283-022-10217-7

A full list of these papers is given by Katz on the website https://u.math.biu.ac.il/~katzmik/infinitesimals.html (accessed May 31,

2022).2 A long list of scholars who have conducted their research without having adopted the non-standard analytical standpoint is proudly exhibited by Katz on his webpage referenced above ("List of Critics"). As can be seen, all that is important is for these "critics" to be identified and rebutted.

This had been a well-known fact amongst Leibniz scholars for some time, see for example Herbert Breger, "Das Kontinuum bei Leibniz", in A. Lamarra (ed.), L'infinito in Leibniz : Problemi e terminologia, Rome, Edizioni dell'Ataneo, 1990, pp 53-67, and in the same volume: George MacDonald Ross, "Are there real Infinitesimals in Leibniz's Metaphysics ?" (pp. 125-141); see also O. Bradley Bassler. "Leibniz on the Indefinite as Infinite", The Review ofMetaphysics, vol. 51, no.

4, 1998, pp. 849-74. 4 Arthur and Rabouin will dedicate a specific study to this, providing several sources in which Leibniz explicitly claimed that lineae infinitae terminatae are contradictory entities.

Preprint. Please quote the published version: Tom Archibald, Richard T. W. Arthur, Giovanni Ferraro, Jeremy Gray, Douglas Jesseph, Jesper Lützen, Marco Panza, David Rabouin and Gert Schubring, "A Question of Fundamental Methodology: Reply to Mikhail Katz and His Coauthors", The Mathematical Intelligencer, 2022, online first: https://doi.org/10.1007/s00283-022-10217-7 infinitesimal in a way that is somewhere between tendentious and incorrect, as Gilain, Lützen, Schubring, and others have shown. 6