Parameterized CLEAN Deconvolution in Radio Synthesis Imaging

The Högbom algorithm (Högbom 1974) is much simpler than the CLEAN framework. In addition to calculating the dirty image, this algorithm only executes the minor cycle in the CLEAN framework. The residual image is updated by subtracting a scaled full PSF containing all sidelobe patterns to remove the contribution of each component. It is generally believed that the closer the residual is to the noise, the closer the model image is to the potential true image. Finding the model is equivalent to minimizing the following objective function,

$$\begin{eqnarray}&&{\chi }^{2}={\parallel {I}_{i-1}^{\mathrm{residual}}-{B}^{\mathrm{dirty}}{I}_{i}^{\mathrm{comp}}\parallel }_{2}^{2},\end{eqnarray} \tag{ 14 }$$

where $\parallel \,{\parallel }_{2}$ is the Euclidean norm. This algorithm is based on the steepest descent method (Bhatnagar & Cornwell 2004). From the viewpoint of function optimization, the first derivative needs to be calculated to find the update direction of each iteration.

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\chi }^{2}}{\partial {p}_{i}}\equiv \displaystyle \frac{\partial {\chi }^{2}}{\partial {I}_{i}^{\mathrm{amp}}}=-2{\left[{I}_{i}^{\mathrm{residual}}\right]}^{T}\displaystyle \frac{\partial {p}_{i}}{\partial {I}_{i}^{\mathrm{amp}}}=-2{I}_{i,{x}_{i},{y}_{i}}^{\mathrm{residual}},\end{eqnarray} \tag{ 15 }$$

where ${p}_{i}={I}_{i}^{\mathrm{amp}}$ and ${I}_{i,{x}_{i},{y}_{i}}^{\mathrm{residual}}$ is the peak of the current residual image located at $({x}_{i},{y}_{i})$ in the ith iterations. The peak search of the CLEAN algorithm is not explicitly fitting, but rather it is equivalent to minimize the objective function (Equation (15)) along the gradient direction  by iteratively peak searching and subtracting the scale PSF. This is achieved through the first three steps of the CLEAN framework. The result is a series of delta functions with the amplitudes and positions of the peaks in each residual image that make up the model image.

It is well known that the magnitude of the Fourier transform of a delta function is a constant  that fills the entire (u, v) plane. This also causes the model's visibility function to be extrapolated out of the measured spatial frequency range. Although a set of delta functions can partly constrain the spatial frequency of the extrapolation, the use of a delta function model to represent extended sources is not physically accurate. This can easily lead to over-interpretation outside the measured spatial frequency range, which is one of the main reasons for using a CLEAN beam to convolve the model.

When deconvolution converges, it is usually assumed that the image mainly contains noise. But when dealing with extended sources, deconvolution usually leaves some signals in the residual image. Therefore the final residual image is often added to the final restored image.Which not only ensures the integrity of the signal in the restored image, but also indicates the level of uncertainty in the restored image.

This algorithm is only performed in the image plane and require a dirty beam image covering twice the extent in each dimension of the image to be deconvolved. Although the calculation is efficient, errors arising from gridding and the fast Fourier transform (FFT) algorithm in the dirty image cannot be corrected.

The shift-and-subtraction operation is very effective for calculating the convolution of point sources. However, it is well known that using the convolution theorem for non-point source convolutions is more efficient. In order to reduce the computational complexity of the situation where a large number of delta functions are required to approximate the potential true image, Clark (1980) proposed a new CLEAN algorithm that uses the convolution theorem and FFTs.

The Clark algorithm consists of a minor cycle and a major cycle. Its minor cycle is similar to that in the Högbom algorithm, but the difference is that a beam patch is used instead of the full beam. A beam patch is composed of the main lobe and the maximum sidelobe of the full beam. In the minor cycle, the residual image is updated as follows,

$$\begin{eqnarray}&&{I}_{i}^{\mathrm{residual}}={I}_{i-1}^{\mathrm{residual}}-{B}^{\mathrm{patch}}{I}_{i}^{\mathrm{comp}},\end{eqnarray} \tag{ 16 }$$

where B^patch is the beam patch approximated from the dirty beam. The criteria for stopping the minor cycle and starting the major cycle are based on the maximum peak outside the beam patch or a certain number of components being collected. The full beam will be used in the major cycle to remove the effects of the set of point sources identified in the minor cycle from the residual image.

$$\begin{eqnarray}&&{I}_{j}^{\mathrm{res}-\mathrm{maj}}={I}_{j-1}^{\mathrm{res}-\mathrm{maj}}-{F}^{-1}\left(F\left({B}^{\mathrm{dirty}}\right)F\left({I}_{i}^{\mathrm{comp}}\right)\right),\end{eqnarray} \tag{ 17 }$$

where B^dirty is the full dirty beam. Here the convolution theorem and the FFT are used to  quickly calculate the convolution. The full beam also corrects some of the errors caused by the beam patch approximation in the minor cycle to some extent. This algorithm is several times faster than the Högbom algorithm.

The Clark algorithm uses the convolution theorem and FFTs to quickly calculate the convolution, and transforms the model image and the dirty beam to the gridded (u, v) plane. These computations greatly accelerate the deconvolution process. At the same time, it also corrects the errors caused by the minor cycle. However, this algorithm still cannot correct the errors caused when calculating the dirty image.

In order to solve this problem, the Cotton–Schwab algorithm (Schwab 1984) was proposed, which is a variant of the Clark algorithm. It performs a major-cycle subtraction on the measured visibility data points,

$$\begin{eqnarray}&&{I}_{j}^{\mathrm{res}-\mathrm{maj}}={F}^{-1}({V}^{\mathrm{measured}}-{{AI}}_{j}^{\mathrm{model}}),\end{eqnarray} \tag{ 18 }$$

where I_j^model is the model image in the jth major cycle, and A maps the model image from the minor cycle onto ungridded measured data points. This algorithm is directly fitted to the measured spatial frequency domain.

$$\begin{eqnarray}&&{\chi }^{2}={\parallel {V}_{i-1}^{\mathrm{residual}}-{{AI}}_{j}^{\mathrm{model}}\parallel }_{2}^{2}.\end{eqnarray} \tag{ 19 }$$

In other words, the model prediction is compared with the original data in the data measurement domain. This significantly reduces aliasing noise and gridding errors. Another obvious advantage of this algorithm is that it can separate CLEAN in each field in the minor cycle, and thus the effects of the components of all fields are removed together (Rau 2010). The framework of this algorithm, which uses minor cycles to determine models in the image domain and major cycles to correct errors in the visibility domain, is effectively the same as the standard framework employed by most CLEAN-based algorithms.

The algorithms mentioned above all have been implemented in major radio synthesis imaging software programs such as CASA.⁵ These scale-free algorithms maintain the idea of the original CLEAN algorithm, which parameterize the radio sky into a series of scaled delta functions. Although experience has shown that these scale-free algorithms also work for extended and complex sources (Taylor et al. 1999; Thompson et al. 2017), as mentioned previously, a widely known problem is that they tend to produce some spurious artifacts such as stripes in the model image when dealing with extended sources. Clark (1982)  provided a heuristic explanation of this issue. When scale-free algorithms remove a point source response (a scaled PSF) from the peak of an extended source, the negative sidelobes of the dirty beam generate new extreme values in the following residual image. When searching for the next component, the peaks will be identified as model components in subsequent cycles. The dirty beam pattern  is therefore in the residual image of each iteration, resulting in the same or similar fringe spacing in the model image  due to the sidelobes of the dirty beam. As the sidelobe level increases, the stripes become more pronounced.

There are many ways to solve this problem. As this problem is caused by the sidelobes of the dirty beam, the general idea is to start with the sidelobes. For example, a small loop gain can be used so that the residual image is reduced by the sidelobes, thereby suppressing the formation of stripes in the model image. As mentioned above, the method (Cornwell 1983) is to add an artificial "spike" at the center of the main lobe of the dirty beam to compress the sidelobes  and in turn reduce/remove the stripes from the model image. Another idea is to directly suppress the stripes in the model image (Palagi 1982; Cornwell 1983). For example, Cornwell (1983) added a smooth regularization term based on model images to the objective function to achieve the suppression of stripes. Steer et al. (1984) attempted to solve this problem by using a group of delta functions as model components above the threshold of the current residuals. A more effective method is to use a scale-sensitive basis to parameterize the radio sky, which is the subject of the next two sections.

In this algorithm, the radio sky is parameterized as a combination of some tapered and truncated parabola with limited and different spatial scales. The component shapes are controlled by the following function f (Cornwell 2008),

$$\begin{eqnarray}&&f(r,s)={\rm{\Phi }}(r)\left(1-{\left(r/s\right)}^{2}\right),\end{eqnarray} \tag{ 21 }$$

where s is the size of the scale used to model the component and ${\rm{\Phi }}(r)$ is a prolate spheroidal wave function (refer to Schwab (1978, 1980) for more information). In this algorithm, the components are determined by finding the global peaks in the modified residual images and by combining the corresponding scale information. Different scales are iteratively found at different locations to form complex features  to approximate any potential true emission. Here, each scale is essentially a convolution of a scaled delta function $\delta (a,l)$ and a normalized scale function ${I}^{\mathrm{norm}-\mathrm{scale}}(a,l,s)$,

$$\begin{eqnarray}&&{I}^{\mathrm{comp}}(a,l,s)={I}^{\mathrm{norm}-\mathrm{scale}}(a,l,s)\delta (a,l).\end{eqnarray} \tag{ 22 }$$

A sub-image of  at a given scale is a convolution between a series of scaled delta functions centered at different positions and this scale function. By combining sub-images at different scales, a model image can be constructed to approximate the potential true emission,

$$\begin{eqnarray}&&{I}^{\mathrm{model}}=\sum _{i=1}^{{N}_{s}}{I}_{i}^{\mathrm{norm}-\mathrm{scale}}(a,l,s){\delta }_{i}(a,l),\end{eqnarray} \tag{ 23 }$$

where N_s is the total number of components that make up the model image, which actually has only s different scales in N_s. In order to  reconstruct compact sources, s = 0 is allowed, that is, ${I}_{s=0}^{\mathrm{norm}-\mathrm{scale}}(a,l,s=0)$ is a delta function and ${I}_{s\gt 0}^{\mathrm{norm}-\mathrm{scale}}(a,l,s\gt 0)$ is a tapered, truncated parabola of scale s. If all the scales are s = 0, then the emission is decomposed into a series of delta functions just as a scale-free model. These sizes and number of scales are determined by the user, which has both flexibility and difficulty in determining the configuration.

In this algorithm, the scale found by the matching filtering is biased to smaller scales before entering the model image and updating the next residual image. A valid empirical formula achieves this bias,

$$\begin{eqnarray}&&{s}_{\mathrm{bias}}=1-0.6s/{s}_{\max },\end{eqnarray} \tag{ 24 }$$

where s_bias is the biased scale, s is the scale found by the matching filtering, and s_max is the user-specified maximum scale. This bias can achieve better reconstruction of extended sources (Cornwell 2008). For example, if there is a strong point source and a small, weak extended emission around it, the matching filtering will detect a large scale, and if no bias is used, a large scale component will be subtracted from it. Then the strong point source in the following residual image will be maintained and a large negative bowl structure around the strong point source will be produced at the same time. The elimination of the negative structure requires many components (iterations).

In scale-free CLEAN algorithms, if the loop gain for extended emission is much greater than 0.1, then component-estimation errors will result, which requires a large number of iterations to correct. Just like a scale-free CLEAN algorithm, the Cornwell–Holdaway algorithm also uses a loop gain to optimize the magnitude of each component in order to correct the errors caused by component estimation. Fortunately, the loop gain of the Cornwell–Holdaway algorithm has a wide range of choices,where a loop gain of 0.5 or more can produce a good reconstructed image.

This algorithm has been implemented in CASA as in tclean task of CASA by deconvolver = "multiscale," and it has been successfully applied in many published data sets (e.g., Rich et al. 2008).

The basic idea of  the Grisen–Moorsel–Spekkens algorithm is consistent with that of the Cornwell–Holdaway algorithm (Greisen et al. 2009); that is, an image is parameterized as multiple scale basic functions. However, scale basis functions have many choices when certain conditions are met (Cornwell 2008). The Grisen–Moorsel–Spekkens algorithm uses a two-dimensional Gaussian function to model the components instead of the tapered, truncated parabola in the Cornwell–Holdaway algorithm.

$$\begin{eqnarray}&&G\left({x}_{n},{y}_{n},{\sigma }_{n}\right)=\displaystyle \frac{1}{2\pi {\sigma }_{n}^{2}}\exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\left(x-{x}_{n}\right)}^{2}+{\left(y-{y}_{n}\right)}^{2}}{{\sigma }_{n}^{2}}\right),\end{eqnarray} \tag{ 25 }$$

where x_n and y_n are the mean values of the x and y axes respectively, and ${\sigma }_{n}$ is the standard deviation. The sky emission is approximated by a set of Gaussian functions whose widths are ${\sigma }_{n}$, which are specified by the user.

$$\begin{eqnarray}&&{I}^{\mathrm{true}}=\sum _{i=1}^{{N}_{s}}{G}_{i}\left({x}_{n},{y}_{n},{\sigma }_{n}\right){\delta }_{i}(a,l)+\epsilon .\end{eqnarray} \tag{ 26 }$$

In the Cornwell–Holdaway algorithm, pre-calculating the cross-terms between the scales specified by the user makes it possible to find the optimal component at all scales in each minor cycle. This is done without recomputing the smoothed residual image in the spatial frequency domain. However, the cross terms in the Cornwell–Holdaway algorithm  are not calculated in this algorithm and only one scale is computed during each minor cycle. Like the Cornwell–Holdaway algorithm,  the Grisen–Moorsel–Spekkens algorithm uses a scale bias term as follows to suppress the spatial scales to reduce errors,

$$\begin{eqnarray}&&{s}_{\mathrm{bias}}=1.0/{s}^{2k},\end{eqnarray} \tag{ 27 }$$

where $k\in \{0.2,0.7\}$. This is equivalent to using a Gaussian area to normalize the residual images (Rau 2010).

This algorithm is available in AIPS,⁶ and has been successfully applied to many published data sets (Thilker et al. 2002; Clarke & Ensslin 2006; Momjian et al. 2006; Kent et al. 2009).

The scale biases in the Cornwell–Holdaway and Grisen–Moorsel–Spekkens algorithms are equivalent to normalizing the smoothed residual images (Rau 2010). The Cornwell–Rau algorithm (Rau & Cornwell 2011) achieves this goal through a Hessian matrix $[{H}^{\mathrm{peak}}]$. Each element of the Hessian matrix ${H}_{s,p}^{\mathrm{peak}}$ represents the sum of the uv-tapered gridded imaging weights,

$$\begin{eqnarray}&&{H}_{s,p}^{\mathrm{peak}}={tr}\left([{T}_{s}][{S}^{-1}{WS}][{T}_{p}]\right)\,\,\forall s,p\in 0,1,\ldots ,{N}_{s}-1.\end{eqnarray} \tag{ 28 }$$

where ${tr}[M]$ is the trace of matrix $[M]$, T_s is an uv-taper function of scale s, $s,p$ denote the indices of the scale basis functions, which are tapered truncated parabolas, S is the sampling function, ${S}^{-1}$ is the inverse of S, W is the weights, and N_s is the number of scale basis functions. The normalization of the Hessian matrix is useful for the case of overlapping components at different spatial scales. This makes estimating the total flux of the component more accurate than the Cornwell–Holdaway algorithm. A measure of the area of the main lobe of the PSF at different spatial scales corresponds to the diagonal elements of this Hessian matrix. When calculating the principal solution, these normalization methods are equivalent to using the diagonal elements of the Hessian matrix. This work (Rau 2010) has carried out detailed development of the theory of CLEAN algorithms through linear algebra.  The Cornwell–Rau algorithm has been combined with multiscale CLEAN to solve the problem of wideband imaging, and has become a commonly used algorithm, such as in the tclean task of CASA by deconvoler = "mtmfs."

Like the Cornwell–Rau algorithm (Rau & Cornwell 2011),  the Offringa–Smirnov algorithm (Offringa & Smirnov 2017) made some changes to the multi-scale algorithm (Cornwell 2008),  which was extends to the wideband case using the joined-channel method, and was verified using real data. The minor cycle of this algorithm is similar to that of the Clark algorithm (Clark 1980), i.e., there is subminor cycle. In the subminor cycle, only one scale is considered and multiple iterations are performed at the current scale. This method no longer seeks the sparsest solution under the enumerated scales with an increasing number of components, which reduces the computational load. A tapered quadratic function and a Gaussian function are implemented as scale basis functions. In order to implement this modified minor cycle, a new scale-bias function was also proposed. With this automatic scale-dependent masking technique, a deeper deconvolution can be performed.

While achieving similar results, the algorithm's minor cycle is more than an order of magnitude faster than the Cornwell–Holdaway algorithm (Cornwell 2008) in the single-frequency mode, and 2–3 orders of magnitude faster than the Cornwell–Rau algorithm (Rau & Cornwell 2011) in multi-frequency mode.

There are other implementations, such as WSCLEAN (Offringa et al. 2014) and the Hybrid Matching Pursuit and sub-space optimization algorithms applied to the wideband widefield cases (Tasse et al. 2018). These multi-scale algorithms use scale basis functions to parameterize their corresponding features at multiple scales to represent the sky emission. Compared to scale-free algorithms, these algorithms can extract extended emission more efficiently. In general, multi-scale algorithms first extract large-scale features, and then gradually find smaller-scale features. This is in constrast to scale-free algorithms, which extract pixel details and scale features are broken into points.

The Wakker–Schwarz multi-resolution algorithm (Wakker & Schwarz 1988) can be classified in the category of multi-scale deconvolution. Unlike the multi-scale algorithms mentioned previously, it first extracts scale features from the lowest resolution and then extracts finer structures from higher resolution images. It considers a single scale in each minor cycle instead of processing all scales in each minor cycle, as performed by the Cornwell–Holdaway algorithm and the Cornwell–Rau algorithm. Another algorithm similar to the Wakker–Schwarz multi-resolution algorithm is the wavelet CLEAN agorithm (Horiuchi et al. 2001), which uses wavelet transforms to achieve multi-resolution reconstruction. However, these two algorithms still perform a series of Högbom minor cycles, which parameterizes the radio sky as a set of delta functions. From this perspective, they are also scale-free CLEAN algorithms.

In general, multi-scale algorithms have greatly improved the performance of extracting extended and complicated sources from radio images (Bhatnagar & Cornwell 2004). However, scale selection is limited to the chosen enumerated scales and there is no theory or method that can help guide the selection of better components. All component scales must be consistent with pre-specified scales. Features that differ from the pre-specified scale will be broken, which requires more iterations to correct. This problem is solved by finding ways to scale adaptively, as well be described in the following section.

This adaptive scale pixel decomposition algorithm (also called the Asp-Clean algorithm) (Bhatnagar & Cornwell 2004) finds the optimal component at the "energy" peak position of the current residual image, which is achieved by blurring the current residual image using several Gaussian functions and then finding the global peak. The number of model parameters monotonically increases with the number of iterations. If all the model parameters are fitted in each minor cycle, the computational complexity will be very large for a complex image containing several hundred components.  Therefore, only one active-set is optimized at a time rather than optimizing all components simultaneously,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{2}\parallel {{\boldsymbol{I}}}^{\mathrm{dirty}}-{{\boldsymbol{B}}}^{\mathrm{dirty}}{{\boldsymbol{I}}}_{i}^{\mathrm{model}}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 32 }$$

where ${{\boldsymbol{B}}}^{\mathrm{dirty}}$ and ${{\boldsymbol{I}}}_{i}^{\mathrm{model}}$ are convolutions, and only some of the parameters in ${{\boldsymbol{I}}}_{i}^{\mathrm{model}}$ can be changed, that is, the part related to the active-set. The number of components in an active-set also determines the dimensions of the search space in each iteration. The active-set of components is adaptively determined in each iteration, including  how many components are used each time.

In  the Bhatnagar–Cornwell algorithm, the length of the parameters' derivatives need to be calculated,

$$\begin{eqnarray}&&{L}_{i}=| {{\rm{\nabla }}}_{i}{\chi }^{2}| ,\end{eqnarray} \tag{ 33 }$$

where ∇_i is the derivative operation of the parameters of the ith component. The method of determining the active-set is simple. First, L_i is computed before each fitting iteration, then a threshold L₀ is used to determine if a component is in the active-set. If L_i is less than L₀, then this component will be removed from the active-set. The threshold L₀ is determined as follows,

$$\begin{eqnarray}&&{L}_{0}=\lambda \sum {I}_{i}^{\mathrm{residual}},\end{eqnarray} \tag{ 34 }$$

where λ is a hyperparameter that affects the size of the active-set, and has an experimental value. The non-orthogonality of the search space of this algorithm solves the problem of model representation in which the pixels of the extended feature are not independent.

The Bhatnagar–Cornwell algorithm finds the optimal model components by minimizing an objective function. However, it is a time-consuming process. In a typical function-optimization process, finding an optimal result often requires a large number of iterative calculations. Zhang et al. (2016a) identified the most time-consuming part of the Bhatnagar–Cornwell algorithm, and henced proposed a faster algorithm known as the Asp-Clean2016 algorithm.

In this algorithm, model components are found by minimizing a different objective function,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{2}\parallel {{\boldsymbol{I}}}_{i-1}^{\mathrm{residual}}-{{\boldsymbol{I}}}_{i}^{\mathrm{comp}}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 35 }$$

where ${{\boldsymbol{I}}}_{i}^{\mathrm{comp}}={{\boldsymbol{I}}}_{{ic}}^{\mathrm{gauss}}({\alpha }_{{ic}},{x}_{{ic}},{y}_{{ic}},{\omega }_{{ic}})$ containing the amplitude ${\alpha }_{{ic}}$, position $({x}_{{ic}},{y}_{{ic}})$ and width ${\omega }_{{ic}}$, which is a fitted Gaussian function that contains a model component whose parameters are obtained using the Levenberg–Marquardt optimization algorithm (Marquardt 1963). After each fitting convergence, the potential components can be extracted from I_ic^Gauss by an analytical method. The width of the model component contained in I_ic^Gauss is calculated as follows,

$$\begin{eqnarray}&&{\omega }_{i}=\sqrt{{\omega }_{{ic}}^{2}-{\omega }_{b}^{2}},\end{eqnarray} \tag{ 36 }$$

where ${\omega }_{i}$, ${\omega }_{{ic}}$ and ${\omega }_{b}$ are the widths of model component I_i^component, I_ic^Gaussian and the Gaussian beam, which is approximated from the PSF, respectively. After obtaining the widths of the current model components, the magnitude ${\alpha }_{i}$ of each component is estimated as follows,

$$\begin{eqnarray}&&{\alpha }_{i}=\displaystyle \frac{{\alpha }_{{ic}}{\omega }_{{ic}}^{2}}{2\pi {\alpha }_{b}{\omega }_{b}^{2}{\omega }_{i}^{2}},\end{eqnarray} \tag{ 37 }$$

where ${\alpha }_{b}$ and ${\alpha }_{{ic}}$ are the amplitudes of the Gaussian beam and I_ic^Gaussian, respectively.

Analytically deconvolving the PSF allows the convolution calculations in the objective function to be removed, which results in a significant reduction in computational effort relative to the Bhatnagar–Cornwell algorithm. The quality of the imaging obtained with both algorithms is similar.

Scale-free algorithms can represent well-separated point emission very well and their computational methods are very effective; however, such algorithms are not optimal for the representation and calculation of extended and complex emission. Although multi-scale and adaptive-scale algorithms already include a zero-scale (i.e., delta functions) for approximating compact emission, their computational methods are not very efficient. The fused-Clean algorithm (Zhang 2018a) was proposed, which combines the Högbom algorithm and the Asp-Clean2016 algorithm.

The fused-Clean algorithm automatically triggers Högbom minor cycles and Asp-Clean2016 minor cycles according to different situations. When the initial scale is less than a pre-set threshold (e.g., 1.2 times the width of the PSF) or small scales occur frequently in recent iterations (e.g., the scales of five components during the previous 10 iterations were less than 1 pixel), then the minor cycle of the Högbom algorithm is triggered. Once the minor cycle of the Högbom algorithm is triggered, it will be executed many times. The number of iterations executed each time increases with the number of triggers, which works well when the following relationship is satisfied,

$$\begin{eqnarray}&&{t}_{{tn}}=100+50\ast ({e}^{0.05\ast {tn}}-1).\end{eqnarray} \tag{ 38 }$$

where t_tn is the number of the minor cycles used by the Högbom algorithm executed at the tnth trigger.

The Asp-Clean2016 algorithm is almost always preferred when there is strong extended emission in the residual image. As the deconvolution deepens, the signal in the residual image gradually decreases, and many scale-free emissions will appear. Then the Högbom algorithm is called to reconstruct the scale-free emission. The scale-adaptive performance of the Asp-Clean2016 model effectively separates emission and noise, whereas the Högbom–Clean algorithm can effectively reconstruct dense emission. The combination of these two algorithms ensures the effective separation of the emission and noise and also accelerates the entire deconvolution process. In general, the fused-Clean algorithm is several times faster than the Asp-Clean2016 algorithm. This fusion represents an incremental improvement to the Asp-Clean2016 algorithm and is a useful step toward making the Asp-Clean algorithm more practically usable.

There are also other scale-adaptive algorithms, such as the Algas-Clean algorithm (Zhang et al. 2016b; Zhang 2018b). This algorithm not only has the ability to perform scale adaptation, but it also achieves accurate reconstruction of the signal through a shape-dependent adaptive loop gain scheme.

CS theory (Candes 2006; Candes & Wakin 2008) tells us that when a signal is sparse, or can be sparsely represented under a given basis function dictionary, the signal can be reconstructed with much less sampling than required by the Nyquist–Shannon theory. Sky brightness is sparse in the image domain for compact sources, and it is also sparse in a basis function dictionary when the sky brightness contains diffuse sources. Most of these methods use L1 norms, which can achieve a more sparse solution, with a certain regularization. Because the Fourier domain and the image domain are perfectly incoherent (Candes et al. 2006), CS can be applied directly to the image domain to reconstruct compact sources. For example, the Partial Fourier based CS deconvolution method (Li et al. 2011) applies the L1 norm in the image domain,

$$\begin{eqnarray}&&\min \parallel {\boldsymbol{I}}{\parallel }_{l1}\qquad s.t.\parallel {MFI}-V{\parallel }_{l2}\leqslant \epsilon ,\end{eqnarray} \tag{ 39 }$$

where M is the measurement matrix. This method can effectively reconstruct compact sources. Wavelet transformations are often used to achieve sparse transforms of diffuse sources in radio astronomy. Many studies are based on wavelet transforms such as the undecimated wavelet transform (Starck et al. 2007), which use different strategies or objective functions to apply CS during the deconvolution of diffuse radio images, such as Li et al. (2011), SARA (Carrillo et al. 2012), PURIFY (Carrillo et al. 2014), (Garsden et al. 2015), MORESANE (Dabbech et al. 2015) and SASIR (Girard et al. 2015).

The emergence of CS-based deconvolution (e.g., Pratley et al. 2018) has provided a new alternative image reconstruction scheme for radio synthesis imaging which is basically different to the schemes used by CLEAN algorithms (Wiaux et al. 2009a). The CLEAN algorithms yield high-quality reconstructed images by introducing scale priors and combining L2 norms, while CS methods use sparse priors that are combined with L1 norms and a regularization to obtain high-quality images for different cases. CS methods can reconstruct rich structures that contain sharper and more details in specific and simple scenarios. However, most of them have been validated only for relatively simple test samples, and they are deemed unstable for application to real data (Offringa & Smirnov 2017).

MEMs (Cornwell & Evans 1985; Narayan & Nityananda 1986; Cornwell et al. 1999; Akiyama et al. 2017b; EHT Collaboration et al. 2019) are technically similar to CS methods in that they often introduce a priori information through regularization. This regularization is implemented by an entropy function, which is a function of the model emission reconstructed from the dirty image. During the reconstruction of an image, the entropy function is maximized by making the model's visibilities fit mostly to the measured visibilities. Many studies have discussed the construction of entropy functions from theoretical and philosophical perspective, and a large number of these algorithms are based on Bayesian statistics. Some common entropy functions $E({I}^{\mathrm{model}})$ have been described by Thompson et al. (2017), such as

$$\begin{eqnarray}&&{E}_{1}=-\sum _{i}\mathrm{log}{I}_{i}^{\mathrm{model}},\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{E}_{2}=-\sum _{i}\displaystyle \frac{{I}_{i}^{\mathrm{model}}}{{\sum }_{i}{I}_{i}^{\mathrm{model}}}\mathrm{log}\displaystyle \frac{{I}_{i}^{\mathrm{model}}}{{\sum }_{i}{I}_{i}^{\mathrm{model}}}.\end{eqnarray} \tag{ 41 }$$

MEM still needs to minimize the errors resulting from the measured visibilities and model visibilities. The difference relative to CS methods is that MEM uses an entropy function instead of L1 norms.

Both MEM and CLEAN have been used to estimate the true sky emission from incomplete noisy measurements. However, CLEAN is more straightforward and easier to implement. There are some differences between these two types of algorithms: MEM employs explicit optimization, which is not the approach used by most CLEAN algorithms. Moreover, MEM requires an initial model, but CLEAN algorithms do not. MEM images tend to be smooth while CLEAN images show more details. Additionally, MEM is not good at point-source reconstruction, but multi-scale and adaptive-scale CLEAN algorithms are capable of processing point sources and extended sources simultaneously.