Upcoming radio interferometers are aiming to image the sky at new levels of resolution and sensitivity, with wide-band image cubes reaching close to the Petabyte scale for SKA. Modern proximal optimization algorithms have shown a potential to significantly outperform CLEAN thanks to their ability to inject complex image models to regularize the inverse problem for image formation from visibility data. They were also shown to be scalable to large data volumes thanks to a splitting functionality enabling the decomposition of data into blocks, for parallel processing of block-specific data-fidelity terms of the objective function. In this work, the splitting functionality is further exploited to decompose the image cube into spatio-spectral facets, and enable parallel processing of facet-specific regularization terms in the objective. The resulting Faceted HyperSARA algorithm is implemented in MATLAB (code available on GitHub). Simulation results on synthetic image cubes confirm that faceting can provide a major increase in scalability at no cost in imaging quality. A proof-of-concept reconstruction of a 15 GB image of Cyg A from 7.4 GB of VLA data, utilizing 496 CPU cores on a HPC system for 68 hours, confirms both scalability and a quantum jump in imaging quality from CLEAN. Assuming slow spectral slope of Cyg A, we also demonstrate that Faceted HyperSARA can be combined with a dimensionality reduction technique, enabling utilizing only 31 CPU cores for 142 hours to form the Cyg A image from the same data, while preserving reconstruction quality. Cyg A reconstructed cubes are available online.

New generation of radio interferometers are envisaged to produce high quality, high dynamic range Stokes images of the observed sky from the corresponding under-sampled Fourier domain measurements. In practice, these measurements are contaminated by the instrumental and atmospheric effects that are well represented by Jones matrices, and are most often varying with observation direction and time. These effects, usually unknown, act as a limiting factor in achieving the required imaging performance and thus, their calibration is crucial. To address this issue, we develop a global algorithm, named Polca SARA, aiming to perform full polarization, direction-dependent calibration and sparse imaging by employing a non-convex optimization technique. In contrast with the existing approaches, the proposed method offers global convergence guarantees and flexibility to incorporate sophisticated priors to regularize the imaging as well as the calibration problem. Thus, we adapt a polarimetric imaging specific method, enforcing the physical polarization constraint along with a sparsity prior for the sought images. We perform extensive simulation studies of the proposed algorithm. While indicating the superior performance of polarization constraint based imaging, the obtained results also highlight the importance of calibrating for direction-dependent effects as well as for off-diagonal terms (denoting polarization leakage) in the associated Jones matrices, without inclusion of which the imaging quality deteriorates.

We propose a Bayesian uncertainty quantification method for large-scale imaging inverse problems. Our method applies to all Bayesian models that are log-concave, where maximum-a-posteriori (MAP) estimation is a convex optimization problem. The method is a framework to analyse the confidence in specific structures observed in MAP estimates (e.g., lesions in medical imaging, celestial sources in astronomical imaging), to enable using them as evidence to inform decisions and conclusions. Precisely, following Bayesian decision theory, we seek to assert the structures under scrutiny by performing a Bayesian hypothesis test that proceeds as follows: firstly, it postulates that the structures are not present in the true image, and then seeks to use the data and prior knowledge to reject this null hypothesis with high probability. Computing such tests for imaging problems is generally very difficult because of the high dimensionality involved. A main feature of this work is to leverage probability concentration phenomena and the underlying convex geometry to formulate the Bayesian hypothesis test as a convex problem, that we then efficiently solve by using scalable optimization algorithms. This allows scaling to high-resolution and high-sensitivity imaging problems that are computationally unaffordable for other Bayesian computation approaches. We illustrate our methodology, dubbed BUQO (Bayesian Uncertainty Quantification by Optimization), on a range of challenging Fourier imaging problems arising in astronomy and medicine.

We develop a novel algorithm for sparse Stokes parameters imaging in radio interferometry under the polarization constraint. The latter is a physical non-linear relation between the Stokes parameters, imposing that the polarization intensity is a lower bound on the total intensity. To solve the joint inverse Stokes imaging problem including this bound, we leverage epigraphical projection techniques in convex optimization and design a primal-dual method offering a highly flexible and parallelizable structure. In addition, we propose to regularize each Stokes parameter map through an average sparsity prior in the context of a reweighted analysis approach (SARA). The resulting approach is dubbed Polarized SARA. We demonstrate on simulated observations of M87 with the Event Horizon Telescope that imposing the polarization constraint leads to superior image quality. The results also confirm that the performance of the average sparsity prior surpasses the alternative state-of-the-art priors for polarimetric imaging.

Radio interferometric imaging aims to estimate an unknown sky intensity image from degraded observations, acquired through an antenna array. In the theoretical case of a perfectly calibrated array, it has been shown that solving the corresponding imaging problem by iterative algorithms based on convex optimization and compressive sensing theory can be competitive with classical algorithms such as CLEAN. However, in practice, antenna-based gains are unknown and have to be calibrated. Future radio telescopes, such as the SKA, aim at improving imaging resolution and sensitivity by orders of magnitude. At this precision level, the direction-dependency of the gains must be accounted for, and radio interferometric imaging can be understood as a blind deconvolution problem. In this context, the underlying minimization problem is non-convex, and adapted techniques have to be designed. In this work, leveraging recent developments in non-convex optimization, we propose the first joint calibration and imaging method in radio interferometry, with proven convergence guarantees. Our approach, based on a block-coordinate forward-backward algorithm, jointly accounts for visibilities and suitable priors on both the image and the direction-dependent effects (DDEs). As demonstrated in recent works, sparsity remains the prior of choice for the image, while DDEs are modelled as smooth functions of the sky, i.e. spatially band-limited. Finally, we show through simulations the efficiency of our method, for the reconstruction of both images of point sources and complex extended sources. MATLAB code is available on GitHub.

In the context of optical interferometry, only undersampled power spectrum and bispectrum data are accessible. It poses an ill-posed inverse problem for image recovery. Recently, a tri-linear model was proposed for monochromatic imaging, leading to an alternated minimization problem. In that work, only a positivity constraint was considered, and the problem was solved by an approximated Gauss-Seidel method. In this paper, we propose to improve the approach on three fundamental aspects. Firstly, we define the estimated image as a solution of a regularized minimization problem, promoting sparsity in a fixed dictionary using either an $\ell_1$ or a weighted-$\ell_1$ regularization term. Secondly, we solve the resultant non-convex minimization problem using a block-coordinate forward-backward algorithm. This algorithm is able to deal both with smooth and non-smooth functions, and benefits from convergence guarantees even in a non-convex context. Finally, we generalize our model and algorithm to the hyperspectral case, promoting a joint sparsity prior through an $\ell_{2,1}$ regularization term. We present simulation results, both for monochromatic and hyperspectral cases, to validate the proposed approach.

In the context of next generation radio telescopes, like the Square Kilometre Array, the efficient processing of large-scale datasets is extremely important. Convex optimisation tasks under the compressive sensing framework have recently emerged and provide both enhanced image reconstruction quality and scalability to increasingly larger data sets. We focus herein mainly on scalability and propose two new convex optimisation algorithmic structures able to solve the convex optimisation tasks arising in radio-interferometric imaging. They rely on proximal splitting and forward-backward iterations and can be seen, by analogy with the CLEAN major-minor cycle, as running sophisticated CLEAN-like iterations in parallel in multiple data, prior, and image spaces. Both methods support any convex regularisation function, in particular the well studied l1 priors promoting image sparsity in an adequate domain. Tailored for big-data, they employ parallel and distributed computations to achieve scalability, in terms of memory and computational requirements. One of them also exploits randomisation, over data blocks at each iteration, offering further flexibility. We present simulation results showing the feasibility of the proposed methods as well as their advantages compared to state-of-the-art algorithmic solvers. Our Matlab code is available online on GitHub.