Multi-GPU maximum entropy image synthesis for radio astronomy

TL;DR

This paper introduces a high-performance multi-GPU implementation of the maximum entropy method for radio astronomy image synthesis, achieving significant speedups over CPU versions while maintaining image quality.

Contribution

It presents a novel multi-GPU version of non-gridding MEM utilizing Peer-to-Peer and Unified Virtual Addressing for efficient large-scale radio interferometry imaging.

Findings

Achieves 1000 to 5000 times speedup over CPU implementations.

Successfully reconstructs complex astronomical data in minutes instead of days.

Demonstrates effectiveness on real ALMA data sets.

Abstract

The maximum entropy method (MEM) is a well known deconvolution technique in radio-interferometry. This method solves a non-linear optimization problem with an entropy regularization term. Other heuristics such as CLEAN are faster but highly user dependent. Nevertheless, MEM has the following advantages: it is unsupervised, it has a statistical basis, it has a better resolution and better image quality under certain conditions. This work presents a high performance GPU version of non-gridding MEM, which is tested using real and simulated data. We propose a single-GPU and a multi-GPU implementation for single and multi-spectral data, respectively. We also make use of the Peer-to-Peer and Unified Virtual Addressing features of newer GPUs which allows to exploit transparently and efficiently multiple GPUs. Several ALMA data sets are used to demonstrate the effectiveness in imaging and to evaluate GPU performance. The results show that a speedup from 1000 to 5000 times faster than a sequential version can be achieved, depending on data and image size. This allows to reconstruct the HD142527 CO(6-5) short baseline data set in 2.1 minutes, instead of 2.5 days that takes a sequential version on CPU.

Tables

Table 1. Table 1: ALMA data sets summary.

Data set	ADS/JAO.ALMA	Beam	SPWs	Channels/SPW	$Z$
	Code
SDP 8.1 Band 7	2011.0.000016.SV	${0.02}^{\prime \prime }\times {0.02}^{\prime \prime }$	4	4	121,322,400
HLTau Band 6	2011.0.000015.SV	${0.04}^{\prime \prime }\times {0.02}^{\prime \prime }$	4	4	96,399,248
HLTau Band 6 SPW 0	2011.0.000015.S	${0.03}^{\prime \prime }\times {0.02}^{\prime \prime }$	1	4	835,360
Antennae North Band 7	2011.0.00003.SV	${1.06}^{\prime \prime }\times {0.68}^{\prime \prime }$	1	1	149,390
HD142527 CO(6-5) Band 9	2011.0.00465.S	${0.23}^{\prime \prime }\times {0.18}^{\prime \prime }$	1	1	107,494

Table 2. Table 2: Average time per iteration (min) and speed factor for single CPU and single GPU versions of MEM, for two data sets and varying image size.

Image Size	Data set	CPU time	GPU time	Speedup
$1024\times 1024$	CO(6-5)	$131$	$0.08$	$1,638$
	HLTau B6w0	$2,645$	$0.60$	$4,408$
$2048\times 2048$	CO(6-5)	$349$	$0.11$	$3,173$
	HLTau B6w0	$10,339$	$2.39$	$4,326$
$4096\times 4096$	CO(6-5)	$360$	$0.13$	$2,769$
	HLTau B6w0	$55,009$	$9.86$	$5,579$

Table 3. Table 3: Occupancy and taxonomy (Gregg and Hazelwood, 2011 ) of the most important kernels using a 1024 thread block size and different image sizes.

Kernel	Taxonomy	Registers per thread	Achieved Occupancy (%)
			1024x1024	2048x2048	4096x4096
KGradPhi	ND	10	77	76	78
KChi2Res	ND	12	81	79	80
KGradChi2	ND	30	100	100	100
KEntropy	ND	13	86	81	79
KGradEntropy	ND	13	87	81	80
KAttenuation	SDH2D/ND	24	89	90	90
KInterpolation	ND	32	92	90	88
KModulation	ND	18	89	89	88
KReduce	SDD2H	13	75	88	96

Equations

$$V(u,v)={\int }_{{\mathbb{R}}^2}A(x,y)I(x,y)e^{-2\pi i(ux+vy)}dxdy$$

$$P(I|V^o)=\frac{P(V^o|I)P(I)}{P(V^o)}$$

$${\hat{I}}_{\text{MAP}}=\arg \underset{I}{\max }P(I|V^o)$$

$${\hat{I}}_{\text{MAP}}=\arg \underset{I}{\max }\Phi (I,V^o)$$

$$\Phi (I,V^o)=\log P(V^o|I)+\log P(I)$$

$$\displaystyle P(V^{o}|I)\propto\displaystyle\prod_{k}\exp\biggl{\{}-{\frac{1}{2}\biggl{|}\frac{V^{m}_{k}(I)-V^{o}_{k}}{\sigma_{k}}}\biggr{|}^{2}\biggr{\}}$$

$$P(I)=\frac{N!}{M^N{\prod }_i{N}_i!}$$

$$\displaystyle\log(P(Ṿ^{o}|I))=$$

$$\log (P(I))\sim S=$$

$$\displaystyle\Phi(I,V^{o};\lambda,G)=\frac{1}{2}\sum_{k}{{\biggl{|}\frac{V^{m}_{k}(I)-V^{o}_{k}}{\sigma_{k}}\biggr{|}^{2}}+\lambda\sum_{i}{\frac{I_{i}}{G}\log{\frac{I_{i}}{G}}}}$$

$$\Phi (I,V^o;\lambda ,G)={\chi }^2(V^o,I)-\lambda S(I;G)$$

$$\displaystyle\chi^{2}=\frac{1}{2}\sum_{w=0}^{W-1}\sum_{c=0}^{C-1}\sum_{k=0}^{Z-1}{{\biggl{|}\frac{V^{m}_{w,c,k}(I)-V^{o}_{w,c,k}}{\sigma_{w,c,k}}\biggr{|}^{2}}}$$

$${\hat{I}}_{\text{MEM}}=\arg \underset{I\ge G}{\min }\Phi (I;V^o;\lambda ,G)$$

$${\sigma }_D=\frac{1}{\sqrt{{\sum }_k\frac{1}{{\sigma }_k^2}}}$$

$$V_{w,c}^F={\mathcal{F}}_{\text{2D}}\{A_{w,c}\cdot I\}$$

$${\hat{V}}^F(u,v)=V^F(u,v)\exp \{2\pi i(u{x}_c+v{y}_c)/M^{1/2}\}$$

$$[\nabla S{]}_i=1+\log \frac{I_i}{G}$$

$$[\nabla\chi^{2}]_{m,n}=\sum_{k=0}^{Z-1}\frac{\text{Re}\bigl{(}V^{R}_{k}e^{2\pi iX_{i}\cdot U_{k}}\bigr{)}}{\sigma^{2}_{k}}$$

$$O=\frac{\text{Number of Active Warps}}{\text{Total Number of Warps}}$$

$$I^D(X)=\sum _k{\omega }_k{V}_k^G\exp (2\pi i{U}_k\cdot X),$$

$$\text{leavevmode}\text{with},{\omega }_k=(1/{\sigma }_k^2)/\sum _l(1/{\sigma }_l^2).$$

$$\text{PSF}(X)=\sum _k{\omega }_k\cos (U_k\cdot X)$$

$${\sigma }_D^2$$

$$\displaystyle=\text{Var}\biggl{(}\sum_{k}{\frac{\omega_{k}}{\sum_{l}{\omega_{l}}}\bigl{(}V^{oR}_{k}\cos(\theta_{k})-V^{oI}_{k}\sin(\theta_{k})\bigr{)}}\biggr{)}$$

$$\displaystyle=\sum_{k}{\frac{\omega_{k}^{2}}{(\sum_{l}{\omega_{l}})^{2}}\biggl{(}\frac{1}{\omega_{k}}\cos^{2}(\theta_{k})+\frac{1}{\omega_{k}}\sin^{2}(\theta_{k})\biggr{)}}$$

$$=\frac{{\sum }_k{\omega }_k}{({\sum }_k{\omega }_k{)}^2}$$

$${\sigma }_D$$

$$\frac{(\Delta x{)}^2}{(\pi /(4\ast \log (2))\ast \text{BMAJ}\ast \text{BMIN}}$$

$$\frac{\partial V_k^R}{\partial I_l}=\frac{\partial V_k^m}{\partial I_l}\sim \frac{\partial }{\partial I_l}\sum _j{I}_j{e}^{-2\pi i{\theta }_{kj}}=e^{-2\pi i{\theta }_{kl}}$$

$$\frac{\partial }{\partial I_l}{\chi }^2$$

$$\displaystyle=\frac{1}{2}\sum_{k}{\frac{1}{\sigma_{k}^{2}}\biggl{(}\frac{\partial\overline{V^{R}_{k}}}{\partial I_{l}}V^{R}_{k}+\overline{V^{R}_{k}}\frac{\partial V^{R}_{k}}{\partial I_{l}}\biggr{)}}$$

$$\displaystyle=\sum_{k}{\frac{1}{\sigma_{k}^{2}}\text{Re}\biggl{(}\frac{\partial\overline{V^{R}_{k}}}{\partial I_{l}}V^{R}_{k}\biggr{)}}=\sum_{k}{\frac{\text{Re}\bigl{(}V^{R}_{k}e^{2\pi i\theta_{kl}}\bigr{)}}{\sigma^{2}_{k}}}$$

Full text

1 Introduction

Current operating radio astronomy observatories (e.g. ALMA, VLA, ATCA) consist of a number of antennas capable of collecting radio signals from specific sources. Each antenna’s signal is correlated with every other signal to produce samples of the sky image $I(x,y)$, but on the Fourier domain (Candan et al., 2000). These samples $V(u,v)$ are called visibilities and comprise a sparse and irregularly sampled set of complex numbers in the $(u,v)$ plane. A typical ALMA sampling data set contains from $10^{4}$ to more than $10^{9}$ sparse samples in one or more frequency channels.

In the case where $V(u,v)$ is completely sampled, Equation 1 states the simple linear relationship between image and data:

[TABLE]

Thus the image can be recovered by Fourier inversion of the interferometric signal (Clark, 1999). In this equation, kernel $A(x,y)$ is called the primary beam (PB) and corresponds to the solid angle reception pattern of the individual antennas and it is modelled as a Gaussian function. If the antennas are dissimilar, it is the geometric mean of the patterns of the individual antennas making up each individual baseline (Taylor et al., 1999).

In the real scenario of collecting noisy and irregularly sampled data, this problem is not well defined (Marechal and Wallach, 2009; Chen, 2011). To approximate the inverse problem of recovering the image from a sparse and irregularly sampled Fourier data a process called Image Synthesis (Thompson et al., 2008) or Fourier Synthesis (Marechal and Wallach, 2009) is used. Current interferometers are able to collect a large number of (observed) samples in order to fill as much as possible the Fourier domain. As an example, Figure 1 shows the ALMA 400 meter short baseline sampling for Cycle 2 observation of the HD142527 plotoplanetary disc. Additionally, the interferometer is able to estimate data variance $\sigma^{2}_{k}$ per visibility as a function of the antenna thermal noise (Thompson et al., 2008).

Many algorithms have been proposed for solving the image synthesis problem and the standard procedure is the CLEAN heuristic (Hogbom, 1974). This algorithm is based on the dirty image/beam representation of the problem (A), which results in a deconvolution problem (Taylor et al., 1999). CLEAN has been interpreted as a matching pursuit heuristic (Lannes et al., 1997), which is a pure greedy-type algorithm (Temlyakov, 2008). Image reconstruction in CLEAN is performed in the image space using the convolution relationship, and it is therefore quite efficiently implemented using FFTs. This algorithm is also supervised. The user could indicate iteratively in which region of the image the algorithm should focus. However, statistical interpretation of resulting images and remaining artifacts are yet to be described by a well founded theory.

MEM is inspired by a maximum likelihood argumentation since interferometer measurements are assumed to be corrupted by Gaussian noise (Sutton and Wandelt, 2006; Cabrera et al., 2008). Reconstructed images with this method have been considered to have higher resolution and fewer artifacts than CLEAN images (Cornwell and Evans, 1985a; Narayan and Nityananda, 1986; Donoho et al., 1992). MEM was the second choice image synthesis algorithm (Casassus et al., 2006), and is mainly used for checking CLEAN bias and images with complex structure (Neff et al., 2015; Coughlan and Gabuzda, 2013; Warmuth and Mann, 2013). However, routine use of MEM with large data sets has been hindered due to its high computational demands (Cornwell and Evans, 1985b; Taylor et al., 1999; Narayan and Nityananda, 1986). For instance, to reconstruct an image of $10^{3}\times 10^{3}$ pixels from the ALMA Long Baseline campaign dataset HL Tau Band 3 (279,921,600 visibilities), MEM computes $7.9\times 10^{16}$ floating-point operations per iteration, approximately. It was also claimed that MEM could not reconstruct point-source like images over a plateau e.g. (Rastorgueva, E. A. et al., 2011). However, we have found (Section 4) this not to be true using a simulated data set. In fact, the experiments show MEM can reconstruct the point source and maintain a smooth plateau.

The MEM algorithm is traditionally implemented with optimization methods (OM) based on the gradient (Cornwell and Evans, 1985b), also called first order methods.

Typical examples of such methods are the conjugate gradient (Press et al., 1992), quasi-newton methods like the modified Broyden’s L-BFGS-B (Nocedal and Wright, 2006), and the Nesterov’s class of accelerated gradient descent (Nesterov, 2004). All of them require the first derivative calculation which has a computation complexity of $O(M\cdot Z)$, where $Z$ is the number of samples and $M$ is the number of image pixels. The gradient computation is the most expensive part per iteration of first order OMs. In consequence, such OMs are equivalent per iteration from the complexity point of view. In particular, we choose to study the GPU implementation for the Conjugate Gradient (CG). This version uses a positive projected gradient due its simplicity for the case of image synthesis and large data ($Z\sim 10^{9}$ and $M\sim 10^{7}$).

A GPU implementation that supports Bayesian Inference for Radio Observations (BIRO) (Lochner et al., 2015) has been proposed before (Perkins et al., 2015). This approach uses Bayesian inference to sample a parameter set representing a sky model to propose visibilities that best match the model. Their model handles point and Gaussian sources. Our approach employs a non-linear optimization problem to directly solve the Bayesian model with a maximum entropy prior using real and synthetic data. However, our method is not currently able to estimate uncertainties by itself as BIRO does.

In this paper we present a high performance, high-throughput implementation of MEM for large scale, multi-frequency image synthesis. Our aim is to demonstrate computational performance of the algorithm making it practical for research in radio-interferometry for today data sets. The main features of the solution are:

GPU implementation: With the advent of larger interferometric facilities such as SKA (Quinn et al., 2015) and LOFAR (van Haarlem et al., 2013), efficient image synthesis based on optimization cannot be delivered by modern multi-core computers. However, we have found that the image synthesis problem, as formulated by MEM, fits well into the Single Instruction Multiple Thread (SIMT) paradigm (Section 2.6), such that the solution can make efficient use of the massive array of cores of current GPUs. Also, our GPU proposal exploits features like the Unified Virtual Addressing (UVA) and the Peer-to-peer (P2P) communications to devise a multiple GPU solution for large multi-frequency data.

• 2.

Parameterized reconstruction: Supervised algorithms could become obsolete in this high-throughput regime in favor of more automatic methods based on a fitting criteria. Although MEM requires a few parameters, it does not require user assistance during algorithm iteration. In this sense, we say MEM is an unsupervised image synthesis algorithm.

• 3.

Non-gridding approach: Gridding is the process of resampling visibility data into a regular grid (Taylor et al., 1999). This task is usually carried out by convolving the data with a suitable kernel. However, the resampling process could produce biased results: loose of flux density (Winkel et al., 2016), aliasing and Gibb’s phenomenon (Briggs et al., 1999). Although this processing reduces the amount of data and the computation time by data averaging, it does to the detriment of further statistical fit. The use of a high performance implementation allows us to process the full dataset without gridding, and still achieve excellent computational performance. However, an interpolation step is still needed to compute model visibilities from the Fourier transform of the image estimate.

• 4.

Mosaic support: Mosaic images in interferometry allow the study of large scale objects in the sky. In this Bayesian approach we fit a single model image to the ensemble of all pointings. The process of image restoration (see Sec.3.3) is then performed on residual images obtained with the linear mosaic formula (Thompson et al., 2008; Taylor et al., 1999).

• 5.

Multi-frequency support: Spectral dependency can introduce strong effects into image synthesis (Rau and Cornwell, 2011). In this implementation we have applied our algorithm to reconstruct ALMA datasets with several channels and spectral windows, but for a single band. Typically, total bandwidth amounts to a maximum of $\sim$10% the central frequency. In these cases, image synthesis can be performed assuming a zero spectral index (Thompson et al., 2008).

2 Method and implementation

This section describes the mathematical formulation of the mono-frequency Maximum Entropy Method (MEM) and the multi-frequency MEM together with a positively constrained conjugate gradient minimization algorithm. Finally, the GPU implementation details are given.

2.1 Data description

Let $V^{o}=\{V^{o}_{k}\}$, $k=0,\ldots,Z-1$, be the (observed) visibility data, let $\{\sigma_{k}\}$ $k=0,\ldots,Z-1$ the estimated deviations of each sample, and let $I=[I_{i}]$, $i=0,\ldots,M-1$. be the image to be reconstructed from $V^{o}$. Notice that $I$ is a regularly sampled, and usually a square image function, whereas $V^{o}$ is a list of sparse points.

Sampled visibilities are grouped in spectral windows. Spectral windows are contiguous spectrum whose frequencies are uniformly spaced and also uniformly divided in channels. Therefore, for the case of multi-spectral, multi-frequency data the observed visibilities are indexed as $V^{o}_{w,c,k}$ where $w$ corresponds to a spectral window, and $c$ corresponds to a channel of the spectral window $w$.

Interferomeric data also contains polarization. We simplify our approach by considering only the case of intensity polarization (I). This simplification is rather general in practice. Radio interferometers such as ALMA, VLA, or ATCA generate data on orthogonal polarizations (either linear XX, YY, or circular LL, RR), which can enter directly in the algorithm as polarization I e.g. (Taylor et al., 1999).

2.2 Bayes formulation and the maximum entropy method

The Maximum Entropy Method (MEM) can be seen as a Bayesian strategy that selects one image among many feasible. Since the data is noisy and incomplete, the solution space is first reduced choosing those images that fit the measured visibilites to within noise level. Among these, the Bayesian strategy selects the one that has a maximum probability of being observed according to a counting rule (Cornwell, 1988; Narayan and Nityananda, 1986; Thompson et al., 2008).

Specifically, let $P(V^{o}|I)$ be the likelihood of observing the data given image $I$ and let $P(I)$ be a prior knowledge of the image based on the chosen counting rule. Then, by Bayes theorem we obtain the a posteriori probability (Pina and Puetter, 1993)

[TABLE]

where $P(V^{o})$ is the normalizing constant. Thus, the Maximum a Posteriori (MAP) image estimate is

[TABLE]

$\hat{I}_{\text{MAP}}$ corresponds to that image with maximum probability of occurring according to the assumed prior knowledge, among those with maximum probability of producing the observed data within noise level. As exponential functions will be used to model the likelihood, it is preferred to work with the logarithm of Equation (3). Dropping term $\log P(V^{o})$, independent of $I$, we obtain

[TABLE]

where

[TABLE]

The likelihood $P(V^{o}|I)$ can be approximated using the fact that visibilities are independent and identically distributed Gaussian random variables corrupted by Gaussian noise of mean zero and standard deviation $\sigma_{k}$. Then, the likelihood can be expressed as

[TABLE]

where $V^{m}(I)$ denotes model visibilities which are functions of the image estimate.

The MEM image prior is assumed to be a multinomial distribution (Gull and Daniell, 1978) of a discrete total intensity that covers the entire image, in analogy with a CCD camera array receiving quantized luminance from the sky (Pina and Puetter, 1993). Let $M$ be the number of image pixels and let $N_{i}=I_{i}/G$ be the quantized brightness (Sutton and Wandelt, 2006) collected at pixel $i$. $G$ is considered as the minimal indistinguishable signal variation (Pina and Puetter, 1993). Then, the prior can be expressed by counting equivalent brightness configurations:

[TABLE]

where $N=\sum_{i}N_{i}$. Under the assumption of a large number of quantas per pixel the Stirling’s approximation ($\log(x!)\sim x\log(x)-x,x\gg 1$), which results in equation (9). There is a subtle difference from derivation found in (Sutton and Wandelt, 2006; Pina and Puetter, 1993) since the quantization scale $G$ is explicit in the entropy (Cabrera et al., 2008). This results in an additional constraint on the sky intensity $I_{i}\geq G$, where $G$ is sufficiently small to reproduce a large number of brightness quantas in signal. Taking logarithms to Equations (6) and (7), we obtain

[TABLE]

Replacing Equations (8) and (9) into Equation (5), we arrive at the following objective function for minimization with $\lambda=1$:

[TABLE]

We recognize Equation (10) as a typical $\chi^{2}$ term plus an entropy regularization term $-\lambda S$ with a penalization factor $\lambda$ as in (11). The previous generalization is a common expression of the objective function for MEM (Narayan and Nityananda, 1986). Other authors argument penalization in equation (10) as an entropy-based distance from $I$ to a prior blank image $G$ (Cornwell and Evans, 1985b).

[TABLE]

For multi-frequency data, the $\chi^{2}$ term is expressed as:

[TABLE]

Finally, the MEM method for image synthesis consists of solving the following non-linear, constrained optimization problem:

[TABLE]

2.3 MEM parameters

Although our implementation of MEM is an unsupervised algorithm it is still dependent on four important parameters that determine properties of the resulting image. These parameters are the entropy penalization factor $\lambda$, the minimal image value $G$, pixel size $\Delta x$, and image size $M^{1/2}$. Even though a complete study on how these parameters affect image properties is beyond the scope of this work, we briefly mention our approach to select reasonable values for them.

In theory, pixel size should satisfy Nyquist criterion, such that image sampling rate should be at least twice as large as the maximum sampled frequency, that is $1/\Delta x\geq 2u_{\text{max}}$. In practice, this is approximated as 1/5 to 1/10 of the Full-Width-Half-Maximum (FWHM) across the main lobe of the dirty beam (see B).

Assuming the primary beam entirely covers the region of interest, it is recommended that the canvas size ($M^{1/2}\Delta x$) has to be larger than the full extent of the main signal. This argumentation is based on the fact that MEM is known to have optimal performance for nearly black images (Donoho et al., 1992).

In the absence of a non-blank prior image, it is assumed that the minimal value a pixel can take ($G$) is at most the thermal noise $\sigma_{D}$ given by (see B)

[TABLE]

However, it is preferred to use a small fraction of $\sigma_{D}$, specially in low signal-to-noise regimes. Notice that when (and if) the restored image is computed, an additional factor of $\sigma_{D}$ is added to the image.

The penalization factor $\lambda\geq 0$ controls the relative importance between data and the entropy term. When $\lambda=0$ the problem becomes a least-square optimization problem. In case of a very small lambda the problem is nearly a least-square with a lower bound constrain ($I\geq G$). But when $\lambda$ increases, less importance is given to data and solution becomes smoother until become a constant equal to $G$. As long as the image synthesis problem has a non-unique solution, each value for $\lambda$ represents a possible solution for the problem. Therefore, MEM allows to explore smoother solutions by changing this single parameter at the cost of degrading residuals ($\chi^{2}$). For example, we start by choosing a very small value of $\lambda\approx 10^{-6}$ which can be considered a near least-square solution. Then, we re-run the program increasing $\lambda$ by a factor of ten, until the user is satisfied. This procedure can be made practical with a fast enough algorithm implementation.

2.4 Objective function evaluation

The minimization problem (13) can be solved by a constrained Conjugate Gradient (CG) algorithm, which repeatedly evaluates the objective function and its gradient to compute the search direction and step size. Evaluation of the entropy term is straightforward, but evaluation of the $\chi^{2}$ term requires some additional processing. The following steps are required for computation of model visibilities $V^{m}$ required for evaluating $\chi^{2}$.

The attenuation image $A_{m,n}$ or primary beam represents the discrete version of the reception pattern of the telescope in Equation 1. This image is generally modeled as a Gaussian and depends on the radio-interferometer (Taylor et al., 1999) scaling linearly with the frequency. For ALMA the FWHM is 21 arcsec at 300GHz for antennas of 12m diameter. 2. 2.

Model visibilities $V^{F}$ on a grid are obtained by applying a 2D Discrete Fourier Transform to the model image. For the case of multi-frequency data, an attenuation matrix $A_{w,c}$ is built for each channel of each spectral window, resulting in a set of model visibilities, $V^{F}_{w,c}$. Thus, the Fourier modulation and the interpolation step are carried out independently for each frequency channel.

[TABLE]

In case of a mosaic measurement, for each pointing in the sky an attenuation image is calculated and centered in the corresponding field of view. 3. 3.

The phase-tracking of the object has a center according to the celestial sphere (Taylor et al., 1999) coordinates. This implies that the image is shifted according to that center. Thus, a modulation factor is applied as follows:

[TABLE]

where $(u,v)$ are the uniformly spaced grid locations, and $(x_{c},y_{c})$ are the direction cosines of the phase-tracking center. In case of mosaic image, data is divided into several pointing in the sky shifting each field of view according to its location. 4. 4.

Model visibilities $V^{m}$ are approximated from $\hat{V}^{F}$ by using bilinear interpolation method. Bilinear interpolation has the effect of local smoothing over a cell grid, thus avoiding overshoot effects at object’s edges (Burger and Burge, 2010). Aliasing is still an issue in this method, but a common way to overcome this problem is by reducing the Fourier grid cell size increasing the number of image pixels. 5. 5.

Residual visibilities $V^{R}=\{V^{R}_{k}\}$ are calculated from $V^{R}_{k}=V^{m}_{k}-V^{o}_{k}$ in order to compute the $\chi^{2}$ term.

It is worth emphasizing that our algorithm computes the error term for each visibility at their exact $uv$ locations, and do not apply any type of Fourier data gridding to reduce the amount of computation.

2.5 Gradient evaluation

Gradient computation is required for first order optimization methods. Notice that before gradient function is called, the objective function is always executed. Thus, residual visibilities $V^{R}$ are available for gradient calculation. As well as the objective function the first computation has to do with the entropy gradient. Therefore, every coordinate of the entropy gradient ($\nabla S$) has a value that follows Equation 14.

[TABLE]

The $\chi^{2}$ gradient term is calculated from the derivative of the Direct Fourier Transform of the image (see C) as shown in Equation (15). The term $X_{i}=(x,y)$ is the image coordinate of the $i$th pixel according to the direction cosines of the phase tracking center, and $U_{w,c,k}=(u,v)$ is the sampling coordinate of the $k$th visibility at channel $c$ of spectral window $s$.

[TABLE]

2.6 SIMT and GPU implementation

The Single Instruction Multiple Thread (SIMT) paradigm is a compute model in which one instruction is applied to several threads in parallel. This means that when a group of threads is ready to run, the SIMT scheduler broadcasts the next instruction to all threads in the group. This is similar but not equal to the Single Instruction Multiple Data (SIMD) paradigm, in which the same instruction is executed over different data, without referring to a particular execution model. This distinction is subtle but important to produce efficient GPU solutions, because an understanding of what is involved in terms of thread execution can make a big difference in the design and performance of a code.

In particular, our approach satisfies the following three general principles for effective and efficient GPU solutions:

Write simple and short kernels to use a small number of registers per thread, increasing warp number per Streaming Multiprocessor. 2. 2.

Write kernels with spatial locality access to increase memory transfer and cache hit rate. 3. 3.

Keep most of the data in global memory to avoid slow transfers from host to device and vice versa.

Although a complete discussion of GPU programming and GPU terminology can be found in (Cheng et al., 2014), we briefly discuss here some of the terms used in this paper. A kernel is the code executed by threads in parallel on a GPU. How and when these threads execute on the GPU depend on how the programmer decomposes the problem into sub-problems that can be solved in parallel on the available GPU resources. The most basic GPU resources are the processing elements or cores. A core can run at most one thread at a time. Cores are grouped into one or more Streaming Multiprocessors (SM), and all threads in an SM execute the same instruction at the same time. From the programming point of view, threads are grouped into blocks, which in turn are organized into a grid of blocks. Different thread blocks can be scheduled to run in different SM, but once a block is assigned to a SM, all threads from that block will execute on that particular SM until they finish. This also means that all block resources, like thread registers and shared memory need to be allocated for the entire block execution. Using too many resources per thread will limit the number of thread blocks that can be assigned to a SM, which in turn will limit the number of active blocks running in the SM.

We have implemented a positively constrained CG algorithm to solve Equation (13). The algorithm’s main iteration loop is kept in host, while compute intensive functions are implemented in GPUs. To minimize data transfer, most vectors and matrices are allocated and kept in device global memory. Only the image estimate is transferred back and forth between host and device. At convergence, model visibilities are also sent back to host to compute final residuals.

As an example, here we show implementation details of the two most compute intensive functions, namely $\Phi$ and $\nabla\Phi$. Algorithm 1 shows the $\chi^{2}$ host function that invokes 1D and 2D kernels to accomplish its goal. Lines 2 to 5 apply correction factors and the Fourier transformation steps. Although not shown here, each kernel invocation has an associated grid on which threads are run. For instance, the attenuation and modulation kernels on line (2) and (4), respectively, use a $M^{\frac{1}{2}}\times M^{\frac{1}{2}}$ kernel grid, while the interpolation kernel at line 5 uses a 1D grid of $T$ threads, where $T$ is the next higher power of 2, greater than or equal to $Z$. All 2D grids are organized in blocks of $32\times 32$ threads, while 1D grids are organized in blocks of $1024\times 1$ threads. The last two steps correspond to the residuals vector computation (line 6) and its global reduction sum (line 7).

Algorithm 2 shows the KChi2Res kernel. First, each thread computes its index into the grid (line 2), and then computes the corresponding contribution to the $\chi^{2}$ term (lines 4 and 5).

Algorithm 3 lists the pseudo-code for the calculation of the gradient from Equation 15, which is the most compute intensive kernel. This kernel uses a $M^{1/2}\times M^{1/2}$ grid. The first lines (2 and 3) compute the index into the 2D grid. Then every thread calculates their corresponding contribution to a cell of the gradient, for every visibility. Recall that sine and cosine functions are executed in the special function units (SFU) of a GPU. Therefore, if Tesla Kepler GK210B has 32 SFU and 192 cores per SM, an inevitable hardware bottleneck occurs.

2.7 Multi-GPU Strategy

The main idea of using several GPUs comes from Equation 12. Since the contribution of every channel to the $\chi^{2}$ term can be calculated independently from any other channel, it is possible to schedule different channels to different devices, and then do the sum reduction in host. Generally, the number of channels is larger than the number of devices and it is therefore necessary to distribute channels as equal as possible to achieve a balanced workload.

Kernels executing on 64-bits systems and on devices with compute capability 2.0 and higher (CUDA version higher than 4.0 is also required) can use Peer-to-peer (P2P) (Cheng et al., 2014) to reference a pointer whose memory pointed to has been allocated in any other device connected to the same PCIe root node. This, together with the Unified Virtual Addressing (UVA) (Cheng et al., 2014), which maps host memory and device global memory into a single virtual address space can be combined to access memory on any device transparently. In other words, the programmer does not need to explicitly program memory copies from one device to another. However, since memory referencing under P2P and UVA has to be done in the same process, OpenMP is used (OpenMP Architecture Review Board, 2015) to create as many host threads as devices to distribute channels.

Figure 2 illustrates the multi-GPU reconstruction strategy for the case of data with four spectral windows and four channels each, on a four device system. A master thread creates a team of four worker threads which are assigned channels in a round-robin fashion with portions of one channel at a time. Channel data is copied from host to the corresponding device by the master thread. Algorithm 4 shows the multi-frequency host function pseudo code that implements this idea. In line 2 a global variable is defined for the $\chi^{2}$ term, and in line 3 the number of threads is set, always equal to the number of devices. Each thread invokes the Chi2 host function (Algorithm 1) with the appropriate channel data, in parallel. Once the function is done, a critical section guarantees exclusive access to workers to update the shared variable. Since the update consists of adding a scalar value to the global sum, we have decided to keep this computation on host.

Computation of $\nabla\chi^{2}$ follows the same logic as before (see Algorithm 5). The major difference now is that the $\nabla\chi^{2}$ is a vector allocated on GPU number zero, and a kernel is used to sum up partial results. Recall that with P2P and UVA, kernels can read and write data from and to any GPU.

3 Experimental Settings

In this section, we test our GPU version of MEM with three synthetic and three ALMA observatory data sets. We use them to demonstrate the effectiveness of the algorithm and to measure computational performance. As we pointed out in the introduction, this paper is not meant to be a comparison between MEM and CLEAN. This task would require an extensive imaging study which is beyond the scope of this work. Nevertheless, we find useful to display CLEAN images as a reference and because it is the most commonly method used today.

3.1 Data-sets

A $1024\times 1024$ simulated image of a $7.5\times 10^{-4}$ (Jy/beam) point source located at the image center and surrounded by a plateau was used to show that MEM can indeed reconstruct isolated point sources. Image pixel size was set to 0.003 arcsec and the plateau’s magnitude was of $3.9\times 10^{-4}$ (Jy/beam). Synthetic interferometer data was created with ft task from CASA, and the HL Tauri $uv$ coverage (see Table 1). Image resolution was measured by the Full-Width-Half-maximum (FWHM) for MEM and CLEAN images. For this particular example, $\lambda=0.1$ was used in MEM.

To further explore image resolution and other properties of MEM, we have created two simulated data sets using task simobserve. This task let us sample an input image using parameters such as antenna configurations, source direction, frequency, to name a few.

Firstly, a $512\times 512$ noiseless image of ten point sources distributed as shown in Figure 3 was simulated at 245 GHz. Point source intensities were randomly assigned normalized values between 0 and 1, and pixel size was set to 0.043 arcseconds. As it can be seen, point sources were located at random positions except for $p_{4}$ and $p_{5}$, which were placed 0.215 arcseconds apart.

Secondly, a $256\times 256$ image of a processed version of galaxy 3C288 (Bridle et al., 1989) with a pixel size set to 0.2 arcseconds, shown in Figure 6(a) was simulated at 84 GHz. The original image is an extended emission radio galaxy at 4.9 GHz with a mix of smooth, low resolution regions together with a small and high-resolution filament. Pre-processing implied clipping to zero all pixel values that were smaller than $10^{-3}$ (Jy/pixel), and adding a 1.0 Jy noise factor to simulated visibilities. Root mean square error (RMSE) of normalized values were measured in a circular region of interest shown in Figure 6(a).

We test our MEM implementation on real data summarized in Table 1. The software CASA (McMullin et al., 2007) was used to reconstruct CLEAN images that were not downloaded from the ALMA Science verification site and also to use the mstransform command to time average specific spectral windows of the HL Tau Band 6 dataset.

The first ALMA data set was the gravitationally lensed galaxy SDP 8.1 on Band 7 with four spectral windows, four channels and a total of 121,322,400 non-flagged visibilities (Tamura et al., 2015). Channel 1 of spectral window 3 was completely flagged, so there were no visibilities from this particular channel. MEM images of $2048\times 2048$ pixels were reconstructed, with $\lambda=0.0$. The CLEAN image was downloaded from the ALMA Science Verification site, and had $3000\times 3000$ pixels.

The second data set was the CO(6-5) emission line of the HD142527 protoplanetary disc, on band 9 with one channel of 107,494 visibilities (Casassus et al., 2015). This small data set was only used for code profiling and measuring GPU occupancy. Larger data sets could not be used to this end because profiling counters (32 bits) overflow with bigger reconstructions.

To test mosaic support the Antennae Galaxies Northern mosaic Band 7 data set was used. This data had 23 fields, 1 channel, and 149,390 visibilities. Once again, images of $512\times 512$ pixels were reconstructed with $\lambda=0.1$ for MEM reconstruction.

Finally, the HL Tau Band 6 long baseline data set that corresponds to the observations of the young star HL Tauri surrounded by a protoplanetary disk (ALMA Partnership et al., 2015; Pinte et al., 2016; Tamayo et al., 2015). Spectral window zero (four channels) of this data set was used to measure speedup factor for varying image size. Data was time averaged on windows of 300 seconds, producing only 835,360 visibilities. To stress multi-GPU compute capacity, a MEM image of $2048\times 2048$ pixels was reconstructed from HL Tau full Band 6 with four spectral windows, four channels each and a total of 96,399,248 visibilities. The CLEAN image was downloaded from the ALMA Science Verification site and had $1600\times 1600$ pixels. For reference, Table 1 lists relevant features of all ALMA data used.

3.2 Computing performance metrics

Computational performance was measured by the speedup factor between one GPU and one CPU (single thread) wall-clock execution time versions of MEM, for three images sizes, namely $1024\times 1024$, $2048\times 2048$, and $4096\times 4096$. Since the number of MEM iterations depended on the data set, and some of the sequential reconstructions took excessively long, the average execution time per iteration, was used as a timing measure. All speedup results were for short-spaced channels data sets. The CPU implementation is based on the conjugated gradient function frprmn() from the Numerical Recipes book (Press et al., 1992). It also employs the Fast Fourier transform, line search minimization, and line bracketing functions from the same source. The GPU platform consisted of a cluster of four Tesla K80, each one with two Tesla GK210B GPUs, with 2496 streaming processors (CUDA cores) and 12 Gbytes of global memory. It is important to highlight that this system had two PCI-Express ports on which two GPU were connected to each port. The CPU version was timed on an Intel Xeon E5-2640 V2 2.0 GHz processor.

Unfortunately, speedup cannot be measured in function of GPU cores. Speed factor is usually measured according to a certain number of processing elements like threads or processors. However, GPU programmers do not have control over the number of cores used in every kernel, but only over the grid dimensions. Grid size affects the streaming multiprocessors (SM) occupancy which is defined as:

[TABLE]

Thus, occupancy measures how efficiently the SM are being used. Any program inefficiencies like inappropriate grid and block sizes, uncoalesced memory access, thread divergence or unbalanced workload will decrease the number of active warps and therefore SM occupancy and speedup. We used the NVIDIA Profiler (NVIDIA Corporation, 2016b) to measure kernel’s occupancy, floating-point operations, compute bound and memory bound kernels.

3.3 Image restoration

Since MEM separates signal from noise (Thompson et al., 2008; Taylor et al., 1999), the CASA package was used to join the MEM model images and their residuals, following the technique of image restoration. Firstly, the MEM model image was convolved with the clean elliptical beam, which effectively changes the MEM image units from Jy/pixel to Jy/beam, and also degrades its resolution. Secondly, a dirty image of each pointing was produced from the visibility residuals, using a user-specified weighting scheme (Briggs et al., 1999). These two last images were added to generate what is called a MEM restored image.

4 Results

4.1 Computational performance

When using CPU timers, it is critical to remember that all kernel launches and API functions are asynchronous, in other words, they return control back to the calling CPU thread prior completing their work. Therefore, to accurately measured the elapsed time it is necessary to synchronize the CPU thread by calling cudaDeviceSynchronize() which blocks the calling CPU thread until all previously issued CUDA calls by the thread are completed (NVIDIA Corporation, 2016a).

Table 2 shows average time per iteration in minutes, for the sequential and single GPU versions of MEM. Notice that for the HL Tau case, execution times increased proportionally with image size. Thus, going from a $1024\times 1024$ to a $4096\times 4096$ image, GPU execution time increased by a factor of 16, approximately. A similar behavior can be observed for the CPU case. However, the results for the CO(6-5) data sets are different, in this case the execution time and image size are not proportional. It must be remembered that this is the smallest data set used and neither CPU nor GPU are stressed to their maximal computational capacity. In any case, time differences are significant. A single iteration for the $4096\times 4096$ HL Tau sequential case took 38.2 days to finish, while the same reconstruction in GPU took only 9.85 minutes. Considering that the algorithm converged in 147 iterations in GPU, it is clear why we could not measure total execution time of the sequential version of MEM. Speedup factors for all cases are presented in Table 2. The smallest and largest speedup factors achieved were 1,638 and 5,579, respectively. The smallest speedup corresponded to the short baseline CO(6-5) with the smallest image size. Full reconstruction took 2.5 days in CPU and only 2.1 minutes in GPU.

Results for multi-GPU reconstruction are displayed in Figure 4. This graph shows average execution time per iteration as a function of number of GPUs, for the two largest data sets. The total number of visibilities of these data sets were similar, but the number of visibilities per channel were different. Also the SDP 8.1 data had one channel with zero visibilities. Both curves display a typical fixed workload timing behavior, where execution time does not decrease linearly as more GPUs are used. Due to the fact that the number of visibilities in each channel varies, some load imbalance among GPUs may be causing this behavior. However, it is well known in parallel computing that fixed workload always has a maximum speedup possible.

For some of the most important kernels, Table 3 lists number of registers per thread and achieved occupancy for the CO(6-5) data set. The KGradChi2 kernel is the only kernel that achieved 100% occupancy for any image size. This is probably due to the fact that it is the most compute intensive of all the kernels and it is able to maintain the GPU busy without many context-switches. Even though the other kernels require a smaller number of registers per thread than the KGradChi2, which means they can have more active warps per SM, they do not achieve maximum occupancy. Another observation which can be made is that occupancy varies differently for different kernels with image size. For instance, the KGradPhi kernel has a consistent 76% to 78% occupancy, but the KReduce kernel occupancy improves with image size, showing that the kernel itself is not big enough to fully exploit the GPU.

The excellent performance results can be explained by the fact that most implemented kernels do not move data to or from device. Table 3 also classifies each kernel according to a taxonomy proposed in (Gregg and Hazelwood, 2011) based on memory-transfer overhead. As it can be seen all kernels except two are of the Non-Dependent (ND) type, which means they do not depend on data transfer between CPU and GPU. The KAttenuation kernel is a Single-Dependent Host-to-Device (SDH2D) type, because it requires data to be moved from host to GPU in order to proceed. But this is only true during the initialization step of the algorithm, after which becomes a ND kernel for the rest of the iterations. For instance, recall Algorithm 2 which is the host function to compute the $\chi^{2}$ term of the objective function. Once visibility data is on GPU, the program proceeds as a pipeline of kernels which never transfer data back to host, until the last reduce kernel is finished. Nevertheless, the data moved from device to host at the end of each iteration is just a scalar variable and only when the conjugate gradient achieves convergence, the final image and model visibilities are transferred to host.

4.2 Image reconstruction

Figure 5 depicts CLEAN and restored MEM image profiles of the isolated point source image. On one hand, CLEAN achieved a higher point source intensity, but on the other hand MEM performed a smoother fit of the plateau. The FWHM was 0.064 and 0.062 arcsec for MEM and CLEAN, respectively. Notice that CLEAN overestimated the point source intensity and MEM underestimated the value. CLEAN tends to fit the point source and plateau overestimating the peak, but MEM is known to lose flux in recovered images (Narayan and Nityananda, 1986).

For the field of point sources (Figure 3), CLEAN delivered uniform resolution across the entire field of view (0.142 arcsecs) and even though MEM ($\lambda=0.005$ ) produced better resolution than CLEAN, values greatly depend on image location. For instance, FWHM for $p_{0}$ and $p_{9}$ were 0.09 and 0.07 arcseconds, respectively, but for point sources near the center like $p_{3}$ and $p_{2}$, FWHM was 0.05. Since the distance between the two center points $p_{4}$ and $p_{5}$ was 0.215 arcsec, resolution at these locations could no be computed in the CLEAN image as there are only 1.5 pixels to fit a Gaussian. MEM instead delivered a FWHM of 0.05 arcseconds for both points. For intensity reconstruction, CLEAN was able to recover over 99.7% of the signal strength in all cases, and MEM recovered differently at different locations. For point sources away from the center, like $p_{0}$ and $p_{1}$ peak recovery was 98% and 99%, respectively, but for point sources near the center, like $p_{2}$ and $p_{3}$, recovery was only 85% and 74%, respectively.

Images for the last simulated data set are displayed in Figures 6(b) and 6(c). It is clearly noticeable that MEM does a better job removing systematic background noise. In particular at the region of interest, RMSE of normalized reconstructions reach 0.0239 for CLEAN versus a closer 0.0159 for MEM. Filament in the middle of the image is better resolved in MEM due to its smaller resolution. This corresponds the same phenomena found previously in the point source field. It is worth mentioning that CLEAN images for the first two data sets were reconstructed without user supervision, but for the extended emission example, there were necessary 10 interactive cycles of 1000 iterations each.

Image reconstruction results for real data are shown in Figure 7. First column displays MEM model images, second column shows MEM restored images, and third column shows CLEAN images. First and second rows corresponds to the SDP 8.1 on Band 7 and HL Tau on Band 6, respectively. As has been said, both are long baseline data sets, and therefore high angular resolution and greater detail are expected. Third row shows the Antennae Galaxies North on Band 7. This dataset has short baselines and is also a mosaic. Model images (first column) show images with only signal as MEM separates it from noise, and in all cases restored images are of lower resolution than model images.

5 Conclusions

We have developed a high performance computing solution for maximum entropy image reconstruction in interferometry. The solution is based on GPU for single channel data and on multiple GPUs for multi-spectral data. The implementation uses a host algorithm that runs the main iteration loop and orchestrates kernel calls. We have decided to write small kernels to improve data locality and minimize thread divergence. Most data is kept on device memory which also minimizes memory moves between host and devices. Overall, we have found that the algorithm renders naturally into a SIMT paradigm and makes it a good candidate for successful GPU implementation.

The resulting code achieves a speedup of approximately 1681 times for our smallest data set, which means that CO(6-5) can be reconstructed in 2.1 minutes, instead of 2.5 days in single CPU. Long baseline HL Tau Band 6, spectral window zero and $1024^{2}$ image size, takes 58 minutes approximately. However for full band, multi-spectral data sets, like SDP 8.1 and HL Tau, non-gridded MEM reconstruction is still a challenging task. For instance, full data-set band 6 HL Tau reconstruction requires 200 minutes per iteration, which is still expensive for routine use with the considered hardware.

Results of image synthesis with MEM has demonstrated that the method removes systematic background noise, has a better resolution, but reaches less intensities at point sources. Moreover, MEM has demonstrate to have less signal variations in flat areas, and having smoother fit than CLEAN on extended structures. This could clearly benefit MEM users by allowing them to identify a finer morphology of the objects, but at the cost of losing signal intensity.

We would like to emphasize that although this work implements the MEM algorithm using the conjugate gradient algorithm, many other first order optimization methods based on gradient can be used. Furthermore, the design permits the use of other well behaved regularization terms to penalize the solution at a similar cost. In consequence we demonstrate that non-linear optimization methods implemented in GPU can be used in practice for research in image synthesis.

Finally, each image synthesis method has its own bias, but making use of many of them could improve the astrophysical analysis of reconstructions. Therefore, a proven high performance MEM implementation could contribute as a tool for image assessment.

The code used in this work is licensed under the GNU General Public License (GNU GPL) v3.0 and is freely available at https://github.com/miguelcarcamov/gpuvmem.

6 Future work

Results are encouraging, but there are still several challenges that need to be addressed. Most of the algorithm’s time is spent on gradient computation, and for nearly-black objects, evaluation on background pixels far from the source carry little or zero contribution to the signal. Therefore, we are currently modifying our code to include a focus-of-attention improvement which will allow the user to select a region on the image where reconstruction should take place, thus reducing substantially the number of free parameters. We are also working on a solution to reconstruct data sets with long-spaced multiple channels, also called multi-frequency synthesis. Finally, we are now in a position to conduct a throughout experimental evaluation of MEM imaging parameters including resolution, signal-to-noise ratio, signal recovery, and associated uncertainties for real and simulated data sets.

Acknowledgements

This paper makes use of the following ALMA data: ADS/JAO.ALMA #2011.0.000015.SV, ADS/JAO.ALMA #2011.0.00465.S, ADS/JAO.ALMA #2011.0.00003.SV, ADS/JAO.ALMA #2011.0.000016.SV, ADS/JAO.ALMA #2013.1.00305.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. Also, the calculations used in this work were performed in the Brelka cluster, financed by Fondequip project EQM140101 and housed at MAD/DAS/FCFM, Universidad de Chile. F. R. Rannou and M. Carcamo were partially funded by DICYT project 061519RF, Universidad de Santiago de Chile. P. Roman acknowledges postdoctoral FONDECYT projects 3140634, Basal PFB-03 Universidad de Chile, FONDECYT grant 1171841, and Conicyt PAI79160119 Universidad de Santiago de Chile. S. Casassus acknowledges support from Millennium Nucleus RC130007 (Chilean Ministry of Economy), and additionally by FONDECYT grant 1171624.

Appendix A Beam size

Provided some degree of $uv$-coverage, a zero-order approximation to the sky signal can be obtained from the dirty image

[TABLE]

where $V^{G}_{k}$ are the gridded visibilities (Briggs et al., 1999).

$I^{D}(X)$ is called the natural-weights dirty map and has units of Jy beam*-1*, where the beam is the solid angle subtended by the FWHM of the Point Spread Function (PSF)

[TABLE]

The beam is usually approximated by fitting an elliptical Gaussian on the PSF shape and considering the FWHM. The parameters of the ellipse are BMAJ: Major diameter, BMIN: Minor diameter, and BPA: anti-counter clockwise angle in degree. This shape is interpreted as a natural resolution for the CLEAN algorithm.

Appendix B Thermal noise

Under the assumption of uncorrelated samples111visibility data on $uv$ points closer than the size of the dishes are correlated we can obtain an analytic expression for the thermal noise $\sigma_{D}$ using the following error propagation from dirty image $I^{D}$

[TABLE]

where $\theta_{k}=U_{k}\cdot X$, $\omega_{k}=\frac{1}{\sigma^{2}_{k}}$. We also assume uncorrelated imaginary and real part of the samples. $\sigma_{D}$ units are Jy beam*-1*. Then, the constant to change units from Jy beam*-1* to Jy pixel*-1* is

[TABLE]

where $\Delta x$ is the pixel side size in radians.

Appendix C $\chi^{2}$ gradient

The gradient of a model visibility is approximated from the DFT of model image.

[TABLE]

Where $\theta_{kl}=U_{k}\cdot X_{l}$, from previous approximation we obtain

[TABLE]