The w-effect in interferometric imaging: from a fast sparse measurement operator to super-resolution

TL;DR

This paper explores how the w-effect in interferometric imaging, especially with large w-modulations, can be modeled efficiently using sparse operators, enhancing super-resolution capabilities and reducing computational costs in radio astronomy imaging.

Contribution

It introduces sparse, fast models of the measurement operator for interferometric imaging that leverage the spread spectrum effect for improved super-resolution, with analysis of their impact on reconstruction quality.

Findings

Sparse models significantly reduce memory and computation.

Large w-modulations enhance super-resolution potential.

Sparse operators maintain acceptable image reconstruction quality.

Abstract

Modern radio telescopes, such as the Square Kilometre Array (SKA), will probe the radio sky over large fields-of-view, which results in large w-modulations of the sky image. This effect complicates the relationship between the measured visibilities and the image under scrutiny. In algorithmic terms, it gives rise to massive memory and computational time requirements. Yet, it can be a blessing in terms of reconstruction quality of the sky image. In recent years, several works have shown that large w-modulations promote the spread spectrum effect. Within the compressive sensing framework, this effect increases the incoherence between the sensing basis and the sparsity basis of the signal to be recovered, leading to better estimation of the sky image. In this article, we revisit the w-projection approach using convex optimisation in realistic settings, where the measurement operator couples the w-terms in Fourier and the de-gridding kernels. We provide sparse, thus fast, models of the Fourier part of the measurement operator through adaptive sparsification procedures. Consequently, memory requirements and computational cost are significantly alleviated, at the expense of introducing errors on the radio-interferometric data model. We present a first investigation of the impact of the sparse variants of the measurement operator on the image reconstruction quality. We finally analyse the interesting super-resolution potential associated with the spread spectrum effect of the w-modulation, and showcase it through simulations. Our C++ code is available online on GitHub.

Equations

$$y(\mathbf{u},w)=\int n(\mathbf{l}{)}^{-1}a(\mathbf{l})e^{-2i\pi w(n(\mathbf{l})-1)}\tilde{x}(\mathbf{l})e^{-2i\pi \mathbf{u}\cdot \mathbf{l}}{\mathrm{d}}^2\mathbf{l},$$

$$\mathbf{y}=\mathsf{G}\mathsf{F}\mathsf{Z}\mathsf{D}\mathbf{x},$$

$$y_{\ell }={\tilde{\mathbf{g}}}_{\ell }^{\top }\mathsf{F}\mathsf{Z}\mathsf{D}\mathbf{x}\text{, with }{\tilde{\mathbf{g}}}_{\ell }={\hat{\mathbf{c}}}_{\ell }\ast {\mathbf{g}}_{\ell },$$

$$\mathbf{y}=\tilde{\mathsf{G}}\mathsf{F}\mathsf{Z}\mathsf{D}\mathbf{x}+\mathbf{n}=\mathbf{\Phi }\mathbf{x}+\mathbf{n},$$

$$\displaystyle\min_{\boldsymbol{x}}\parallel{\boldsymbol{\Psi}}^{\dagger}{\boldsymbol{x}}\parallel_{1}\text{ s.t. }\left\{\begin{array}[]{l}\parallel{\boldsymbol{y}}-{\boldsymbol{\Phi}}{\boldsymbol{x}}\parallel_{2}\leq\epsilon,\\ {\boldsymbol{x}\geq\boldsymbol{\mathsf{0}}},\end{array}\right.$$

$$\underset{\mathbf{x}}{\min }\parallel {\mathbf{\Psi }}^{†}\mathbf{x}{\parallel }_1+{\iota }_{{\mathbb{R}}_{+}^{\mathbb{N}}}(\mathbf{x})+{\iota }_{\mathcal{B}}(\mathbf{\Phi }\mathbf{x}),$$

$${\iota }_{\mathcal{C}}(\mathbf{z})\stackrel{\Delta }{=}\{\begin{array}{rl}0& \mathbf{z}\in \mathcal{C}\\ +\infty & \mathbf{z}\notin \mathcal{C}\end{array}$$

$$\underset{\mathbf{x}}{\min }{g}_1(\mathbf{x})+g_2(\mathbf{x})\cdots +g_n(\mathbf{x}),$$

$$\Big{(}\operatorname{\boldsymbol{\mathcal{H}}}_{\tau_{\ell}}(\boldsymbol{r}_{\ell})\Big{)}_{i}\overset{\Delta}{=}\left\{\begin{array}[]{cl}\displaystyle r_{i,\ell}&\qquad|r_{i,\ell}|>\tau_{\ell}\\ \displaystyle 0&\qquad|r_{i,\ell}|\leq\tau_{\ell}\\ \end{array}\right.\quad\forall i,$$

$$\mathbf{y}=\tilde{\mathsf{G}}\mathsf{F}\mathsf{Z}\mathsf{D}\mathsf{U}\mathbf{x}+\mathbf{n}=\mathbf{\Phi }\mathbf{x}+\mathbf{n},$$

$$\mathrm{SNR}=20{\log }_{10}\frac{||\mathbf{x}|{|}_2}{||\mathbf{x}-\bar{\mathbf{x}}|{|}_2},$$

$$\mathrm{MR}=\frac{1}{N}\sum _{i=1}^N\min (\frac{x_i}{{\bar{x}}_i},\frac{{\bar{x}}_i}{x_i}),$$

$$\mathrm{SNR}=20{\log }_{10}\frac{||\mathbf{z}|{|}_2}{||\mathbf{z}-\bar{\mathbf{z}}|{|}_2},$$

Full text

††pubyear: 2017††pagerange: The $w$-effect in interferometric imaging: from a fast sparse measurement operator to super-resolution–References

1 Introduction

radio interferometric imaging has been subject of numerous developments in the recent years, with the rise of the new-generation radio telescopes, namely the Square Kilometre Array (SKA)111http://www.skatelescope.org/ and its pathfinders such as the Low-Frequency Array (LoFAR)222http://www.lofar.org/, the upgraded Karl G. Jansky Very Large Array (VLA)333https://science.nrao.edu/facilities/vla. These instruments are characterised with extremely high resolution and sensitivity along with their capabilities of mapping large regions of the radio sky. To meet such specifications and hence, deliver the expected science goals, an urgent need to revolutionize radio interferometric imaging arises. Newly proposed imaging techniques in the last decade (Wiaux et al., 2009a; Wenger et al., 2010; McEwen & Wiaux, 2011; Li et al., 2011; Carrillo et al., 2012, 2014; Dabbech et al., 2012; Dabbech et al., 2015; Garsden et al., 2015; Onose et al., 2016; Pratley et al., 2016; Onose et al., 2017), have been built almost exclusively based on sparse regularisations, and they have shown to outperform the conventional methods, namely CLEAN and its variants (Högbom, 1974; Wakker & Schwarz, 1988; Cornwell et al., 2008; Pratley & Johnston-Hollitt, 2016). These approaches solve the two-dimensional (2D) radio interferometric imaging problem, in the absence of propagation and receiver errors and within small fields of view. Under these simplified conditions, the problem consists of the recovery of the radio image from its partial and noisy Fourier samples. The imaging problem is ill-posed due to the incompleteness of the Fourier sampling, resulting from the finite number of antennas, and is challenging per se.

In the context of wide-field radio interferometric imaging, the non-coplanarity of the radio interferometric array, its long baselines and its probed large field-of-view give rise to the so-called $w$-effect, which originates from the third dimensionality of the Fourier space where visibilities are measured. This results in a modulation of the imaged radio sky with a chirp-like phase term parametrised by the $w$ component of the baseline in the direction of sight, also referred to as the $w$-term. Consequently, the radio interferometric imaging problem becomes three-dimensional, requiring a massive increase of memory requirements and computational cost. Yet, ignoring this effect in the imaging step can hamper significantly the quality of the reconstructed image. In the literature, several approaches have been proposed to correct for the $w$-effect, namely $w$-projection (Cornwell et al., 2008) and $w$-stacking (Offringa et al., 2014). These techniques have demonstrated a reasonable compromise between memory requirements and computational cost on the one hand, and the quality of the recovered image on the other hand with respect to the traditional approaches such as faceting (Cornwell & Perley, 1992). The latter exploits the fact that the $w$ component is smaller close to the phase center. The imaged radio sky is thus decomposed into facets, where on each facet, the $w$-effect is considered flat, allowing for the use of a single point spread function (PSF) in the deconvolution step. The $w$-projection approach relies on the use of $w$-dependent convolution kernels in the Fourier domain. The data are obtained by projecting the sky Fourier transform from the plane $(u,v,w=0)$ to the correct $(u,v,w)$ plane. This implies that the two-dimensional Fourier transform can be used, hence preserving the advantages of FFTs. Whereas the $w$-stacking approach operates rather on the image domain; visibilities are gridded on discrete $w$-layers and the correction is performed on the image domain.

Recently, $w$-projection within the context of convex optimization has been introduced (Wiaux et al., 2009b; Wolz et al., 2013). In Wiaux et al. (2009b), the authors presented a theoretical study of the spread spectrum effect promoted by the $w$-terms within the compressive sensing theory. The latter states that the unknown signal can be accurately recovered from a few number of measurements provided that the sensing basis and the signal’s sparsity basis are incoherent. In this context, the authors have presented a seminal study of how a constant $w$-modulation increases the incoherence between the sensing basis and the signal sparsity basis. As a result, the quality of the imaged radio sky is significantly enhanced. Wolz et al. (2013) demonstrated the spread spectrum effect for varying $w$ values, considering that the data arise from sparse versions of the $w$-terms in the Fourier domain, in order to lower memory requirements. In this article, we revisit the $w$-projection for realistic settings within the same framework. We study the image reconstruction quality when the original visibilities, arising from exact non-sparse $w$-kernels, are modelled with sparse approximations of the $w$-kernels. These models give rise to errors on the measurement operator. We provide a first study of their limits and effects on the quality of the image recovery as well as their efficiency in terms of memory requirements and accelerated computational time. Though, the general tendency is to reduce or neglect the $w$-terms due to their resulting massive computational and memory demands, we provide a proof of concept that the $w$-modulation has a large potential to promote super-resolution, recovering images beyond the resolution of the radio interferometer limited by its largest projected baseline.

The remainder of the paper is organised as follows. Section 2 describes the radio interferometric imaging problem in the presence of the $w$-term. In Section 3, we recall the minimisation problem in radio interferometric imaging that we solve using PURIFY, a software for radio interferometric imaging (Carrillo et al., 2014; Pratley et al., 2016), into which we have incorporated $w$-projection. In Section 4, we describe two sparsification strategies for approximating the measurement operator and also study the effect of sparsification on the image reconstruction quality. In Section 5, we analyse and showcase the potential of the $w$-term to promote super-resolution. Finally, we state our conclusions in Section 6.

2 Radio interferometric imaging problem

radio interferometers probe the electric field $\vec{E}$ coming from the sky, each acquired measurement, called visibility, corresponds to the cross-correlation of the sensed signals by an antenna pair. Assuming a monochromatic radiation, the vectorial distance separating two antennas $(p,q)$ and denoted by $\vec{\bf b}_{pq}$ is described by its coordinates in units of the wavelength $(u,v,{w})$, with $w$ is the coordinate in the direction of sight and ${\boldsymbol{u}}=(u,v)$ are the relative coordinates on its perpendicular plane. The radio interferometric measurement $y({\boldsymbol{u}})$ is related to the sky surface brightness $\tilde{x}$ through:

[TABLE]

where ${\boldsymbol{l}}=(l,m)$ are the coordinates of a source in a tangent plane to the celestial sphere and $n({\boldsymbol{l}})=\sqrt{1-|{\boldsymbol{l}}|^{2}}$ its coordinate on the line of sight, with $|{\boldsymbol{l}}|=\sqrt{l^{2}+m^{2}}$. We explicitly introduce the $w$-modulation term as a function $c(w,{\boldsymbol{l}})=e^{-2i\pi{w}(n({\boldsymbol{l}})-1)}$ into (1). The function $a({\boldsymbol{l}})$ is the illumination function of the antenna and depends on the source position. Other direction dependent effects (DDEs), describing unknown instrumental errors and ionospheric perturbations, can be incorporated in the same manner as the illumination function $a({\boldsymbol{l}})$, even if not done so explicitly here. Note that, the term $c(w,{\boldsymbol{\mathsf{;}}})$ is also a DDE as it depends on the source position through $n({\boldsymbol{l}})$, yet unlike the mentioned DDEs it is known, thus easier to model accurately although at the expense of large memory requirements. From now on, we call the sky intensity $x({\boldsymbol{l}})\equiv a({\boldsymbol{l}})\tilde{x}({\boldsymbol{l}})$; that is the modulated sky surface brightness with the illumination function. Here, we assume that the illumination function is known and invariant. However, this is not always the case in practice; $a({\boldsymbol{l}})$ underlies time-dependent variations, and consequently it can be described as a DDE requiring calibration.

When observing a narrow field of view $|{\boldsymbol{l}}|^{2}\ll 1$ (see Wiaux et al., 2009b) and also when the array is co-planar ${w}=0$, the ${w}$-term corresponds to a flat function $c(w,{\boldsymbol{l}})=1,\forall{\boldsymbol{l}}$. In this case, each visibility reduces to a Fourier component of the sky intensity image. These approximations, however, break down in the case of wide-field imaging; modern radio interferometers will probe large FoVs and their antennas span hundreds and even thousands of kilometres. Futhermore, the radio sky is rather spherical (McEwen & Scaife, 2008; McEwen & Wiaux, 2011). Nevertheless, to remain in the context of a 2D imaging, in the present study, we consider a relatively small field-of-view (FoV) with non-negligeable $w$-terms. Under this assumption, we adopt a first order approximation of the ${w}$-term that is a linear chirp $c(w,{\boldsymbol{l}})=e^{i\pi{w}{\boldsymbol{l}^{2}}}$. The $w$-term is thus considered as a norm-preserving phase modulation of the sky intensity map $x({\boldsymbol{l}})$ with a linear chirp of rate ${w}$. This modulation has been well studied in Wiaux et al. (2009b); Wolz et al. (2013) where the authors have demonstrated that the chirp modulation induces the spread spectrum effect, that is in favour of a better image reconstruction quality of the sparsity-based approaches. Note that the spread spectrum effect is not restricted to the $w$-modulation, in fact, all DDEs inherently induce similar effects due to their convolution nature in the Fourier domain.

2.1 The de-gridding operator

In order to promote the use of the Fast Fourier Transform (FFT), radio measurements must lie on a regular grid (usually with a size of a power of 2). Therefore, they are modelled from uniformly sampled Fourier components via linear operations, where the continuous (i.e., off the grid) measurements are interpolated using convolution kernels in the Fourier domain. In absence of the ${w}$-term and the noise, the radio interferometric measurement equation in (1) can be recast in a discrete setting as follows

[TABLE]

where the data ${\boldsymbol{y}}\in\mathbb{C}^{M}\equiv\{y_{\ell}\equiv y({\boldsymbol{u}}_{\ell},w_{\ell})\}_{1\leq\ell\leq M}$ are complex measurements of size $M$. The signal to be recovered ${\boldsymbol{x}}\in\mathbb{R}_{+}^{N}$ is positive and of size $N$. The operator $\boldsymbol{\mathsf{F}}$ is the Fourier transform operator. The matrix ${{\boldsymbol{\mathsf{G}}}}{\in\mathbb{C}}^{(M,\alpha^{2}N)}$ is sparse. Its rows ${{\boldsymbol{g}}}_{\ell}\in\mathbb{C}^{\alpha^{2}N}$ are convolution kernels with $P\ll N$ non-zero elements, called de-gridding kernels, each centred around its corresponding continuous spacial frequency point ${\boldsymbol{u}}_{\ell}=(u_{\ell},v_{\ell})$. The diagonal matrix ${\boldsymbol{\mathsf{D}}}\in{\mathbb{R}}^{(\alpha^{2}N,\alpha^{2}N)}$ is correcting for the convolution in the Fourier domain; when applied to the sky, it corresponds to point-wise division of $\boldsymbol{x}$ with the inverse Fourier transform of the de-gridding kernel. The operator $\boldsymbol{\mathsf{Z}}$ is an oversampling operator by a factor $\alpha^{2}$ (that is a factor $\alpha$ on each dimension often set to 2), allowing for a fine discrete uniform sampling in the Fourier domain via zero-padding of the image to be recovered.

The de-gridding operator, not only allows for modelling the continuous measurements, but also corrects for aliasing. In fact, it attenuates the sky-image brightness at the boundaries of the region of interest. The choice and the support (i.e., number of non-zero elements) of the convolution kernels is crucial. Ideally one would use the Sinc function for an accurate interpolation. Adopting this function translates into a multiplication of the imaged sky with a box function in the image domain, thus cancelling aliasing artefacts from sources outside the imaged region of the radio sky. However, this function is not band limited i.e., it has an infinite support in Fourier, thus leading to heavy computational cost comparable to the use of the Direct Fourier Transform. Therefore, band limited convolution kernels with a compact support in the Fourier domain are usually preferred. In radio interferometry, prolate-spheroidal wave functions (PSWF) are widely adopted (Schwab, 1984; Sault et al., 1995) with a support size of $6\times 6$ pixels. Kaiser-Bessel functions are also used (Offringa et al., 2014) often with a support size $8\times 8$ pixels. Unlike PSWFs, the Kaiser-Bessel kernels present the advantage of having an analytic expression that can be evaluated easily and accurately. Furthermore, they can be adopted with lower support size. In general, both functions lead to comparable performance (see Greisen, 1998; Pratley et al., 2016, for further discussion).

2.2 The measurement operator

In this section, we introduce the $w$-term into the discrete setting of (2). In the image domain, the ${w}$-modulation can be expressed as a multiplication of the sky image with the ${w}$-term. The latter can be approximated by a linear chirp ${\boldsymbol{c}_{w}}=\{e^{i\pi{{w}}\|{\boldsymbol{l}}\|^{2}}\},\forall~{}\boldsymbol{l}$ pixel positions. Equivalently, in the Fourier plane, the modulation corresponds to the convolution of the Fourier transform of the linear chirp with that of the signal. Since the $w$-rate is baseline-dependent, each sensed measurement $y_{\ell}$ is the coefficient at the spacial frequency ${\boldsymbol{u}}_{\ell}$ of the signal resulting from the convolution of the Fourier transform of the linear chirp $\hat{\boldsymbol{c}}_{\ell}$, characterised by the rate $w_{\ell}$, with that of the sky image $\hat{\boldsymbol{x}}$. This means that the resulting sensed signal affected by the $w$-modulation is not the 2D conventional convolution but rather the outcome of a 2D convolution with a visibility-dependent kernel.

For a narrow FoV denoted by $L$, the band limit of a linear chirp with a rate ${w}_{\ell}$, is deduced from its instantaneous frequency and is approximated as $B_{{w}_{\ell}}={|{w}_{\ell}|L}/{2}$ (Wiaux et al., 2009b). The largest bandwidth of the modulated signal $B^{\prime}_{\rm max}$ is given by $B^{\prime}_{\rm max}=\displaystyle\max_{{1\leq\ell\leq M}}{B_{w_{\ell}}}+{\parallel{\boldsymbol{u}_{\ell}}\parallel_{2}}$. Note that, in absence of $w$-modulation, the bandwidth of the signal to be recovered is $B=\displaystyle\max_{1\leq\ell\leq M}\|\boldsymbol{u}_{\ell}\|_{2}$, which constitutes the band limit of the radio interferometer. The true radio sky is characterized with an infinite bandwidth. Therefore, a correct modelling of the visibilities requires accommodating the bandwidth of the unknown signal, at least, up to the largest bandwidth $B^{\prime}_{\rm max}$, that is accessible by the highest $w$-modulation. This translates in choosing a pixel-size $\Delta l^{\prime}=1/2B^{\prime}_{\rm max}$ in the image domain, that is smaller than the radio interferometer resolution (i.e. $\Delta l=1/2B$). Consequently, the imaging step consists in estimating the unknown signal at a resolution above the instrument band limit. Yet, in practice a convolution of the deconvolved image of the sky with a low-pass filter, that is smoothing kernel corresponding to the instrument resolution, is often applied. Hence, all the spatial Fourier modes above the instrument’s band limit are set to zero. When accounting for the $w$-term, each measurement $y_{\ell}$ reads:

[TABLE]

where $\ast$ denotes linear convolution. The estimated image $\boldsymbol{x}$ is of size $N^{\prime}$, with a bandwidth $B^{\prime}_{\rm max}$. The vector ${\hat{\boldsymbol{c}}}_{\ell}\in\mathbb{C}^{{N}^{\prime\prime}}$ denotes the Fourier transform of the chirp modulation of size ${N}^{\prime\prime}=\alpha^{2}N^{\prime}$ and support $K_{{w}_{\ell}}<N^{\prime\prime}$. $\tilde{\boldsymbol{g}}_{\ell}\in\mathbb{C}^{{N}^{\prime\prime}}$ are convolution kernels coupling the de-gridding interpolation kernels and the ${w}_{\ell}$-terms, with a support $N_{\ell}=P+K_{{w}_{\ell}}$. The radio measurements are hence modelled as:

[TABLE]

where $\tilde{\boldsymbol{\mathsf{G}}}$ is the $w$-projection operator, whose rows correspond to ${\tilde{\boldsymbol{g}}}_{1\leq\ell\leq M}\in\mathbb{C}^{N^{\prime\prime}}$, and ${\boldsymbol{n}}\in\mathbb{C}^{M}$ is additive white Gaussian noise, with known statistics, modelling the instrumental noise. The operator ${\boldsymbol{\Phi}}\equiv{\tilde{\boldsymbol{\mathsf{G}}}}{{\boldsymbol{\mathsf{F}}}{\boldsymbol{\mathsf{Z}}}{\boldsymbol{\mathsf{D}}}}$ is the measurement operator, incorporating the ${w}$-modulation. The problem of recovering the signal $\boldsymbol{x}$ from the measurements $\boldsymbol{y}$ is ill-posed due to the incomplete sampling and the noise i.e, it does not present a unique solution. Therefore, additional prior information on the signal to be recovered $\boldsymbol{x}$ is usually imposed to better constrain the imaging problem and consequently obtaining a good approximate of the unknown image.

3 Sparse image reconstruction

Numerous works have demonstrated the applicability of compressive sensing theory to solve the radio interferometric imaging problem. Within this framework, the unknown signal is assumed to be sparse or compressible in a data representation space ${\boldsymbol{\Psi}}$ (e.g. wavelets, curvelets, discrete cosine transform). This is generally the case for natural signals. The theory of compressive sensing states that the exact recovery of the unknown signal is possible from a small number of measurements $M$ that is below the Shannon-Nyquist sampling rate, provided that the sensing operator $\boldsymbol{\Phi}$ and the signal’s sparsity basis ${\boldsymbol{\Psi}}$ are incoherent (Candès, 2006; Donoho, 2006); the mutual coherence defined as the largest cross-correlation between the columns of ${\boldsymbol{\Psi}}$ and ${\boldsymbol{\Phi}}$.

In radio interferometry, the Fourier basis constitutes the sensing basis. Furthermore, assuming sparsity-by-synthesis of the signal $\boldsymbol{x}$ in a dictionary $\boldsymbol{\Psi}$, that is $\boldsymbol{x}=\boldsymbol{\Psi}\boldsymbol{\alpha}$, where the synthesis coefficients vector $\boldsymbol{\alpha}$ is sparse, the mutual coherence of the sparsity and the sensing basis is defined as the largest modulus of the Fourier coefficients of the basis functions constituting the sparsity basis $\boldsymbol{\Psi}$. For instance, considering a signal of size $N$, sparse in its domain, its sparsity basis $\boldsymbol{\Psi}$ is the Dirac basis and the mutual coherence of the Fourier basis and the Dirac basis is given by $1/\sqrt{N}$. For large $N$, the mutual coherence goes to zero. Therefore, the Dirac basis is totally incoherent with the Fourier basis. For more complex signals, sparsity is assumed in more sophisticated dictionaries such as wavelets. The latter can be sufficiently incoherent with the Fourier basis. Note that, the choice of the adequate sparsity basis is important and is application-dependent. In radio astronomy, several dictionaries have been adopted and have shown to be relevant, namely wavelets such as the isotropic undecimated wavelet transform adopted in Li et al. (2011); Dabbech et al. (2012); Dabbech et al. (2015); Garsden et al. (2015), the Daubechies wavelets concatenated with the Dirac basis in Carrillo et al. (2012); Wolz et al. (2013); Carrillo et al. (2013); Carrillo et al. (2014); Onose et al. (2016); Pratley et al. (2016); Onose et al. (2017).

The $w$-term, as a norm-preserving modulation, spreads the energy of the Fourier coefficients of the sparsity dictionary’s basis functions over the neighbouring Fourier modes, reducing the amplitudes of these Fourier coefficients. Consequently, the mutual coherence of the sparsity dictionary and the sensing basis is significantly decreased. Therefore, intrinsically, the $w$-term reinforces the embedding of the radio interferometric imaging problem within the compressive sensing theory, a detailed theoretical study on the matter is provided in Wiaux et al. (2009b).

3.1 Imaging with convex optimisation

To solve the inverse problem given in (4), we adopt sparsity-by-analysis promoting prior as a reguliser. This imposes that the projection of the unknown signal in a sparsity basis $\boldsymbol{\Psi}$ is sparse, that is $\boldsymbol{\alpha}=\boldsymbol{\Psi}^{\dagger}\boldsymbol{x}$ and the analysis vector $\boldsymbol{\alpha}$ being sparse. The adopted sparsity basis is that proposed in Carrillo et al. (2012), $\boldsymbol{\Psi}$ is over-complete and is chosen to be the concatenation of nine basis; the Dirac basis and the eight Daubechies wavelet basis. The minimisation problem reads:

[TABLE]

where the bound $\epsilon$ is the $\ell_{2}$ norm of the noise $\boldsymbol{n}$. To impose sparsity on $\boldsymbol{\Psi}^{\dagger}\boldsymbol{x}$, intuitively one should minimise the $\ell_{0}$ norm of $\boldsymbol{\Psi}^{\dagger}\boldsymbol{x}$, however this leads to NP-hard problem. Therefore, the $\ell_{0}$ norm is often relaxed by adopting the $\ell_{1}$ norm, that is sparsity-promoting and convex. Re-weighted schemes of the $\ell_{1}$ norm have been proposed in the literature to better approximate the $\ell_{0}$ norm (Candès et al., 2008). Furthermore, they have shown considerable improvements of the image recovery quality (Carrillo et al., 2012, 2014).

The constrained minimisation problem (5) can be efficiently solved using convex optimisation, furthermore, it can be written as follows

[TABLE]

where, the constraints defined in (5) are imposed using the indicator function $\iota_{\mathcal{C}}$, defined on a non-empty set ${\mathcal{C}}$ as

[TABLE]

Minimising the function $\iota_{\mathcal{\mathbb{R_{+}^{N}}}}$ imposes the positivity and reality constraints of the signal to be recovered $\boldsymbol{x}$. The function $\iota_{\mathcal{B}}$ applied to $\boldsymbol{\Phi}\boldsymbol{x}$ is the data fidelity term. When minimised, it constrains the residual to be in the $\ell_{2}$ ball $\mathcal{B}$. The latter is a convex set, characterised by the noise level $\epsilon$ such that $\mathcal{B}=\{\boldsymbol{z}\in\mathbb{C}^{M}:\|\boldsymbol{z}-\boldsymbol{y}\|_{2}\leq\epsilon\}$.

The problem (6) fits well within proximal splitting methods. These are solving for optimisation problems of the form

[TABLE]

where the functions $g_{i\in\{1\dots n\}}$ are proper, convex and lower semi-continuous. Each function is minimised individually. Solutions to differentiable functions, such as the $\ell_{2}$ norm, are obtained using their gradient. While, minimising non-smooth functions involve the use of their proximity operators, in particular the $\ell_{1}$ norm is minimised through its proximal operator, that is the soft-thresholding. We direct the reader to Combettes & Pesquet (2009) for a good review. These methods have gained wide interest in the recent years thanks to the splitting of the functions, hence their scalability to large-scale problems. In particular, several proximal splitting methods have been proposed to solve (6) for radio interferometry. Carrillo et al. (2012) used the Douglas-Rachford splitting method (Combettes & Pesquet, 2007), Carrillo et al. (2014) adopted the simultaneous-direction method of multipliers (SDMM) (Setzer et al., 2010). More recently, Onose et al. (2016) have proposed two highly parallelisable and distributed algorithmic structures based on the alternating direction method of multipliers (ADMM) (Boyd et al., 2011) and the Primal-Dual (PD) algorithm (Komodakis & Pesquet, 2015), allowing for full splitting of the data and the linear operator into blocks.

3.2 Radio interferometric imaging in PURIFY

In our study, we use the software PURIFY444http://basp-group.github.io/purify/ (Carrillo et al., 2014; Pratley et al., 2016), that is a package written in C++ for radio interferometric imaging and using $\ell_{1}$ regularisation-based solvers, in particular, the SARA algorithm originally proposed in Carrillo et al. (2012). PURIFY supports standard de-gridding kernels (Pratley et al., 2016) in particular Kaiser-Bessel convolution kernels, with optimized parameters suggested in Fessler & Sutton (2003) and the prolate spheroidal wave kernels as described in Schwab (1984). We have implemented the $w$-terms as described in (3) leading to the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$ detailed in Section 2. Variants of the measurement operator are also included where sparse models of the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$ through adaptive sparsification strategies of its convolution kernels are adopted, these are detailed in Section 4.

PURIFY calls the Sparse OPTimisation (SOPT) software package555http://basp-group.github.io/sopt/, a collection of solvers of the convex problem (6), namely the simultaneous-direction method of multipliers (SDMM) algorithm proposed by Carrillo et al. (2014) and the proximal alternating direction method of multipliers proposed by Onose et al. (2016). The PD-based splitting algorithm proposed in the latter article will be included in a future release of SOPT. While SDMM is known to be computationally heavy since it involves matrix inversion, both algorithmic structures proposed in Onose et al. (2016) are parallelisable and allow for an efficient distributed implementation, through the splitting of the data into blocks666In particular, PD-based splitting algorithm has shown to be very flexible allowing for full data and operators splitting, resulting in lower computational burden and memory requirements. It is also prone to further increased scalability by using randomized updates; where only randomly selected data blocks are processed within each iteration.. In the present study, we have adopted the ADMM algorithm as a solver. Note that, the block splitting feature will be included in a future release of PURIFY.

4 Sparse representations of the measurement operator

The $w$-projection approach allows for an accurate correction of the $w$-term, accounting for the exact $w$-modulation through its discrete Fourier transform. The resulting $w$-terms, however, are usually non-sparse resulting in a high-dimensional, non-sparse measurement operator. In fact, in the presence of high $w$-rates their corresponding chirp kernels $\hat{\boldsymbol{c}}$ do not decay rapidly to zero. Moreover, even in the presence of small $w$-rates, their corresponding $\hat{\boldsymbol{c}}$ can present large amounts of very small coefficients yet non-zero valued, and so they are rather compressible but not strictly sparse, the same applies for the resulting $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$. This is computationally demanding as it involves a heavy measurement operator ${\boldsymbol{\Phi}}$. In practice, Cornwell et al. (2008) suggest to consider truncated $w$-terms such that coefficients below $10\%$ of the maximum coefficient in terms of modulus are set to zero. Moreover, the $w$-rates are further sampled on a discrete small set of values. These considerations result in faster operators, though they might introduce large errors on the model of the measurement operator, hence hampering the reconstruction quality and the dynamic range in particular.

We shall define the Fourier-transformed Chirp operator $\hat{\boldsymbol{\mathsf{C}}}$ whose rows ${\hat{\boldsymbol{c}}}_{1\leq\ell\leq M}\in\mathbb{C}^{N^{\prime\prime}}$ are the Fourier transforms of the linear chirp modulations with the rates $w_{{1\leq\ell\leq M}}$, centred at the zero frequency. In the context of discrete measurements on a uniformly sampled grid, no de-gridding kernels are involved and the $w$-projection operator is given by the Fourier transformed $w$-terms ${\hat{\boldsymbol{c}}}_{\ell}$, each centred at its corresponding $\boldsymbol{u}_{\ell}$ point. In this case, Wolz et al. (2013) show that the $w$-projection operator, when expressed as a sparse matrix, results in low effective computational time and memory requirements for imaging, provided that it is remains accurate.

In this study, we present two approaches for sparse representations of the measurement operator. The first strategy is based on the sparsification of the Fourier-transformed Chirp operator $\hat{\boldsymbol{\mathsf{C}}}$ before the convolution with the de-gridding operator ${\boldsymbol{\mathsf{G}}}$ and resulting in a sparse $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$. The second approach consists in the direct sparsification of the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$. In the first approach, thanks to the reduced number of the non-zero elements of $\hat{\boldsymbol{\mathsf{C}}}$, the process of the row-wise convolution of the de-gridding operator ${\boldsymbol{\mathsf{G}}}$ and $\hat{\boldsymbol{\mathsf{C}}}$ is accelerated since the number of non-zero elements in each $w$-term is significantly reduced. This is not the case in the direct sparsification of the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$, since the convolution is performed using the non sparse $w$-terms, though in general this approach yields more sparsity of the measurement operator as it is shown in this section.

To determine a sparse representation of the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$, we adopt the sparsification technique proposed in Wolz et al. (2013). The idea is to apply an adaptive hard-thresholding on each row ${\boldsymbol{r}}_{\ell}$ of an operator $\boldsymbol{\mathsf{R}}\in\mathbb{C}^{M\times N}$, that is a component-wise operation. The hard-thresholding operator denoted by $\operatorname{\boldsymbol{\mathcal{H}}}$ is defined as follows

[TABLE]

with the constraint of loosing a fixed energy fraction $\gamma$ across all the rows. The energy of a vector $\boldsymbol{r}_{\ell}\in\mathbb{C}^{M}$ is defined as the sum of the squared modulus of its elements $E^{\rm{total}}_{\ell}\equiv\sum_{i=1}^{N}|r_{\ell,i}|^{2}$. For a preserved energy $E^{\rm{sparse}}_{\ell}=(1-\gamma)\times E^{\rm{total}}_{\ell}$, where $0\leq\gamma\leq 1$, the significant non-zero elements of $\boldsymbol{r}_{\ell}$ are determined by computing a threshold $\tau_{\ell}$ using a bisecting method such that elements contributing to the energy $E^{\rm{lost}}_{\ell}=E^{\rm{total}}_{\ell}-E^{\rm{sparse}}_{\ell}$ are set to zero. This approach allows to adaptively reduce the number of significant elements for any given operator while not limiting the process to a specific support size. A sparse approximation of the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$ is obtained by applying this sparsification technique on the rows ${\hat{\boldsymbol{c}}}_{{1\leq\ell\leq M}}$ of $\hat{\boldsymbol{\mathsf{C}}}$ or the rows $\tilde{\boldsymbol{g}}_{1\leq\ell\leq M}$ of $\tilde{\boldsymbol{\mathsf{G}}}$.

The present study consists in studying the effect of adopting sparse approximations of the $w$-projection operator in the image recovery step. This results in model errors of the radio interferometric measurements. We therefore provide limit values of the sparsification levels considered in the sparse approximations of the $w$-projection operator, for additive noise $\boldsymbol{n}$ characterised by two input signal-to-noise ratios. In these limit cases, degradation of the image reconstruction quality is observed; the errors of the adopted backward model are no longer buried in the noise. We shall recall that in Wolz et al. (2013), the authors have shown that sparse Chirp operators do enhance the image reconstruction quality. Exact sparse $w$-projection operators have been employed for imaging; in both the forward model, adopted for radio interferometric data simulation following (10), and the backward model, considered in the image recovery using the approach described in (5).

4.1 Simulations settings

We study the explained strategies for a sparse measurement operator through realistic simulations of radio interferometric data. A realistic $uv$-coverage is simulated using the Meqtrees software (Noordam & Smirnov, 2010) from the antenna configuration of the Australian Square Kilometre Array Pathfinder (ASKAP)777http://www.atnf.csiro.au/projects/askap/index.html (see Fig. 1, left panel). The coverage is obtained for a total observation time of 1 hour and time spacing $\delta{\rm t}=5$ minutes, pointing at declination $-10d0m0s$ and right-ascension $0d0m0s$. The total number of $\boldsymbol{u}$ points is $M=7560$. The $w$ components generated in this setting are extremely small leading to chirp kernels that can be easily approximated by a Dirac for nearly $90\%$ of the $uv$-coverage. Therefore, in our tests, we consider artificially large $w$ values, which we generate as zero-mean uniform distribution $w$ values - hence not correlated to the $uv$ points distribution -, such that the maximum $w$-rate $|w_{\rm max}|=w_{f}\times w^{*}$, where $w^{*}=2B/{L^{2}}$ is the $w$-rate resulting in a $w$-modulation having the bandwidth $B$ of the signal. Hence, the parameter $w_{f}$ translates the fraction of the maximum $w$-modulation with respect to $w^{*}$. In Fig. 1, we showcase an example of uniformly distributed $w$ values for $w_{f}=0.1$ versus the naturally generated $w$ values, where the maximum $w$-rate corresponds to $w_{f}=8.5\times 10^{-4}$.

The considered ground-truth sky image is an image of a galaxy cluster ($\rm{GC}$), obtained using the FARADAY tool (Murgia et al., 2004) (see Fig. 2, right panel). The $\rm{GC}$ image is a collection of very bright compact sources corresponding to galaxies and a central diffuse faint emission that corresponds to a radio halo. Beside its diversity as an image, in terms of structure, it is also characterized with a high dynamic range of order $10^{6}$. We consider a low resolution version of the GC image, that is of size $N=128\times 128$ pixels with its highest pixel value scaled to 1 and its bandwidth is equal to the radio interferometric band limit, that is $B=\displaystyle\max_{1\leq\ell\leq M}\parallel{\boldsymbol{u}_{\ell}}\parallel_{2}$. The observed field of view is narrow with $L=0.21$ degrees and is determined by the simulated baseline coverage and the size of the ground-truth image. We generate radio interferometric visibilities at the frequency 1GHz, in MATLAB. The measurements are corrupted with additive white Gaussian noise, with input signal-to-noise ratio (iSNR) of 30dB and 40dB. For the de-gridding kernels, we adopt the Kaiser-Bessel interpolation kernels and fix their support size to $4\times 4$ pixels. Note that, to accommodate the full $w$-modulation kernels up to the bandwidth $B^{\prime}_{\rm max}$ as described in subsection 2.2, we apply an up-sampling operator to the sky image, which consists in performing a zero-padding in the Fourier domain. Simulated observations are thus given by:

[TABLE]

where $\boldsymbol{\mathsf{U}}$ is the up-sampling operator such that the up-sampled image $\boldsymbol{\mathsf{U}}\boldsymbol{x}$ is of size $N^{\prime}$ (determined with respect to the highest $w$-modulation considered, see subsection 2.2). The measurement operator is $\boldsymbol{\Phi}=\tilde{\boldsymbol{\mathsf{G}}}{\boldsymbol{\mathsf{F}}}{\boldsymbol{\mathsf{Z}}}{\boldsymbol{\mathsf{D}}}{\boldsymbol{\mathsf{U}}}$, and its adjoint operator is ${\boldsymbol{\Phi}}^{\dagger}=\tilde{\boldsymbol{\mathsf{U}}}{\boldsymbol{\mathsf{D}}}\tilde{\boldsymbol{\mathsf{Z}}}{\boldsymbol{\mathsf{F}}^{\dagger}}{\tilde{\boldsymbol{\mathsf{G}}}^{\dagger}}$, where the operator $\tilde{\boldsymbol{\mathsf{Z}}}$ is undoing the zero-padding in the image space and $\tilde{\boldsymbol{\mathsf{U}}}$ is the down-sampling operator.

For image recovery, we adopt the ADMM algorithm solving for the minimisation problem (5). For further details on the ADMM algorithm, we direct the reader to the work of Onose et al. (2016). We consider the following parameters; the parameter controlling the proximal operator of the $\ell_{1}$ norm is set to $\kappa=10^{-3}$ and the stopping criteria which are the maximum number of iterations $I_{\rm max}=2500$, the variation on the solution $\beta=5\times 10^{-6}$. The bound $\epsilon$ on the ${\ell}_{2}$-ball is set from the $\chi^{2}$ distribution of the noise as $\epsilon^{2}=(2M+4\sqrt{M})\sigma^{2}/2$, where $\sigma^{2}/2$ is the variance of both the real and imaginary part of the noise (Carrillo et al., 2012; Onose et al., 2016) and on which we allow an upper bound that is $\epsilon^{\prime}=1.005\times\epsilon$. To assess the reconstruction quality, we adopt two metrics. These are the signal-to-noise ratio (SNR) and the minimum ratio metric (MR), defined as

[TABLE]

and

[TABLE]

where $\bar{\boldsymbol{x}}$ denotes the estimated signal of the ground-truth signal $\boldsymbol{x}$. Compared to the SNR, the MR metric can be more sensitive to pixel-wise errors, in particular for high dynamic range images such as the $\rm{GC}$ image, where the $\rm{SNR}$ metric tends to be dominated by the brightest structures, containing the highest energy of the image.

4.2 Sparsification of the Fourier-transformed Chirp operator

We study the sparsification technique to each row $\hat{\boldsymbol{c}}_{\ell}$ of the operator $\hat{\boldsymbol{\mathsf{C}}}$ preserving the energy $E^{\rm sparse}_{\ell}=(1-{10}^{n_{\rm sparse}})\times E^{\rm{total}}_{\ell}$ with $n_{\rm sparse}\in\{-6,-5,-4.5,-4,-3.5,-3,-2\}$ denoting the sparsification level. Since the chirp kernel is a norm-preserving modulation, $E^{\rm{total}}_{\ell}=1$, $\forall\ell$. We consider various $w$-rates within 5 different ranges, such that $\forall~{}{1\leq\ell\leq M},~{}|w_{\ell}|\leqslant w_{f}\times w^{*}$, with $w_{f}\in\{0.1,0.2,0.3,0.4,0.5\}$. To quantify the sparsity of an operator, we adopt the sparsity ratio (SR) metric that is defined as the ratio between the cardinality of an operator $\boldsymbol{\mathsf{R}}$ and the total number of its elements; ${\rm SR}={\rm{card}(\boldsymbol{\mathsf{R}})}/({N^{\prime\prime}\times M})$.

In Fig. 5, we show the sparsity ratio (SR) of the operator $\hat{\boldsymbol{\mathsf{C}}}$ as a function of the sparsification levels $n_{\rm sparse}$. It is clear that increasing the sparsification levels results in a significant decrease of the operator’s SR for all the considered ranges of the $w$-rates, even though the studied energy loss is, at most, two orders of magnitude smaller than the original energy. Naturally, the sparsity ratio of the $w$-projection $\tilde{\boldsymbol{\mathsf{G}}}$ is of the same order as that of the Chirp operator $\hat{\boldsymbol{\mathsf{C}}}$, given that the de-gridding kernels are of fixed small support $P=4\times 4$.

In Fig. 5, we display the reconstruction quality of the estimated images $\bar{\boldsymbol{x}}$, in terms of SNR and MR as functions of the sparsification levels $n_{\rm sparse}$. Results are shown for the different $w_{f}$ considered, describing the different amplitudes of the $w$ values. For additive noise on the generated visibilities characterised by ${\rm iSNR}=30$dB, the SNR is not affected up to sparsification level smaller than $n_{\rm sparse}^{*}=-3$, independently of the $w$ range. Similar behaviour is observed by the MR metric. Note that, with respect to the sparsity of the $w$-projection operator, obtained with the least energy loss corresponding to $n_{\rm sparse}=-6$, the sparsity ratio up to the level $n_{\rm sparse}=-3$ has decreased by at least $50\%$ for $w_{f}=0.5$ and $83\%$ for $w_{f}=0.1$. This constitutes a significant decrease in memory requirements, leading to a fast application of the measurement operators. For sparsification level $n_{\rm sparse}=-2$ and energy loss on each row $E^{\rm{lost}}_{\ell}=10^{-2}$, both SNR and MR decrease significantly. This is due to the fact that model errors of radio interferometric measurements introduced by the sparsification of the operator $\hat{\boldsymbol{\mathsf{C}}}$ hence the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$, are larger than the noise level on the measurements. For $\rm iSNR=40dB$, the sparsification does not hamper the reconstruction quality up to sparsification level $n_{\rm sparse}^{*}=-3.5$. Note that, while this effect is not reflected in the SNR plot (Fig. 5, upper panel), as the metric is also stable up to $n_{\rm sparse}=-3$, the MR metric (Fig. 5, lower panel) is however stable up to $n_{\rm sparse}^{*}=-3.5$, for higher levels a significant decrease is observed (e.g. $50\%$ decrease for $w_{f}=0.1$).

We have also considered the sparsification of the Chirp operator $\hat{\boldsymbol{\mathsf{C}}}$ for the natural $w$-rates of the ASKAP $uv$-coverage. Most of the $w$ values are extremely small, resulting in flat chirps. Nevertheless, the reconstruction quality is affected for high sparsification levels (see Fig. 5). In fact, for ${\rm iSNR}=30,40~{}\rm dB$ a decrease of the SNR by almost $\rm 2dB$ and $\rm 4dB$, respectively, is observed for the sparsification level $n_{\rm sparse}=-2$, where $83\%$ decrease of the SR of the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$ is reached. This is consistent with our finding at the smallest simulated $w$-rates (e.g. $w_{f}=0.1$). The drop in the image quality is mainly exhibited as artefacts surrounding the faint diffuse central structure.

Large errors on the model of radio interferometric measurements introduced by applying the sparsification of the operator $\hat{\boldsymbol{\mathsf{C}}}$ and consequently the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$ result in estimated images with high artefacts. In fact, the convolution kernels of $\tilde{\boldsymbol{\mathsf{G}}}$, being compact and hence not taking into account energy coming from distant Fourier modes, give rise to over-fitting of the probed Fourier modes. This effect can be shown through the inspection of Fig. 6, where for the sparsification level $n_{\rm sparse}=-3$, the recovered images present strong artefacts consisting in large number of spurious point sources around the faint diffuse central structure. The effect is also reflected on the residual images when inspecting their histograms. More precisely, through the measures of skewness (this measures the lack of symmetry of the histogram) and kurtosis (this quantifies if the distribution is heavy-tailed or light-tailed relative to a normal distribution, in the latter case $\rm kurtosis=3$), we found that the residual images with high sparsification level depart from a Gaussian noise. Moreover, for such high sparsification levels, we noticed that ADMM did not reach the bound $\epsilon$ on the $\ell_{2}$ ball that is set with respect to the noise statistics. Naturally, when model errors are larger than the noise level, solving the constraint version of the minimisation problem (5) results in data over-fitting due to the inaccurate bound on the data fidelity term.

In terms of computational time, the sparsification of the chirp kernels resulting in sparse measurement operators leads to significant acceleration of the image reconstruction step while preserving its quality. For instance, almost $45\%$ drop of the computational time is observed for tests with $w_{f}=0.5$ when considering a sparsification level $n_{\rm sparse}^{*}=-3.5$, i.e., a preserved energy $E^{\rm sparse}=0.9997$ with respect to the computational time for sparsification level $n_{\rm sparse}=-6$ with energy $E^{\rm sparse}=0.999999$. In general, for all the studied $w$ ranges, a decrease of at least $40\%$ of the computational time is observed for the same sparsification level $n_{\rm sparse}^{*}=-3.5$.

4.3 Direct sparsification of the $w$-projection operator

Adopting similar settings to the previous section, we study the sparsification technique to each row $\tilde{\boldsymbol{g}}_{\ell}$ of the operator ${\tilde{\boldsymbol{\mathsf{G}}}}$ preserving the energy $E^{\rm sparse}_{\ell}=(1-{10}^{n_{\rm sparse}})\times E^{\rm{total}}_{\ell}$ with $n_{\rm sparse}\in\{-6,-5,-4,-3.5,-3,-2.5,-2\}$ denoting the sparsification level. We consider various $w$-rates within 5 different ranges, such that $\forall~{}{1\leq\ell\leq M},~{}|w_{\ell}|\leqslant w_{f}\times w^{*}$, with $w_{f}\in\{0.1,0.2,0.3,0.4,0.5\}$. We also consider the original, un-scaled $w$ values given by the baseline distribution.

In Fig. 9, we show the sparsity ratio of the operator $\tilde{\boldsymbol{\mathsf{G}}}$, where a drop is clearly visible for increasing sparsification levels. For instance, considering $w_{f}=0.5$, the sparsity ratio for $n_{\rm sparse}=-6$ is 0.38 and decreases by $57\%$ for $n_{\rm sparse}=-3$. Also in the case $w_{f}=0.1$, the sparsity ratio for $n_{\rm sparse}=-6$ is 0.08, and decreases by $87\%$ for $n_{\rm sparse}=-3$. In Fig. 9, first and second panels, the results are obtained with ${\rm iSNR}=30$dB. For the $w$ ranges considered, the SNR values are unaffected by the sparsification of the operator $\tilde{\boldsymbol{\mathsf{G}}}$ up to $n_{\rm sparse}=-2.5$. The MR is however stable up to $n_{\rm sparse}^{*}=-3$. Note that, the values of the SNR at level $n_{\rm sparse}\geq-2.5$ are biased by the reconstruction of the brightest sources in the GC image. In this case, the MR metric allows for a better judgement of the reconstruction quality. For high sparsification level $n_{\rm sparse}=-2$, a significant drop of both SNR and MR is observed. It is worth mentioning that the high SNR values obtained with $n_{\rm sparse}=-2.5$ with respect to smaller levels are promoted by the over-fitting effect due to sparsification. In fact, this effect allows for a better recovery of the higher spacial-frequency content in the image coming from the brightest compact sources. Yet, this comes at the expense of the recovery of the faint extended emission, which is mainly contributing to the low spacial-frequency content. Results of the sparsification with noise level on the visibilities of $\rm{iSNR}=40$dB are shown in Fig. 9, third and fourth panels. Similar behaviour is obtained and the quality of the image reconstruction is preserved up to the sparsification level $n_{\rm sparse}^{*}=-3.5$.

We apply this sparsification strategy for the natural $w$-rates. In this test case, the sparsification of the operator $\tilde{\boldsymbol{\mathsf{G}}}$ can also act on the de-gridding kernels; for extremely small $w$ values where the chirp kernel is reduced to a Dirac. Note that sparsifying the Chirp operator $\hat{\boldsymbol{\mathsf{C}}}$ as described in the previous section preserves the de-gridding kernels as it is performed before the convolution of the two kernels. Interestingly, when it is the case, we notice that the reconstruction quality is not prone to degradation up to sparsity levels $n_{\rm sparse}^{*}=-4$ for ${\rm iSNR}=40$dB and $n_{\rm sparse}^{*}=-3$ for $\rm iSNR=30$dB while decreasing the SR of the operator $\tilde{\boldsymbol{\mathsf{G}}}$ by $37\%$ and $48\%$ respectively (see Fig. 9). For higher sparsification levels a significant drop in the reconstruction quality is observed, that is translated in the image as large artefacts hampering the recovery of the faint structures, and in particular the non-recovery of the faint and extended central emission.

Regarding the computational time, it is generally significantly reduced with respect to the tests corresponding to the minimal sparsification level $n_{\rm sparse}=-6$. A decrease by around $30\%$ is noticed for $w_{f}=0.1$ and at least $40\%$ for $w_{f}>0.1$ when adopting the sparsification level $n_{\rm sparse}^{*}=-3.5$. No noticeable increase for the test case adopting the original $w$-rates.

When compared with the previous sparsification strategy applied on the Chirp operator, the direct sparsification of the $w$-projection operator leads to lower SR of its sparse approximation obtained by the highest and quality preserving sparsification levels; $n_{\rm sparse}^{*}=-3$ for iSNR$=30$dB and $n_{\rm sparse}^{*}=-3.5$ for iSNR$=40$dB. We also notice a decrease in the computational time for these levels. This suggests that the direct sparsification strategy not only lowers memory requirements but also the computational time as well. In particular, direct sparsification of the $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$ with energy loss fraction $\gamma=10^{-4.5}$ has shown to preserve image reconstruction quality for ${\rm iSNR}$ up to 40dB on simulated visibilities. Therefore, we consider this case as our fiducial setting for a sparse measurement operator in the next section.

4.4 Revisiting the spread spectrum effect on image reconstruction quality

We highlight the general quality enhancement of the reconstructed images by adopting the $w$-projection algorithm in the framework of convex optimisation for increased $w$-rates. To demonstrate the versatility of our study, we have extended our simulation settings by further adopting the image of the galaxy M31 of size $128\times 128$ pixels. We have also considered a simulated $uv$-coverage from the antenna configuration of the Karoo Array Telescope (MeerKAT)888http://www.ska.ac.za/gallery/meerkat/, a precursor of the SKA. The coverage is obtained, 1 hour observation time and time interval $\delta{\rm t}=20$ minutes, pointing at declination $0d0m0s$ and right-ascension $10d0m0s$. The total number of $\boldsymbol{u}$ points is $M=6048$. Both, the image of M31 and the MeerKAT coverage are displayed in Fig. 10. We apply direct sparsification of the $w$-projection operator with a sparsification level $n_{\rm sparse}=-3.5$ in the backward model. Results are displayed in Fig. 11, where the evolution of SNR of the reconstructed M31 (left panel) and GC (right panel) images are shown as a function of the $w$-rates characterised by $w_{f}$. The output SNR improves with the increasing ranges of $w$-rates for all settings. That is in agreement with the previous studies of the spread spectrum effect.

5 Super-resolution promoted by the $w$-term

In absence of $w$-modulation, the new imaging approaches based on compressive sensing and sparse representations, have shown remarkable improvement of the angular resolution of the reconstructed images when compared to the output restored images of CLEAN, these are limited by the instrumental PSF. Thanks to the sparsity-based prior and positivity imposed on the image to be recovered, this class of approaches has shown the ability to go beyond the instrument resolution limit. In presence of $w$-modulation, the resulting spread spectrum effect (Wiaux et al., 2009b; Wolz et al., 2013) allows for probing radio interferometric measurements containing not only information on the sensed Fourier modes of the imaged sky but also information coming from the surrounding Fourier modes induced by the convolution effect of the $w$-modulation. In other words, this means that a radio interferometric measurement at a probed spacial frequency $\boldsymbol{u}$ contains energy that consists, not only of the Fourier component of the imaged sky at that specific spacial frequency $\boldsymbol{u}$, but also of a linear combination of the Fourier components of the imaged sky from the neighbouring spacial frequencies. This suggests that high frequency content of the signal, that is beyond the band limit of the interferometer, can be probed. Intuitively, this might yield to the recovery of super-resolved images, with bandwidth $B^{\prime}_{\rm max}$ larger than the bandwidth of the array, that is given by the maximum projected baseline (see subsection 2.2).

In this section, we study the potential of the $w$-term to allow for super-resolution on two test cases; the first consists of the study of super-resolution on the galaxy cluster image and the second the study of super-resolution allowing for the detection of separated point sources below confusion noise. The latter is a direct effect of the band limit of the radio interferometer. We re-emphasize that, the spread spectrum effect is not exclusive to the $w$-modulations; in fact it is fundamentally induced by DDEs as they are modelled through convolutional kernels in the Fourier domain, similarly to the $w$-term. Yet, unlike the $w$-terms, they are unknown. Therefore, they need to be estimated along with the unknown radio image, resulting in non-convex blind deconvolution problem.

5.1 Super-resolution on the Galaxy Cluster image

We simulate observations using the GC image of size $N=128\times 128$ as a ground-truth image and a realistic simulated $uv$-coverage of ASKAP antenna configuration for a total observation time of 1 hour, pointing at declination $-10d0m0s$ and right-ascension $0d0m0s$ with time spacing $\delta{\rm t}=10$ minutes, leading to $M=3780$ visibilities. Obtained visibilities are contaminated with additive white Gaussian noise with $\rm{iSNR}=40$dB following (10). Note that in this study, we consider simulations of the ground-truth image sensed with a radio interferometric configuration characterised with a band limit smaller than its full bandwidth $B$. Therefore, we generate $w$-modulations so that $B^{\prime}_{\rm max}=\displaystyle\max_{{1\leq\ell\leq M}}{B_{w_{\ell}}}+{\parallel{\boldsymbol{u}_{\ell}}\parallel_{2}}$ does not exceed $B$, hence no up-sampling operator is required.

We perform a family of tests by down-scaling the simulated $uv$-coverage with factors $B_{f}\in\{0.1,0.2,0.3,0.4,0.5\}$. For every test, we set the probed bandwidth $\tilde{B}=B_{f}\times B$, where $B$ is the image full bandwidth. We also generate random $w$ values within ranges characterized by $w_{f}\in\{0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1\}\times w_{f}^{*}\times w^{*}$, where $w^{*}$ is the maximum $w$-rate inducing a chirp modulation with a bandwidth $B_{c}=B$. The value $w_{f}^{*}$ is the $w$-fraction allowing to probe the non-sensed band of the Fourier plane by the down-scaled $uv$-coverage with a factor $B_{f}$. In this case, $w_{f}^{*}$ induces a chirp modulation with a bandwidth $B_{c}=(1-B_{f})\times B$.

Let $\boldsymbol{z}$ denote the Fourier components of the sky $\boldsymbol{x}$ at the non-sensed high frequency modes such that ${\boldsymbol{z}}=[\Re({{\boldsymbol{\mathsf{M}}}{\boldsymbol{\mathsf{F}}}{\boldsymbol{x}}});\Im({{\boldsymbol{\mathsf{M}}}{\boldsymbol{\mathsf{F}}}{\boldsymbol{x}}})]$. $\boldsymbol{\mathsf{M}}\in\mathbb{R}^{N\times N}$ is a diagonal matrix with binary entries selecting the Fourier modes that are part of the non-probed band and $\boldsymbol{z}\in\mathbb{R}^{2N\times N}$, is the concatenation of the real and imaginary parts of the selected Fourier components ${\boldsymbol{\mathsf{M}}}{\boldsymbol{\mathsf{F}}}\boldsymbol{x}$. To quantify the super-resolution effect promoted by the $w$-modulation, we use the SNR metric on $\boldsymbol{z}$ as follows:

[TABLE]

where $\bar{\boldsymbol{z}}=[\Re({{\boldsymbol{\mathsf{M}}}{\boldsymbol{\mathsf{F}}}\bar{\boldsymbol{x}}});\Im({{\boldsymbol{\mathsf{M}}}{\boldsymbol{\mathsf{F}}}\bar{\boldsymbol{x}}})]$ and $\bar{\boldsymbol{x}}$ is the reconstructed image of the sky.

In Fig. 12, we display the evolution of the SNR for the Fourier components $\boldsymbol{z}$ belonging to the non-sensed band in the Fourier domain, as a function of the range $w_{f}$ characterising the $w$-rates. Clearly, the increase of the values of $w$ results in the improvement of the SNR for the various sensed bandwidths characterised by $B_{f}\in\{0.1,0.2,0.3,0.4,0.5\}$; more than 5dB increase of the SNR is observed with the highest $w_{f}=w_{f}^{*}\times w^{*}$, that is allowing to probe the full bandwidth of the image, when compared to no $w$-modulation ($w_{f}=0$). Moreover, the visual inspection of Fig. 13 confirms the super-resolution effect promoted by large $w$-modulations, in particular when inspecting compact sources present in the GC image. In general, the SNR improvement is more noticeable for smaller sensed bandwidths. Note that, for $B_{f}\geq 0.4$, the frequency content of the image under scrutiny at the non-sensed Fourier modes is very small. Interestingly, in absence of the $w$-modulation, the positivity and sparsity priors allow for a small super-resolution effect, that is reflected in the obtained SNRs; 3.3dB and $4.96$dB for $B_{f}=0.4$ and $B_{f}=0.5$, respectively.

5.2 Point-source confusion limit

Radio measurements are prone to point source confusion due to the band limit of the radio interferometer. This can be reflected in the so-called confusion noise, that is a noise-like background consisting of large number of very faint and slightly above the noise level point-like sources and which are below the resolution limit, thus unresolved. The confusion can also occur in the case of bright unresolved point sources lying on a diffuse background, such case can represent a supernova remnant or protostars in a nebula environment. In this section, we study the ability of the $w$-term to promote super-resolution in particular the ability to separate two unresolved point sources, hence decreasing the confusion limit due to the band limit of the radio interferometer.

To do so, as a ground-truth image, we simulate an image of size $N=128\times 128$ consisting of two Gaussian sources with $\sigma=1.5$ pixels on a diffuse background (see Fig. 15, top left). The diffuse background is modelled as a 2D Gaussian source with $\sigma=50$ pixels multiplied with a random Gaussian field and convolved with a highly asymmetric filter. Observations are generated using ASKAP antenna configuration for the same total observation time and direction in the sky described in Sect. 4 and with various time intervals $\delta{\rm t}=10,20,30$ minutes yielding $uv$-coverages with number of measurements $M=3780,1890,1260$ respectively, and which are denoted by settings $\rm A$, $\rm B$ and $\rm C$ respectively. Note that the three $uv$-coverages present the same bandwidth, the difference consists in the lower sampling rates for the higher time intervals.

We vary the distance separating the peaks of the two Gaussian sources of interest and adopt distances $d\in\{6,8,10,12,14\}$ pixels. To produce the band limit confusion scenario, we choose to down-sample the coverages with a factor $B_{f}=0.2$. Hence, in theory, the resolution obtained by such sampling is $\theta\simeq 6$ pixels at half width maximum of the PSF’s primary lobe. The diffuse additive background adds complexity to resolve the two Gaussian sources since it mainly consists of low spacial frequency content, hence its energy dominates the sensed low Fourier modes. We also consider additive white Gaussian noise with ${\rm iSNR}=40$dB on the visibility data. Given this setting, we perform a set of tests while generating random $w$ values within the ranges $w_{f}\in\{0,0.1,0.2,0.5,1\}$. Note that, while $w_{f}=1$ is not a realistic setting, the objective of the study is to show the potential of the spread spectrum effect in promoting super-resolution. Measurements are simulated with the three different coverages explained above.

Results are shown in Fig. 14, where we display the cross-section of the reconstructed images at the central column in order to show the separability of the two Gaussian sources. For each plot, the cross-section is displayed for the three settings. In absence of $w$-modulation ($w_{f}=0$), for the separating distances $d=6,8,10$ pixels, the two sources are completely unresolved for all the adopted coverages. For distance $d=12$ pixels, the coverage $\rm A$ and B, having the largest number of measurements, allow for slightly resolving the two sources. For distance $d=6$ pixels, in all settings, including the presence of $w$-modulations, the two sources are unresolved, yet the higher the $w$-modulation, the more the energy coming from the two sources is concentrated. Furthermore, we notice that for the distances $d=8,10,12$ pixels, the two sources are resolved starting from $w$-modulations with $w_{f}=1,0.5,0.1$ respectively. In general, the more important the $w$-rates, the better is the resolution, this is also highlighted in the reconstructed images in Fig. 15, where estimated images are displayed for $d=12$ pixels and coverage $\rm A$, with $M=3780$ measurements. For separating distance $d=14$ pixels, the two sources are resolved in all the settings.

Increasing the number of measurements can allow for slightly better resolution as seen in the cross-section plots for $d=12$ pixels and $w=0$ (Fig. 14: fourth row); both coverages A and B with the high number of points allow to slightly resolve the two sources w.r.t the the coverage C. Yet, the $w$-modulation promotes higher resolution for relatively small $w$-rates, more significantly, as shown in the plot of the same distance with $w_{f}=0.1$ even with small number of measurements (coverage C). In this case, for all the adopted coverages, the two sources under scrutiny are well resolved. Generally, this test case of point source confusion limit demonstrates how powerful the $w$-term is, in terms of small-scale reconstruction as it allows to resolve sources with much higher resolution from fewer visibilities.

6 Summary and Conclusions

In this study, we have revisited the $w$-modulation in the context of convex optimisation, where it is incorporated in the measurement equation as a sparse thus fast measurement operator. We employ adaptive sparsification strategies on the $w$-projection operator based on the energy of its convolution kernels, and study their effects on the image reconstruction quality. Two main strategies have been investigated. The first consists in the sparsification of the Fourier-transformed chirp kernels prior to the convolution with the de-gridding kernels. The second consists in the direct sparsification of the rows of $w$-projection operator $\tilde{\boldsymbol{\mathsf{G}}}$, these are the convolution kernels coupling the exact chirp kernels and the de-gridding kernels. We find that the latter strategy is prone to more sparsity of the measurement operator while preserving the image reconstruction quality. Furthermore, it is robust when no or negligible $w$-modulations are considered. This suggests that a sparser variant of the de-gridding matrix can be adopted. In general, as long as errors in the model of the visibilities introduced by the sparse approximation of the measurement operator are below the noise level, the reconstruction quality is preserved. In particular, for $\rm iSNR$ up to $\rm 40dB$, sparsification of the measurement operator (following both strategies) with energy loss fraction per row of order $\gamma=10^{-4}$, is robust. Moreover, the sparsity ratio of the sparse measurement operator is reduced by at least a factor 2 in most cases, resulting in significantly lower memory requirements and computational time.

Our C++ code is available online on GitHub, http://basp-group.github.io/purify/. As a future work, we plan to apply these sparsification schemes to real radio interferometric data; where we take advantage of the recent block splitting algorithmic structure proposed in Onose et al. (2016). The ability to split the data and the measurement operator into blocks and process them in a distributed manner is very promising for the applicability of $w$-projection in the context of large-scale data. Therefore, we plan to investigate efficient data block splitting strategies taking into account the $w$-terms. We also consider investigating non-constrained formulation of the radio interferometric imaging problem when adopting very sparse $w$-projection operators, yielding large errors on the backward model, possibly up to the noise level.

We have presented a first study of a very interesting potential of the spread spectrum effect that is super-resolution. We showcased super-resolution through simulated band limited observations and point-source confusion limit experiments. Inaccessible information on the images to be recovered are resolved by adding the $w$-component into the imaging process. We demonstrate how the $w$-term accesses Fourier modes beyond the band limit of the instrument. Important $w$-modulations allow the resolution of small-scale structures inaccessible by the band limit of the radio interferometer, in addition to an extensive enhancement of the image quality. Since the bandwidth of the $w$-modulation depends on the $w$ amplitude and the probed FoV, such potential is highly expected for radio-interferometers characterised with very long baselines and probing large field-of-views, namely the SKA. On the other hand, it is not restricted to the $w$-term. It is inherent to all DDEs, independently of their origin, thanks to their convolution nature in the imaging process. This suggests that our study is relevant for future software developments considering multiple DDEs.

Acknowledgements

We warmly thank Federica Govoni and Matteo Murgia for providing the galaxy cluster simulated radio map analysed in this paper. This work was supported by the UK Engineering and Physical Sciences Research Council (EPSRC, grants EP/M008843/1, EP/M011089/1); the UK Science and Technology Facilities Council (STFC, grant ST/M00113X/1); the Centre for All-Sky Astrophysics (grant CE11E0090) and the gSTAR national facility at Swinburne University of Technology, funded by Swinburne and the Australian Government’s Education Investment Fund.