© The Authors 2025

1 Introduction

Currently, radio interferometry has become an indispensable tool for astronomical observation, and it plays an important role in many fields of astronomy and astrophysics. It is well known that only part of the complete uv coverage can be measured by the interferometric array. The imaging process is a typical, ill-posed, linear inverse problem. Astronomers routinely employ the Hogbom CLEAN algorithm (Högbom 1974) along with its variants for radio interferometric imaging, such as multi-scale CLEAN (Cornwell 2008, 2009; Rich et al. 2008) and adaptive-scale CLEAN (Bhatnagar & Cornwell 2004; Zhang et al. 2016; Zhang 2018; Zhang et al. 2020). Although the CLEAN family effectively reduces side-lobe interference caused by incomplete uv coverage, the CLEAN family often produces sub-optimal imaging quality.

In recent years, the theory of compressed sensing (CS; Donoho 2006; Li et al. 2020) has provided a new solution for radio interferometric imaging. By utilizing the sparsity of the signal, CS can reconstruct high-quality images at low sampling rates. However, the L₀ -norm optimisation problem in CS is non-deterministic polynomial time hard (NP-hard). Aiming to resolve this problem, researchers have proposed a variety of alternative methods to obtain sub-optimal solutions. Among them, the widely used method is convex optimisation for solving the L₁-norm minimisation problem or reweighted L₁-norm minimisation problem (Condat 2014; Bubeck et al. 2015; Xu & Noo 2022).

Based on convex optimisation techniques, many CS-based approaches have been developed in radio interferometric imaging, such as the isotropic undecimated wavelet transform (IUWT; Li et al. 2011; Garsden et al. 2015) and sparsity-averaging reweighted analysis (SARA; Carrillo et al. 2012, Carrillo et al. 2014; Onose et al. 2016, 2017; Dabbech et al. 2022; Terris et al. 2023; Wilber et al. 2023a). The IUWT algorithm adopts an isotropic undecimated wavelet basis as the sparse basis according to the characteristics of astronomical images and then solves the L1 -norm minimisation problem by the fast iterative shrinkagethresholding algorithm. Experimental results show that its reconstruction quality is superior to the Högbom-CLEAN and multiscale CLEAN algorithms. The SARA algorithm uses a dictionary composed of a Dirac basis and multiple wavelet bases as the sparse basis, and then solves the reweighted L₁ -norm minimisation problem by a variety of fast and parallelisable algorithms, such as the primal-dual forward-backward (PDFB) one, the alternating direction method of multipliers (ADMM), and the preconditioned primal-dual (PPD; Onose et al. 2016, 2017) method. The simulation results have shown that the reconstruction quality of the SARA is superior to that of the total variation (TV) and IUWT algorithms (Carrillo et al. 2012). In view of the excellent performance of the SARA algorithm, researchers have extended it to the polarized imaging and wide-band imaging in radio interferometry, namely the Polarized SARA (Birdi et al. 2018, 2020) and HyperSARA algorithms (Abdulaziz et al. 2019; Thouvenin et al. 2023a,b).

In addition, with the rapid development of artificial-intelligence technology, researchers have begun to explore the use of deep neural networks (DNNs) for end-to-end reconstruction in radio interferometry (Gheller & Vazza 2022). Furthermore, an artificial intelligence for regularisation in radio-interferometric Imaging (AIRI) approach is presented (Wilber et al. 2023b; Terris et al. 2023), which simultaneously inherits the advantages of convex optimisation and deep learning. Recently, astronomers have proposed a residual-to-residual DNN series for a high-dynamic-range imaging (R2D2) algorithm (Dabbech et al. 2024; Aghabiglou et al. 2024), which adopts a DNN series architecture based on a residual-to-residual method. It uses a DNN to directly learn to generate updates from the current image estimate and the residual dirty image.

It is worth noting that traditional compressed sensing methods approximately replace the L₀ norm by the L₁ norm or weighted L₁ norm, resulting in a bias problem. Specifically, the L₀ norm only cares about whether an element is non-zero rather than its actual magnitude. However, the L₁ norm punishes the elements based on their size. During the reconstruction process, small non-zero elements are easily compressed to zero, and the estimated values of large non-zero elements are also smaller than their actual values, resulting in a deviation between the reconstructed signal and the true sparse signal. Essentially, the bias stems from the compromise made in the accuracy of the sparsity measurement in order to achieve the solvability of the optimisation problem owing to the fact that convex optimisation has a global optimal solution (Zou & Hastie 2005; Candes et al. 2008; Jia et al. 2018).

To ameliorate the bias problem, Yang et al. (2025) proposed a L_q proximal gradient (LPG) algorithm with 0 < q < 1. However, the L_q norm (Marjanovic & Solo 2012) is still a biased approximation of the L₀ norm, and the derivative of the L_q norm at the zero point does not exist. In contrast, the smoothly clipped absolute deviation (SCAD; Fan & Li 2001; Wen et al. 2018) has asymptotic unbiasedness under certain conditions. Moreover, it is continuously differentiable and can efficiently find the optimal solution in the optimisation process. To efficiently obtain an accurate solution in radio interferometry, we propose a non-convex, sparse regularisation method dubbed improved SCAD (ISCAD) in this paper. Our proposed approach employs the continuous SCAD penalty to approximate the L0 norm and leverages an improved proximal-gradient algorithm to effectively address the non-convex optimisation problem. Compared with the existing proximal gradient algorithms (Wen et al. 2018), the main improvements of the proposed algorithm lie in the following two aspects. On one hand, we adopted an adaptive, non-monotonic step-size (ANMS) strategy (Liu et al. 2022) to enhance the convergence rate of the proximal gradient algorithm for a non-convex optimisation problem. On the other hand, we adaptively updated the regularisation parameter using the prior information of the image.

The paper is structured as follows. Section 2 briefly introduces typical radio interferometric imaging methods based on CS. Section 3 describes a non-convex sparse regularisation method based on SCAD in detail. In Sect. 4, we describe several simulation experiments carried out to assess the performance of the proposed method. Section 5 provides the final conclusion.

2 Radio interferometric imaging methods based on compressed sensing

A radio interferometer generates the visibilities by measuring the correlation between the signals received by a pair of spatially separated antennas. Based on the Van Cittert–Zernike theorem (Thompson et al. 2017), under the simplifying assumptions of a coplanar array the relationship between the signal, x, and the measured visibilities, y, is given by

$$y(u)={\displaystyle \int x}(l)A(l)e^{-i2\pi l·u}dl,$$(1)

where u = (u, v) represents the baseline vector normalised by the wavelength, l = (l, m) is the coordinate on the image plane, and A(l) is the primary beam that limits the observed field of view. After the discretisation of Eq. (1), a linear measurement model can be obtained as follows:

$$y\text{ = }\Phi x\text{ + }n\text{,}$$(2)

where the vectors y ∈ C^M and x ∈ R^N stand for the visibilities and sampled intensity signal separately, n ∈ C^M represents the noise vector, and the matrix Φ = GFZ ∈ C^M×N represents the linear mapping from the image domain to the visibility domain. The matrix G ∈ C^M×nN includes compact support interpolation kernels and models the direction-dependent effects (DDEs), F ∈ C^nN×nN represents an n-oversampled discrete Fourier operator, and Z ∈ R^nN×N denotes an oversampling and scaling operator for pre-compensated interpolation. It is unfortunate that Eq. (2) is an ill-posed inverse problem (Rau et al. 2009). At present, the CS technology has been successfully applied to radio interferometric imaging, and the state-of-the-art CS-based methods in radio interferometry are the SARA and AIRI.

According to CS theory, the signal, x, needs to be sparse or compressible on the orthogonal basis, ψ, namely x = ψα, where α denotes the decomposition of the signal on the sparse basis. While the measurement matrix satisfies the restricted isometry property (RIP; Candes et al. 2011; Shen & Hu 2019), the signal can be recovered by solving the following minimisation problem:

$$\underset{\alpha }{\text{min}}{\Vert \alpha \Vert }_0\text{subject to\hspace{0.17em}}{\Vert y-\Phi \Psi \alpha \Vert }_2\le \epsilon \text{, }$$(3)

where ||·||₀ represents the L₀ norm of the vector, α, and ε is the upper bound of the L₂ norm of the noise vector. It should be noted that Eq. (3) is NP-hard and intractable.

The SARA algorithm proposed by Carrillo et al. (2012) adopts the concatenation of a Dirac basis and eight Daubechies wavelet bases (Db1-Db8) as a sparsity dictionary in radio interferometric imaging. With respect to a single sparse basis, the average of multiple sparse bases can promote the sparsity of the signal. And the SARA algorithm replaced Eq. (3) by the reweighted L₁ minimisation problem; this is defined as

$$\underset{x}{\text{min}}{\Vert W{\Psi }^{†}x\Vert }_1\text{subject to\hspace{0.17em}}{\Vert y-\Phi x\Vert }_2\le \epsilon \text{, }$$(4)

where $x\in R_{+}^N$ denotes the reality and positivity prior on x, and W stands for the diagonal matrix with positive weights. In this work, PPD (Onose et al. 2017) was employed to solve Eq. (4).

Based on the unconstrained version of the SARA algorithm, namely uSARA, the AIRI algorithm (Terris et al. 2023) trains a DNN as the denoiser to replace the proximal operator in the regularisation step of the uSARA algorithm:

$${\tilde{x}}_{(k+1)}=D\left({\tilde{x}}_{(k)}-\gamma \nabla f\left({\tilde{x}}_{(k)}\right)\right),$$(5)

where D refers to the DNN denoiser, γ is the step size, and ∇f is the gradient of the data-fidelity term $f\left(x\right)=\frac{1}{2}{\Vert y-\Phi x\Vert }_2^2$. Experimental results show that the AIRI algorithm outperforms the traditional SARA algorithm in terms of both reconstruction accuracy and speed (Wilber et al. 2023b).

We note that traditional CS algorithms generally approximate the NP-hard L₀-norm minimisation problem by the L₁-norm or reweighted one, which introduces certain deviations and affects the accuracy of the reconstruction results. In order to remedy the bias issue, Yang et al. (2025) proposed a radio-interferometric-imaging method based on L_q norm. The leastsquares problem with the L_q norm is given by

$$\underset{x\in {ℝ}_{+}^N}{\min }\frac{1}{2}{\Vert y-\Phi x\Vert }_2^2+\lambda {\Vert {\Psi }^{†}x\Vert }_q^q,$$(6)

where λ refers to the regularisation parameter, ${\Vert \cdot \Vert }_q^q$ stands for L_q norm, and the dictionary used in LPG is the same as that in SARA. Equation (6) is efficiently solved by the proximal gradient algorithm.

3 Non-convex, sparse regularisation based on smoothly clipped absolute deviation

Theleast-squares problem with the SCAD term can be expressed as

$$\underset{x\in {ℝ}_{+}^N}{\min }\frac{1}{2}{\Vert y-\Phi x\Vert }_2^2+P_{\lambda }\left({\Psi }^{†}x\right),$$(7)

where the dictionary ψ used is the same as that in SARA and LPG. P_λ(·) represents the SCAD penalty, which is defined as follows (Fan & Li 2001):

$$P_{\lambda }(z)=\{\begin{array}{ll}\lambda \left|z\right|\hfill & \left|z\right|<\lambda \hfill \\ \frac{2a\lambda \left|z\right|-z^2-{\lambda }^2}{2(a-1)}\hfill & \lambda \le \left|z\right|<a\lambda ,\hfill \\ (a+1){\lambda }^2/2\hfill & \left|z\right|\ge a\lambda \hfill \end{array}$$(8)

where a > 2 is a constant and usually set to 3.7 (Fan & Li 2001). From Eq. (8), it can be seen that when the coefficient, z, is small, SCAD behaves similarly to the L₁ norm, imposing strong penalties on small coefficients to make them zero. As the coefficient increases, SCAD takes on a smooth quadratic function form, avoiding excessive penalties for large coefficients. Figure 1 compares different penalty functions. It should be remarked that the L_q norm bridges the gap between the L₀ norm and the L₁ norm. Moreover, SCAD exhibits a closer approximation to the L₀ norm than the L_q norm.

The proximal operator is of vital significance in developing efficient and scalable algorithms for numerous optimisation problems, especially for non-convex and non-smooth optimisation problems. For a penalty function, P_λ(·), the corresponding proximal mapping operator is given by

$${\text{prox}}_{\lambda }\left(r\right)\arg \underset{x}{\min }{P}_{\lambda }(z)+\frac{1}{2}{\Vert z-r\Vert }^2.$$(9)

For the SCAD penalty, the corresponding proximal operator is given by

$${\text{prox}}_{\lambda }\left(t\right)=\{\begin{array}{ll}sign(t)\max \left\{\right|t|-\lambda ,0\}\hfill & \left|t\right|\le 2\lambda \hfill \\ \frac{(a-1)t-a\lambda sign(t)}{(a-2)}\hfill & 2\lambda <\left|t\right|\le a\lambda ,\hfill \\ t\hfill & \left|t\right|>2\lambda \hfill \end{array}$$(10)

where sign(·) represents the sign function. For the sake of facilitating comparison, the penalty functions and the corresponding proximal operators of L₀ norm, L₁ norm, and L_q norm are listed in Table 1.

Figure 2 shows the outputs of proximal operators corresponding to different penalty functions. One can notice that the proximal operator of L₀ norm is discontinuous, and the proximal operator of L₁ norm is continuous but has deviations from that of L₀ norm. The proximal operator of L_q norm is discontinuous and is less biased than that of the L₁ norm. In contrast, the proximal operator of SCAD not only remains continuous, it also results in unbiased estimates of large coefficients.

To effectively address the non-convex optimisation problem, we adopt an improved proximal-gradient algorithm (Wen et al. 2018). The reconstructed image is given by

$${\tilde{x}}_{(k+1)}=\Psi {\text{prox}}_{\lambda }\left({\Psi }^{†}\left(s_{(k)}+{\alpha }_{k+1}{\Phi }^{†}\left(y-\Phi s_{(k)}\right)\right)\right),$$(11)

where prox_λ(·) stands for the proximal operator defined in Eq. (10), and the step size α_κ+1 ∈ (0, 2/L) is generally a constant, where L is the Lipschitz constant of ∇f. s_(k) denotes the temporary vector and is given by

$$s_{(k)}={\tilde{x}}_{(k)}+\frac{t_{(k)}-1}{t_{(k+1)}}\left({\tilde{x}}_{(k)}-{\tilde{x}}_{(k-1)}\right),$$(12)

where $t_{\left(k+1\right)}=\left(1+\sqrt{1+4\left(t_{\left(k\right)}{}^2\right)}\right)/2$ denotes the acceleration factor, and its initial value t(0) is usually set to one.

In addition, we adopted two strategies to accelerate the convergence speed of the algorithm. The first strategy is a gradientrestart one. It means that if ${\left(s_{\left(k+1\right)}-{\tilde{x}}_{\left(k+1\right)}\right)}^{\text{T}}\left({\tilde{x}}_{\left(k+1\right)}-{\tilde{x}}_{\left(k\right)}\right)>0$, the acceleration factor t(k+1) is reset to one, thereby eliminating oscillations (O’donoghue & Candes 2015).

Another strategy is an ANMS one. Liu et al. (2022) designed an ANMS strategy to increase the convergence rate of fast iterative shrinkage-thresholding algorithm (FISTA) for non-smooth and convex minimisation problems. In this paper, we applied the ANMS strategy to enhance the convergence rate of the proximal gradient algorithm for the non-convex optimisation problem. In the ANMS strategy, the acceleration factor is rewritten as

$$t_{(k+1)}=(1+\sqrt{1+4{\theta }_{(k)}(t_{(k)})})/2,$$(13)

where θ_(k) = α_(k)/α_(k+1); the step size, α_(k+1), is adaptively updated by

$$\text{If\hspace{0.17em}}\langle \nabla f\left({\tilde{x}}_{(k+1)}\right)-\nabla f(s_{(k)}),{\tilde{x}}_{(k+1)}-s_{(k)}\rangle >\frac{{\mu }_0}{{\alpha }_{(k)}}{\Vert {\tilde{x}}_{(k+1)}-s_{(k)}\Vert }_2^2$$

$${\alpha }_{(k+1)}={\mu }_1\frac{{\Vert {\tilde{x}}_{(k+1)}-s_{(k)}\Vert }_2^2}{\langle \nabla f\left({\tilde{x}}_{(k+1)}\right)-\nabla f\left(s_{(k)}\right),{\tilde{x}}_{(k+1)}-s_{(k)}\rangle },$$

else

$${\alpha }_{(k+1)}=\left(E_{(k)}+1\right){\alpha }_{(k)}$$(14)

where 〈·〉 represents the inner product, µ₀, µ₁ are two constants satisfying 0 <µ₁ < µ₀ < 1, and $E_{\left(k\right)=\raisebox{1ex}{$w$}\!\left/ \!\raisebox{-1ex}{$k^{1.1}$}\right.}$ is the control parameter with a constant w, which can be adjusted adaptively to control the growth amplitude of the step size.

The selection of the regularisation parameter, λ, is of vital importance to the performance of the algorithm. Generally, it is set based on experience and keeps constant during the iterative process. However, the empirical selection method poses challenges to achieving an optimal equilibrium between the data fidelity term and the regularisation term. To address this issue, this paper proposes an adaptive selection method, which makes full use of the information of the data fidelity term and the regularisation term to update the regularisation parameter (Yang et al. 2024). The regularisation parameter is updated by

$${\lambda }_{(k+1)}=\ln \left(\varsigma \frac{{\Vert y-\Phi {\tilde{x}}_{(k)}\Vert }_2^2}{{\Vert {\Psi }^{†}{\tilde{x}}_{(k)}\Vert }_1}+1\right),$$(15)

where the factor ζ is a constant.

Finally, the computational complexity of different algorithms was analysed. For the SARA algorithm solved by PPD, the computational complexity is O(nN log nN) + O(dN) + O(MN) per iteration, where n denotes the oversampling factor, and d denotes the number of blocks into which the visibility data is divided (Onose et al. 2017). For the LPG and ISCAD algorithms, the computational complexity is the same and is O(nN log nN) + O(N) per iteration. It can be seen that the computational complexity of the LPG and ISCAD algorithms per iteration is evidently lower than that of the SARA algorithm.

It should be noted that, with respective to convex optimisation, the specificities of the proposed algorithm mainly lie in the differences in the approximation accuracy of sparsity and the optimisation characteristics. Compared with convex optimisation, the proposed method utilises piecewise-designed SCAD penalty to be closer to the L0 norm, thereby reducing the bias problem. However, non-convex optimisation faces challenges such as local optima and low convergence efficiency. Therefore, the proposed method adopts restart and ANMS strategies to address these issues. From the perspective of application scope, the proposed algorithm is more suitable for handling high sparsity or structured sparse signals.

Fig. 1 Comparison of different penalty functions.

Fig. 2 Outputs of proximal operators corresponding to penalty functions.

Fig. 3 (a) VLA baseline coverage. (b) Test image from 30Dor, presented on a log10 scale.

4 Simulation results

We conducted numerical simulation experiments based on the Very Large Array (VLA; Perley et al. 2011) and Square Kilometre Array (SKA; Dewdney et al. 2009). We compared the performance of the proposed ISCAD method with respect to classical imaging methods such as the SARA¹, AIRI², and LPG³ methods. All programs were executed in the same hardware and software environment: MATLAB R2022a running on a personal computer equipped with the 2.5GHz Intel i9-12900H processor and 16 GB memory. In addition, we optimised and set the parameters of SARA, AIRI, and LPG according to Onose et al. (2017), Terris et al. (2023), and Yang et al. (2025). The DNN denoiser of AIRI in this article is the same as that trained on optical astronomy images in the Terris et al. (2023). In the LPG method, the value of q in the L_q norm is 0.8 (Yang et al. 2025). The main parameters of ISCAD are set as follows: α = 3.7, µ₀ = 0.99, µ₁ = 0.95, w = 1, ζ ∈ [10^–6,10^–3].

The signal-to-noise ratio (S/N) and S/N_log were used to measure the reconstruction accuracy. They can be expressed as follows:

$$\text{S/N}=20{\log }_{10}\left(\frac{{\Vert x\Vert }_2}{{\Vert x-\tilde{x}\Vert }_2}\right),$$(16)

$${\text{S/N}}_{\log }=20{\log }_{10}\left(\frac{{\Vert {\log }_{10}\left(\frac{x}{\epsilon }+{\text{I}}_N\right)\Vert }_2}{{\Vert {\log }_{10}\left(\frac{x}{\epsilon }+{\text{I}}_N\right)-{\log }_{10}\left(\frac{\tilde{x}}{\epsilon }+{\text{I}}_N\right)\Vert }_2}\right),$$(17)

where x and $\tilde{x}$ denote the original image and the reconstructed image separately, and ε refers to the estimated noise level, and IN is the unit vector with the length N.

In the simulations, we assumed a small field of view and the absence of DDEs such that the measurement operator reduces to a Fourier operator. We employed an oversampled Fourier transform, F, with n = 4 and a matrix, G, that interpolates frequency data to average a nearby uniformly distributed frequency via 8 × 8 Kaiser-Bessel interpolation kernels (Fessler & Sutton 2003).

In the initial experiment, We used the VLA baseline coverage generated by the CASA and CASACORE (Onose et al. 2016, 2017). Figure 3a depicts the VLA baseline coverage containing 894240 uv points. Meanwhile, Fig. 3b illustrates the test image with a pixel size of 256 × 256, presented on a log10 scale. This test image was derived from the 30 Doradus (30Dor)⁴ in the Large Magellanic Cloud, and has brightness temperature values spanning the interval [0.01, 1].

When the input S/N (iS/N) of the additive Gaussian noise is set at 35dB, Fig. 4 illustrates the reconstruction results obtained by the SARA, AIRI, LPG, and ISCAD algorithms, respectively. In Fig. 4, the left column exhibits the reconstructed images on a log10 scale. the middle column displays the error images $\Delta x=x-\tilde{x}$ on a linear scale, while the right column represents the real part of the residual images also in a linear scale. An inspection of Fig. 4 reveals that the inversion result of SARA manifests large residual errors and distinct artefacts for both the cobweb-like structure and the background. In comparison to SARA, AIRI, and LPG yield more favourable inversion results, with a notable reduction in both errors and artifacts. Furthermore, ISCAD produces a reconstructed image characterised by scarcely perceptible errors and artefacts in the cobweb-like structure and the background.

In order to quantitatively assess the performance of diverse approaches, we computed the S/N, S/N_log, and running time. The relevant data are presented in Table 2. The results demonstrate that, with respect to S/N and S/N_log, ISCAD surpasses the SARA, AIRI, and LPG approaches. Besides this, the running time of ISCAD is shorter compared to that of the other three approaches.

We also quantitatively assessed the influence of noise interference on the reconstruction results. Analogously to Fig. 4, the measured visibilities degraded by Gaussian noise of varying levels were utilised for image reconstruction. Figure 5 shows the S/N and S/N_log performance of the algorithms under different noise levels. Each value presented in Fig. 5 is derived from averaging 100 simulation results. As can be seen from Fig. 5, even though the noise level changes, ISCAD surpasses the SARA, AIRI, and LPG methods in terms of reconstruction quality. This validates that ISCAD has stronger robustness against noise interference compared to the SARA, AIRI, and LPG methods.

To further verify the performance of the ISCAD approach, two other experiments were carried out using the simulated SKA Array. The SKA baseline coverage generated by CASA and CASACORE is presented in Figs. 6a and 7a; these consist of 1 447 950 and 578 358 points, respectively (Onose et al. 2016, 2017). The test images were sourced from the W28 supernova remnant⁵ and Messier 106 (Shimwell et al. 2022). The two test images displayed on a log10 scale are shown in Figs. 6b and 7b, respectively.

When the Gaussian noise with iS/N=35dB is added, the reconstruction results for the W28 and Messier 106, obtained through the SARA, AIRI, LPG, and ISCAD algorithms, are presented in Figs. 8 and 9. It can be noted that the image reconstructed by SARA have large residual errors, especially for the bright source region. In contrast, AIRI results in a visually better image with a smaller residual error. Additionally, LPG and ISCAD produce comparable results, which are marginally superior to that of AIRI in the bright source region.

For the reconstruction result in Figs. 8 and 9, the corresponding quantitative indicators such as the S/N, S/N_log , and computational time are separately presented in Tables 3 and 4. It can be noted that in terms of reconstruction quality, ISCAD performs slightly better than AIRI and LPG, and distinctly better than SARA. Meanwhile, the calculation speed of ISCAD is evidently faster than that of the SARA, AIRI, and LPG algorithms. Furthermore, we conducted a detailed analysis of computation time for LPG and ISCAD. In Table 3, the time of the singlemodel update step for LPG and ISCAD is 3.1863 s and 3.0041 s, respectively, and the number of iterations for LPG and ISCAD is 100 and 80, respectively. In Table 4, the time of the single model update step for LPG and ISCAD is 1.7418 s and 1.6472 s, respectively, and the number of iterations for LPG and ISCAD is 90 and 65, respectively. Therefore, under the same computing resources, since ISCAD converges slightly faster than LPG, the calculation time of ISCAD is slightly shorter than that of LPG.

Furthermore, the convergence performance of ISCAD is analyzed. Given that Eq. (7) is non-convex, it is difficult to offer a theoretical proof of its global convergence. In this context, we present empirical evidence to demonstrate the favourable convergence performance of ISCAD. Figure 10 indicates the convergence performance of ISCAD under different noise levels, where the objective function is $\frac{1}{2}{\Vert y-\Phi \tilde{x}\Vert }_2^2+P_{\lambda }\left({\Psi }^{†}\tilde{x}\right)$. From Fig. 10, we can see that the objective function of ISCAD can rapidly decrease to a stable value for both W28 and Messier 106, demonstrating the excellent convergence performance of ISCAD.

Fig. 4 Reconstruction results of 30Dor. (a, d, g, j) Results on a log10 scale reconstructed by SARA, AIRI, LPG, and ISCAD, respectively. (b, e, h, k) Error images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively. (c, f, i, l) Real part of the residual images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively.

Fig. 5 (a) S/N and (b) S/Nlog performance of algorithms under different noise levels.

Fig. 6 (a) SKA baseline coverage. (b) Test image with a pixel size of 1024 × 1024 from the W28 supernova remnant, shown on a log10 scale.

Fig. 7 (a) SKA baseline coverage. (b) The test image with a pixel size of 512 × 512 from the Messier 106, shown on a log10 scale.

Fig. 8 Reconstruction results of W28. (a, d, g, j) Results on a log10 scale reconstructed by SARA, AIRI, LPG, and ISCAD, respectively. (b, e, h, k) Error images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively. (c, f, i, l) Real part of the residual images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively.

Fig. 9 Reconstruction results of Messier 106. (a, d, g, j) Results on a log10 scale reconstructed by SARA, AIRI, LPG, and ISCAD, respectively. (b, e, h, k) Error images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively. (c, f, i, l) Real part of the residual images in linear scale for SARA, AIRI, LPG, and ISCAD, respectively.

Fig. 10 (a, b) Convergence performance of ISCAD algorithm for W28 and Messier 106, respectively.

5 Conclusions

The CS methods based on convex optimisation techniques such as SARA and AIRI have already demonstrated outstanding performance in radio interferometric imaging. However, convex optimisation approximately replaces the L₀ norm by the L₁ norm or reweighted L₁ norm, which gives rise to the bias problem.

Different from convex optimisation, the proposed method adopts the continuous SCAD penalty function to approximate the L₀ norm and efficiently solves the non-convex optimisation problem by using an improved proximal-gradient algorithm. The improved proximal-gradient algorithm introduces a restart strategy and an adaptive non-monotonic step-size strategy, which greatly accelerates the convergence speed of the algorithm. At the same time, the prior information of the image is used to adaptively update the regularisation parameter in the iterative process. Based on the VLA and SKA arrays, we conducted numerical simulation experiments and compared the proposed ISCAD method with state-of-the-art imaging methods such as SARA, AIRI, and LPG. The results show that ISCAD performs better than other methods in terms of reconstruction quality and computational speed.

Data availability

The code underlying this article is publicly available in https://github.com/yangxiaoch209/Radio-interferometric-Imaging.

Acknowledgements

This work is supported by the Specialized Research Fund for State Key Laboratory of Solar Activity and Space Weather. This work has used the data analysis platform of the DAocheng Radio Telescope in Chinese Meridian Project Phase II.

References

All Tables

All Figures

	Fig. 1 Comparison of different penalty functions.
In the text

	Fig. 2 Outputs of proximal operators corresponding to penalty functions.
In the text

	Fig. 3 (a) VLA baseline coverage. (b) Test image from 30Dor, presented on a log10 scale.
In the text

	Fig. 4 Reconstruction results of 30Dor. (a, d, g, j) Results on a log10 scale reconstructed by SARA, AIRI, LPG, and ISCAD, respectively. (b, e, h, k) Error images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively. (c, f, i, l) Real part of the residual images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively.
In the text

	Fig. 5 (a) S/N and (b) S/Nlog performance of algorithms under different noise levels.
In the text

	Fig. 6 (a) SKA baseline coverage. (b) Test image with a pixel size of 1024 × 1024 from the W28 supernova remnant, shown on a log10 scale.
In the text

	Fig. 7 (a) SKA baseline coverage. (b) The test image with a pixel size of 512 × 512 from the Messier 106, shown on a log10 scale.
In the text

	Fig. 8 Reconstruction results of W28. (a, d, g, j) Results on a log10 scale reconstructed by SARA, AIRI, LPG, and ISCAD, respectively. (b, e, h, k) Error images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively. (c, f, i, l) Real part of the residual images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively.
In the text

	Fig. 9 Reconstruction results of Messier 106. (a, d, g, j) Results on a log10 scale reconstructed by SARA, AIRI, LPG, and ISCAD, respectively. (b, e, h, k) Error images on a linear scale for SARA, AIRI, LPG, and ISCAD, respectively. (c, f, i, l) Real part of the residual images in linear scale for SARA, AIRI, LPG, and ISCAD, respectively.
In the text

	Fig. 10 (a, b) Convergence performance of ISCAD algorithm for W28 and Messier 106, respectively.
In the text