© ESO, 2016

1. Introduction

Radio interferometers measure the spatial coherence of the electric field (Thompson et al. 2001) in the Fourier domain. Owing to incomplete sampling in the Fourier plane, the telescope’s point spread function (PSF) has wide-spread sidelobes. The true sky brightness is convolved with this PSF, which limits the imaging to no better than a few 100:1 in the dynamic range.

The CLEAN algorithm (Högbom 1974) and its variants (Clark 1980; Schwab & Cotton 1983) have so far been the most popular algorithms in use to remove the effects of the PSF. Most of these algorithms model the sky brightness as a collection of delta functions, with an implicit assumption that the sky is composed of several well-separated point sources. As such, this approach works well for fields dominated by unresolved sources. With the increase in sensitivity of modern radio telescopes, low-level extended emission is often detected in many fields, particularly at low frequencies. Accurate deconvolution of extended emission is an important problem, especially with the sensitivity of new telescopes. Recently, scale-sensitive methods have been developed to mitigate artefacts induced by the original CLEAN algorithm when used to deconvolve extended emission (Bhatnagar & Cornwell 2004; Cornwell 2008; Rau & Cornwell 2011). These scale-sensitive algorithms have been shown to have better performance than scale-insensitive algorithms in recovering extended emission.

The Multi-Scale CLEAN (MS-Clean) algorithm (Cornwell 2008) uses a matched-filtering technique to find the components. It decomposes a sky image into a set of fixed scales, such as tapered and truncated paraboloids. The Asp-Clean algorithm was proposed by Bhatnagar & Cornwell (2004); here the scales are determined adaptively through a fitting procedure and heuristics are developed to limit the set of scales (active set) that need to be fitted during any given iteration. Compared to MS-Clean, the Asp-Clean algorithm gives better imaging performance, but at the cost of a significant increase in computing time.

In this paper, we propose an efficient implementation that significantly reduces the computational cost while keeping the imaging performance comparable to the original Asp-Clean algorithm (referred to as Asp-Clean2004 in the rest of the text). This is possible because the PSF sidelobes of modern wide-band telescopes with excellent coverage of the uv-plane for continuum imaging are significantly reduced. This allows us to approximate the PSF as a single Gaussian. Then the convolution in the objective function of the iterative fitting (using Levenberg-Marquardt (LM) minimization, Marquardt 1963) can be removed. These convolutions are performed through expensive fast Fourier transformation (FFT)-based approach in the Asp-Clean2004 algorithm. With typically large numbers of minimization iterations required, FFT-based convolution was responsible for the high computing cost of the Asp-Clean2004 implementation. Our current approach, referred to as Asp-Clean2016 in the text below, can be thought of as an extreme case of the approximations used in the Clark CLEAN (Clark-Clean) algorithm (Clark 1980), where the full PSF is approximated by a PSF-patch that only includes the highest sidelobe. As the highest sidelobe becomes increasingly smaller with modern wide-band telescopes, Asp-Clean2016 ignores the sidelobes of the PSF entirely. We show that this still leads to convergence and the total runtime is significantly reduced, even if it is at the expense of a slightly larger number of iterations.

We recall the basics of the Asp-Clean2004 algorithm in Sect. 2, while the Asp-Clean2016 algorithm is described in Sect. 3. Numerical experiment and comparisons are presented in Sect. 4, and we summarize our results in Sect. 5.

2. Asp-Clean2004 algorithm

In the image reconstruction problem, the dirty image I^dirty can be expressed as (1)where B is the PSF, I^true is the true image, I^noise is the noise convolved by the PSF and ⋆ is the convolution operator. One goal of imaging algorithms is to remove the effects of the PSF. In practice, a model with a limited number of components is used to approximate the true image, (2)where the model image , is the ith image component, ϵ is the error between the true image and the model image, and L is the number of components used to approximate the true image. In scale-insensitive algorithms, is a delta function. The true sky, including any extended emission, is thus modelled as a collection of delta functions, which leads to residuals at levels much higher than the thermal noise limit of modern telescopes, as shown by various authors (including Bhatnagar & Cornwell 2004). The primary limitation of the imaging performance of scale-insensitive algorithms is that their models cannot represent the correlation between pixels containing extended emission. This is in conflict with the fundamental assumption underlying scale–insensitive algorithms that each pixel in the image must be independent. The Asp-Clean2004 algorithm tries to solve the deconvolution problem by adaptively determining the scales in various parts of the image by fitting the largest possible scale locally using the following procedure:

• 1.

  Determine the initial values of the scale by smoothing theresidual image with a few Gaussian beams, then identify the peakamong the smoothed images and use the widthof the Gaussian as the initial guess.

• 2.

  Refine the initial values with the LM minimization algorithm to determine the best model, , that fits the data locally.

• 3.

  Update the residual image as .

• 4.

  Repeat steps 1–3 till one of the stopping criteria is satisfied. The stopping criteria can be the total number of iterations or a estimated noise threshold.

As mentioned above, step 2 makes the component scales adaptive (i.e. they are not fixed and pre-determined), while the component scales are fixed and pre-determined in MS-Clean and its variants. The minimization algorithm in step 2 minimizes the objective function χ², which is given by (3)where ∥ ∥ ₂ is the Euclidean norm. In Eq. (3), there is a convolution operation between B and . In a typical function minimization, the objective function is evaluated many times for each fitted scale. The FFT-based convolution used in the Asp-Clean2004 algorithm dominates the computational cost. For an N-pixel image, the fast convolution includes two FFTs, one iFFT, and N multiplication, where the computational complexity of an FFT is about Nlog ₂N (Gonzalez & Woods 2010). As such, the computational complexity of a convolution is N(3log ₂N + 1) for an N-pixel image. For M evaluations of the objective function, the computational complexity is MN(3log ₂N + 1) when solving for a single component. While the Asp-Clean2004 reconstruction gives good results, it is limited by the high computational cost.

3. Asp-Clean2016: An efficient implementation of Asp-Clean2004

As is well known,a two-dimensional Gaussian function g₁(a₁,x₁,y₁,ω₁) that is convolved with another two-dimensional Gaussian function g₂(a₂,x₂,y₂,ω₂) results in a new two-dimensional Gaussian function g₃(a₃,x₃,y₃,ω₃), where a_i,x_i,y_i and ω_i are the amplitude, position (x_i,y_i), and width of Gaussian function g_i, respectively. The Fourier convolution theorem states that if the parameters of the convolved result (e.g. g₃) and a Gaussian function (e.g. g₁) are known, the parameters of another Gaussian function (e.g. g₂) can be analytically computed.

We therefore use the objective function (4)where is a Gaussian containing a model component whose parameters (amplitude a_ic, location x_ic,y_ic and width ω_ic) needs to be determined by optimization (using the LM minimization). After converging to the solution of , the underlying component that models the emission is determined by analytically deconvolving the PSF (approximated as a single Gaussian) as (5)where ω_i, ω_ic, ω_b are the widths of , and the main lobe of the PSF, respectively. The amplitude α_i that corresponds to the is calculated by the equation (6)where α_b and α_ic are the amplitudes of the Gaussian beam approximated from the PSF and , respectively. As in other algorithms, a loop gain is used before the component is added to the model image. In Eqs. (5)and (6), the α_b and ω_b of the PSF are fixed, and changes of α_i and ω_i follow from any adjustments of α_ic and ω_ic. Here we would like to point out similarities between our approximation of the PSF by a single Gaussian (effectively ignoring all levels of the PSF sidelobes) and use of an approximate PSF where only the highest near sidelobe of the PSF is included in a similar step as in Clark-Clean. The difference is that our algorithm uses this approximation to the extreme in the estimation phase of the model components, while Clark-Clean uses it in the subtraction step in a minor cycle to determine the next highest peak in the residual image.

Initial values are important to compute an optimal component. In many tests, we found that good initial values that are closer to the optimal component can significantly reduce the iteration counts of the minimization algorithm to find the optimal component. In other words, good initial values can significantly reduce the computational time. As in the Asp-Clean2004 algorithm, good initial values can be obtained by smoothing the current residual image first and then finding the parameters that correspond to the global peak among the smoothed residual images. Smoothing is computed by a convolution, which has the effect of increasing the computational cost. However, when deconvolving PSF-sized features, the peak of the last residual image and the scale of the main lobe of the PSF can be directly used as the initial values of the parameters. Then the simple method can be triggered to reduce the computational cost. In the Asp-Clean2016 algorithm, the simple method is used only when the area of the optimal component in the last iteration is smaller than 5% of the area of the initial guess in the last iteration. After this, the smoothing method is used again until the condition of using the simple method is met. In practice, there are a large number of small scales in image decomposition. The modified scheme not only retains the reconstruction quality, but also reduces the computational cost.

4. Numerical experiment and comparison

We tested the Asp-Clean2016 algorithm and compared it with the Asp-Clean2004 algorithm using simulated L-band EVLA data in the B-configuration¹, with a bandwidth of 1 GHz and 32 channels. The observation time was six hours. All simulations were made using the CASA software². The robust-weighted PSF is displayed in Fig. 1. The maximum negative sidelobe of the PSF was −0.068 (the peak value of the PSF is normalized to 1). The full width at half maximum (FWHM) of the main lobe of the PSF was about 2′′, while the resolution of the images was 1′′.

Fig. 1 Robust-weighted PSF in the range −0.068 to 1.0 with the logarithmic scaling (CASA scaling power cycles =−1.4).

Fig. 2 Results of the M 31 image from the Asp-Clean2016 algorithm: a) the original image; b) the dirty image with a robust weighting function; c) the reconstructed model image; and d) the corresponding residual image.

We tested the Asp-Clean2016 algorithm on the M 31 image with robust weighting applied. The reconstructed results are displayed in Fig. 2. The original image is shown in Fig. 2a, while the dirty image is displayed in Fig. 2b. Artificial Gaussian noise was added to the simulated visibilities, which resulted in a noise RMS of 5×10^-5 Jy in the dirty image. The loop gain used in the reconstruction was 0.35. Comparing the original image in Fig. 2a, the reconstructed image, which has about 900 components (as shown in Fig. 2c), can reconstruct various scales well.

Fig. 3 Models and residuals of the M 31 image from the Asp-Clean2004 and Asp-Clean2016 algorithms: a) the Asp-clean2004 model image; b) the Asp-Clean2016 model image; c) the Asp-Clean2004 residual image; and d) the Asp-Clean2016 residual image. The model images and the residual images are displayed in the same respective data ranges.

We found that the reconstructed results of Asp-Clean2016 and Asp-Clean2004 are comparable (e.g. see Fig. 3). In Fig. 3 the fundamental reason for the difference between the model (a) and (b) is that image reconstruction in interferometric imaging is non-unique. In other words, infinite numbers of image models will fit the data equally well. Two different algorithms will rarely produce the same set of image model components. However, they both should fit the data well. In Tables 1 and 2 we show that the amplitude range and the total reconstructed model flux from the Asp-Clean2016 algorithm is closer to the original image. However, the Asp-Clean2004 residual image is more noise-like and includes fewer correlated structures. From the histogram (see Fig. 4) of the residual images, we can see that the residual noise profiles of Asp-Clean2004 and Asp-Clean2016 are very close to each other and also very close to the input Gaussian noise. This shows that the new algorithm improved the computing efficiency without significant degradation in reconstruction quality. In all, the results of the Asp-Clean2016 results are comparable to those of the Asp-Clean2004 algorithm, but the Asp-Clean2016 algorithm had a faster convergence rate. Table 3 shows that the Asp-Clean2016 algorithm converged about 20 times faster in the test.

To compare the performance of the Asp-Clean2004 algorithm and the Asp-Clean2016 algorithm with different weightings, we simulated the M 31 image with uniform, natural, and robust weightings, respectively. The image sizes were the same: 256×256 pixels. The loop gain was 0.35 and other parameters were the same. The numerical results are included in Table 3. Different weightings in the spatial frequency domain result in different PSFs in the spatial domain, which in turn affects structures in the dirty image. However, the Asp-Clean2016 requires approximating the PSF as a Gaussian function. Obviously, the approximation is more effective when the PSF includes weaker sidelobes. Thus, for the Asp-Clean2016 algorithm, the convergence speed with a uniform weighting was faster than with a natural weighting under the same parameters. This is also partly because the performance of the minimization algorithm differs between dirty images. Table 3 shows that the total runtime of Asp-Clean2004 is 20 or more times longer than that of Asp-Clean2016. The mean one-iteration runtime ratios, which are defined as the ratio of total runtime and total numbers of iterations, are higher than 40 in this test. In other words, the mean one-iteration runtime of Asp-Clean2004 is 40 times longer than that of Asp-Clean2016. When reconstructing the 512×512 M 31 image, the Python code for Asp-Clean2016 took several minutes using a typical computer (Intel(R) Core(TM) i7-3770 CPU @3.4 GHz, 4.00 GB RAM). After it is implemented with C++, the runtime is expected to be much shorter.

Fig. 4 Histogram of the residual images. “FittedOri” is a Gaussian fitted to the “Original” image noise. “Original” is the added image noise, which is the inverse Fourier transformation of the visibility noise. “Asp04” is the Asp-Clean2004 residual image, and “Asp16” is the Asp-Clean2016 residual image.

5. Summary

In this paper, we approximated the PSF with a single Gaussian, which allowed us to analytically compute the optimal components after a Gaussian was fitted to the current residual image. This greatly reduced the total runtime compared to the 2004 implementation of the Asp-Clean algorithm. We also showed through simulations that there was no significant degradation in the imaging performance when using such an approximate PSF. The approximation we made will only become less of a problem with data produced by future telescopes with tentimes more antennas than current telescopes. This will further reduce the PSF sidelobes, which justifies this analytical approximation of the PSF. In addition, we used a different scheme to determine the initial values of the optimal Gaussian components, which were then passed to the minimization algorithm. In the 2004 implementation, initial values were found by smoothing of the residual image with a few Gaussians of different sizes. The heuristic used in the 2016 implementation can switch between smoothing and the simple method to further reduce the computational cost. Tests showed that Asp-Clean2016 converges faster than Asp-Clean2004, especially for larger images. The current implementation of Asp-Clean2016 is in Python and the CASA package. The Python source code is available on the internet³. Its implementation in the C++ in the CASA package is currently underway.

Acknowledgments

We would like to thank the people who develop Python and CASA, which provided an excellent development and simulation environment. L. Zhang thanks for support from the NRAO Graduate Student Internship program. The work was also supported by the National Basic Research Program of China (973 program: 2012CB821804 and 2015CB857100), the Nation Science Foundation of China (11103055) and the West Light Foundation of the Chinese Academy of Sciences (RCPY201105).

References

All Tables

All Figures

	Fig. 1 Robust-weighted PSF in the range −0.068 to 1.0 with the logarithmic scaling (CASA scaling power cycles =−1.4).
In the text

	Fig. 2 Results of the M 31 image from the Asp-Clean2016 algorithm: a) the original image; b) the dirty image with a robust weighting function; c) the reconstructed model image; and d) the corresponding residual image.
In the text

	Fig. 3 Models and residuals of the M 31 image from the Asp-Clean2004 and Asp-Clean2016 algorithms: a) the Asp-clean2004 model image; b) the Asp-Clean2016 model image; c) the Asp-Clean2004 residual image; and d) the Asp-Clean2016 residual image. The model images and the residual images are displayed in the same respective data ranges.
In the text

	Fig. 4 Histogram of the residual images. “FittedOri” is a Gaussian fitted to the “Original” image noise. “Original” is the added image noise, which is the inverse Fourier transformation of the visibility noise. “Asp04” is the Asp-Clean2004 residual image, and “Asp16” is the Asp-Clean2016 residual image.
In the text