A Nonlinear Acceleration Method for Iterative Algorithms

TL;DR

This paper introduces a nonlinear acceleration technique that enhances the convergence speed and stability of various iterative algorithms used in signal processing and sparse recovery tasks.

Contribution

The paper proposes a novel nonlinear acceleration method that improves convergence and stability of multiple iterative algorithms, including IRLS, IMAT, SL0, ADMM, and Chebyshev Acceleration.

Findings

Enhances convergence speed of iterative algorithms.

Stabilizes and extends the stability range of existing methods.

Applicable to a variety of algorithms in signal processing.

Abstract

Iterative methods have led to better understanding and solving problems such as missing sampling, deconvolution, inverse systems, impulsive and Salt and Pepper noise removal problems. However, the challenges such as the speed of convergence and or the accuracy of the answer still remain. In order to improve the existing iterative algorithms, a non-linear method is discussed in this paper. The mentioned method is analyzed from different aspects, including its convergence and its ability to accelerate recursive algorithms. We show that this method is capable of improving Iterative Method (IM) as a non-uniform sampling reconstruction algorithm and some iterative sparse recovery algorithms such as Iterative Reweighted Least Squares (IRLS), Iterative Method with Adaptive Thresholding (IMAT), Smoothed l0 (SL0) and Alternating Direction Method of Multipliers (ADMM) for solving LASSO problems family (including Lasso itself, Lasso-LSQR and group-Lasso). It is also capable of both accelerating and stabilizing the well-known Chebyshev Acceleration (CA) method. Furthermore, the proposed algorithm can extend the stability range by reducing the sensitivity of iterative algorithms to the changes of adaptation rate.

Equations

$$E\triangleq I-G$$

$$G^{-1}$$

$$\Rightarrow G_{{}_{<K>}}^{-1}=I+\sum _{n=1}^K{E}^n$$

$$G_{{}_{<K+1>}}^{-1}=I+E(G_{{}_{<K>}}^{-1}).$$

$$G^{-1}$$

$$\Rightarrow G_{{}_{<K,\lambda >}}^{-1}=I+\sum _{n=1}^K{E_{\lambda }}^n$$

$$x_k=\lambda (x_0-G(x_{k-1}))+x_{k-1}$$

$$P_0:\underset{s}{\min }{\Vert s\Vert }_0;As=b,P_1:\underset{s}{\min }{\Vert s\Vert }_1;As=b$$

$$f_{\sigma }(s)=N-\sum _{n=1}^N{F}_{\sigma }(s[n]);F_{\sigma }(s[n])=\exp (-\frac{{|s[n]|}^2}{2{\sigma }^2})$$

$$x_k$$

$$T(x)$$

$$x_1=\frac{2}{A+B}{x}_0,{\lambda }_k=(1-\frac{{\rho }^2}{4}{\lambda }_{k-1}{)}^{-1};\text{for}k>1:$$

$$x_k=(x_1+x_{k-1}-\frac{2}{A+B}G(x_{k-1})-x_{k-2}){\lambda }_k+x_{k-2}$$

$$e_k[n]={\alpha }_{k-1}[n]\times e_{k-1}[n];e_{k-1}[n]\ne 0$$

$${\hat{x}}_3-x=\alpha ({\hat{x}}_2-x)={\alpha }^2({\hat{x}}_1-x).$$

$$x_{{}_{NL}}=\frac{{\hat{x}}_3\times {\hat{x}}_1-{\hat{x}}_2^2}{{\hat{x}}_3+{\hat{x}}_1-2{\hat{x}}_2}.$$

$$x_{{}_{NL}}=x+e_{{}_{NL}};e_{{}_{NL}}\triangleq \frac{e_3\times e_1-e_2^2}{e_3+e_1-2e_2}.$$

$${\sigma }_1\triangleq {\hat{x}}_3+{\hat{x}}_1-2{\hat{x}}_2,{\sigma }_0\triangleq {\hat{x}}_2+{\hat{x}}_0-2{\hat{x}}_1,$$

$$x_{{}_{MNL}}=\{\begin{array}{l}\frac{{\hat{x}}_3\times {\hat{x}}_1-{\hat{x}}_2^2}{{\hat{x}}_3+{\hat{x}}_1-2{\hat{x}}_2};|{\sigma }_0|\le |{\sigma }_1|\\ \frac{{\hat{x}}_2\times {\hat{x}}_0-{\hat{x}}_1^2}{{\hat{x}}_2+{\hat{x}}_0-2{\hat{x}}_1};|{\sigma }_0|>|{\sigma }_1|\end{array}.$$

$$|e_1|<|e_2|<|e_3|.$$

$$e_{{}_{NL}}=\frac{(\pm {\alpha }^{-1})(\pm \alpha )\times e_2^2-e_2^2}{(\pm {\alpha }^{-1}+(\pm \alpha )-2)e_2}=\frac{(1-1)}{\pm {\alpha }^{-1}+(\pm \alpha )-2}{e}_2$$

$$e_{{}_{NL}}=\frac{({\alpha }^{-1})(-\alpha )-1}{{\alpha }^{-1}-\alpha -2}{e}_2=\frac{-2}{{\alpha }^{-1}-\alpha -2}{e}_2$$

$$\Rightarrow |\frac{e_{{}_{NL}}}{e_3}|=\frac{2}{|{\alpha }^2+2\alpha -1|}=\frac{2}{|(\alpha +1{)}^2-2|}$$

$$e_{{}_{NL}}=\frac{(-{\alpha }^{-1})(\alpha )-1}{-{\alpha }^{-1}+\alpha -2}{e}_2=\frac{-2}{-{\alpha }^{-1}+\alpha -2}{e}_2$$

$$\Rightarrow |\frac{e_{{}_{NL}}}{e_3}|=\frac{2}{|{\alpha }^2-2\alpha -1|}=\frac{2}{|(\alpha -1{)}^2-2|}$$

$$|e_1|>|e_2|>|e_3|.$$

$$|{\alpha }_1|=\alpha ,|{\alpha }_2|=(1+\delta )\alpha ;\delta >0.$$

$$|{\alpha }_1|,|{\alpha }_2|<1\Rightarrow \alpha <1,\delta <\frac{1}{\alpha }-1.$$

$$e_{{}_{NL}}=\frac{{\alpha }^{-1}\times (1+\delta )\alpha -1}{{\alpha }^{-1}+(1+\delta )\alpha -2}{e}_2=\frac{\delta \alpha }{(\alpha -1{)}^2+\delta {\alpha }^2}$$

$$\Rightarrow |\frac{e_{{}_{NL}}}{e_3}|=\frac{\delta }{(1+\delta )((\alpha -1{)}^2+\delta {\alpha }^2)}\le 1$$

$$e_{{}_{NL}}=\frac{-{\alpha }^{-1}\times -(1+\delta )\alpha -1}{-{\alpha }^{-1}-(1+\delta )\alpha -2}{e}_2=\frac{\delta \alpha }{-(\alpha +1{)}^2-\delta {\alpha }^2}$$

$$\Rightarrow |\frac{e_{{}_{NL}}}{e_3}|=\frac{\delta }{(1+\delta )((\alpha +1{)}^2+\delta {\alpha }^2)}<1;\forall \delta >0.$$

$$e_{{}_{NL}}=\frac{-{\alpha }^{-1}\times (1+\delta )\alpha -1}{{\alpha }^{-1}-(1+\delta )\alpha -2}{e}_2=\frac{-(\delta +2)\alpha }{-((\alpha +1{)}^2-2+\delta {\alpha }^2)}$$

$$\Rightarrow |\frac{e_{{}_{NL}}}{e_3}|=\frac{(\delta +2)\alpha }{(1+\delta )((\alpha +1{)}^2-2+\delta {\alpha }^2)}>1.$$

$$e_{{}_{NL}}=\frac{{\alpha }^{-1}\times -(1+\delta )\alpha -1}{-{\alpha }^{-1}+(1+\delta )\alpha -2}{e}_2=\frac{-(\delta +2)\alpha }{(\alpha -1{)}^2-2+\delta {\alpha }^2}$$

Full text

††journal: Journal of Elsevier Signal Processing

1 Introduction

Ideally, the solution of a problem is obtained either in a closed form or by an analytical approach. However, this perspective is not applicable to most cases since a closed form does not necessarily exist and even if it does, it might be too complicated to achieve. One way to approach such problems is using iterative algorithms. Despite being applicable to many cases, these algorithms have their own disadvantages. For instance, they become more complex as the number of iterations increases. Besides, their convergence/stability should be considered as well.

In order to accelerate iterative algorithms, many different methods have been proposed. Polynomial acceleration techniques are used to iteratively solve a large set of linear equations [1, 2]. The Chebyshev Algorithm (CA) for example, is a polynomial acceleration method that has been introduced to speed up the convergence rate of frame algorithms. Conjugate Gradient (CG) methods are amongst the most useful algorithms for solving optimization problems and can be simply adapted to accelerate nonlinear iterative methods such as CG-Iterative Hard Thresholding (IHT) [1, 3, 4]. Accelerating methods are mostly proposed in order to increase the convergence rate of iterative algorithms based on their target signals. As a result, each accelerating method is only capable of accelerating a limited number of iterative algorithms.

In this paper, a Non-Linear (NL) acceleration method is used to increase the convergence rate of any iterative algorithms. The proposed method is also capable of stabilizing some diverging algorithms. Previously, a similar idea but with a different point of view was used to accelerate an Iterative Method (IM) [5, 6] for non-uniform missing samples recovery problem regarding 1-D Low-Pass (LP) signals[7]. Before that, Aitken used this method to accelerate the rate of converging sequences only [8].

The NL method is capable of increasing the convergence rate of optimization algorithms. Iterative methods are widely used in gradient-based optimization algorithms such as AdaGrad for high dimensional sparse gradients [9], RMSprop for non-stationary and real-time scenarios [10] and AdaMax for a combined case of online and sparse gradient-based optimization problems [11].

Least Absolute Shrinkage and Selection Operator (LASSO) have been popularized as a regularized Least-Squares Estimation (LSE) problem [12, 13] which induces some sparsity to the LS solution. Group-Lasso was introduced to allow predefined groups of covariates to be selected into or out of a model together [14]. LSQR was proposed to solve sparse linear equations and sparse least squares; it can be used to improve LASSO solving algorithms in the case of ill-conditioned measurement matrices [15]. There are lots of algorithms for solving the Lasso problem such as Alternating Direction Method of Multipliers (ADMM) which can iteratively solve LASSO problems family[16, 17].

In this paper, after stabilizing the NL method, we extend it to accelerating image recovery algorithms as well as sparse recovery methods. We then study its interesting capability of stabilizing diverging algorithms.

In a nutshell, the present study consists of (1) improving previous works in order to accelerate and stabilize the IM with higher values of relaxation parameter, (2) accelerating iterative image recovery algorithms and (3) applying the NL method in order to accelerate iterative sparse recovery algorithms. The convergence of the proposed method is analyzed in three different categories of sub-linear, linear and super-linear converging sequences based on the sign changes of the errors in three successive estimations. These statements are confirmed by simulations of various iterative algorithms.

This paper is organized as follows: In Section 2 we briefly review some signal recovery algorithms as well as the CA. The NL algorithm is studied in Section 3. In Section 4 the simulation results are reported. Finally, in Section 5, we will conclude the paper.

2 Preliminaries

In this section, iterative algorithms are considered as a broad group of problem solving approaches and some of them are reviewed.

The IM was first proposed to compensate for the distortion caused by non-ideal interpolation. By defining G as a distortion operator it is desired to find G-1 to compensate for its distortion. The error operator could be defined as

[TABLE]

where $I$ is the identity operator. Hence we can write

[TABLE]

where $G^{-1}_{{}_{<K>}}$ is the $k^{th}$ order estimation of $G^{-1}$. It is clear that

[TABLE]

The convergence rate of the IM can be controlled by defining a relaxation parameter such as $\lambda$ in

[TABLE]

which can be recursively implemented by the equation below:

[TABLE]

where $x_{k}$ is the $k^{th}$ estimated signal. It has been proved that the IM leads to the pseudo-inverse solution and that the convergence (in the sense of stability and speed) can be controlled by tuning the relaxation parameter ($\lambda$) [18]. The IM is suitable for reconstructing band-limited signals and by choosing a proper $G$,it can be used as a non-uniform missing sample recovery algorithm [19].

Most signals are not band-limited. However, they can be sparse in some other domains. Sparse recovery is a broad problem in the literature of signal recovery. Assuming that a given signal is sparse in a specific domain, it can be perfectly reconstructed even with a limited number of observations. The main problem in sparse recovery is the minimization of an $\ell_{0}$ semi-norm minimization ($P_{0}$ problem) [20]. Because of non-convexity, the $P_{0}$ problem is usually substituted with an $\ell_{1}$ norm minimization problem ($P_{1}$ problem). , an $\ell_{1}$ norm minimization ($P_{1}$ problem) is usually substituted for the $P_{0}$. It has been shown that under some conditions regarding the signal sparsity number and the observation matrix, the solution of P0 can be obtained by solving $P_{1}$ [20, 21]. We have

[TABLE]

where $A$ and $b$ are the fat observation matrix and the observed signal, respectively.

The method of Iteratively Reweighted Least Squares (IRLS) is used to iteratively solve approximated $P_{0}$ with a weighted LS problem [22, 23]. Another approach to sparse recovery is approximating the $\ell_{0}$ semi-norm by a smooth function. Smoothed $\ell_{0}$ (SL0) method is an iterative sparse recovery method which can be used to approximate the $\ell_{0}$ semi-norm with a smooth function [24] such as

[TABLE]

where $s$ is a sparse signal with the length $N$ and $s[n]$ is its $n^{th}$ component. It can be seen that $f_{0}(s)=\|s\|_{0}$; by this approximation, the problem can be reduced to an ordinary optimization problem. As a result, its minimum can simply be found using simple minimization methods such as Steepest Descent (SD) method. It should be noted that by assigning a very small value to $\sigma$ the algorithm is trapped in a local minimum [25]. In order to escape the local minimum, the algorithm is run for a moderately larger $\sigma$ and after some iterations, the estimated signal is used for initializing the next run of the algorithm with a smaller value of $\sigma$ (it can be reduced by a decreasing factor such as Sigma Decreasing Factor ($SDF$)). This process lasts until the algorithm converges. In order to satisfy the observation constraints, after each gradient step, a projection step is required, as shown in Alg.1.

Another approach to solve sparse recovery problems is modifying inverse algorithms. In order to use the IM for sparse recovery it first needs to be properly modified. Iterative algorithms such as the IHT –which guarantees suboptimal signal recovery with robustness against additive noise [26, 27]- and Iterative Method with Adaptive Thresholding (IMAT) [28], use thresholding (in a specified transform domain) as an approach to sparse recovery. IMAT can lead to faster and better reconstruction performance compared to the IHT. Besides, the IHT needs prior knowledge regarding the sparsity number while it is not necessary for IMAT [29, 30]. IMAT can be interpreted as a modification for the IM. This can be realized by using a threshold function after each iteration of the IM, as opposed to low pass filtering (2)

[TABLE]

where $T(.)$ is the thresholding function, $Tr(.)$ and $ITr(.)$ are respectively a transformation operator and its inverse. $Tr(.)$ needs to be properly chosen in order to transform the signal to its sparsity domain. One common approach is to reduce the exponential threshold function in each iteration using the equation: $T_{k}=T_{0}\mathrm{e}^{-\alpha k}$ where $T_{0}$ and $\alpha$ are hyper-parameters used to control the threshold values in each iteration.

In non-uniform sampling, $G$ is defined as the sampling operator. Appropriate modification of $G$ can lead to improving the performance of IMAT. A modified version of IMAT is IMAT with Interpolation (IMATI) in which IMAT is improved by using a sampling operator followed by an interpolator or a smoothing function as the distortion operator, as shown in Fig.1.

IMAT and IMATI algorithms can also be simply used for image reconstruction and also 3-D signal recovery [29].

In order to increase the speed of frame algorithms, the CA method can be used. It is represented by the following equations

[TABLE]

where $A,B>0$ are the frame bounds which can control the convergence of the algorithm. Hence, inappropriate selection of these parameters can result in divergence.

3 Theory and Calculation

In this section, the importance of the convergence rate and the stability of iterative algorithms are discussed. Even though these two subjects are generally inconsistent, the NL method and its modification are introduced in order to both speed up and stabilize the iterative algorithm.

Assuming $\hat{x}_{k}[n]$ is the $k^{th}$ estimation of the desired signal at time index $n$, the corresponding recovery error is given by $e_{k}[n]=\hat{x}_{k}[n]-x[n]$. Generally, $e_{k}[n]$ can be written as a proportion of $e_{k-1}[n]$,

[TABLE]

where $\alpha_{k}[n]$ is the coefficient of proportionality. To be concise and in order not to lose generality, we consider three typical estimated signals, such as $\hat{x}_{1}$, $\hat{x}_{2}$ and $\hat{x}_{3}$, for a specific time index. By assigning the same value to the first two successive $\alpha_{k}$’s (i.e., $\alpha_{1}=\alpha_{2}=\alpha$), we can write

[TABLE]

By computing $x$ from (3) the following Non-Linear (NL) formula is obtained:

[TABLE]

Considering $\hat{x}_{i}=x+e_{i}$ for $i=1,2,3$, it can be deduced from the NL formula that

[TABLE]

Due to the sign alternations of $e_{i}$’s, it can be shown that in some cases $|e_{{}_{NL}}|$ can be larger than $|e_{3}|$ -which is probably smaller than both $|e_{2}|$ and $|e_{1}|$-. In other words, the error obtained while using the NL algorithm is larger than the errors of the existing estimations. In order to analyze the convergence of the NL method, it should be noted that there exist $2^{3}=8$ possible cases for the sign alternations of $e_{i}$’s. Assuming that $sign(e_{2})=1$, this number can be reduced to $2^{2}=4$. Furthermore, there are two cases based on the relative changes of $|\alpha_{1}|$ and $|\alpha_{2}|$. The first case is linear convergence (for $|\alpha_{1}|=|\alpha_{2}|$) while we consider the second case to serve as the two cases of sub-linear and super-linear convergence (for $|\alpha_{1}|\neq|\alpha_{2}|$).

In order to make sure that the estimation does not diverge when divided by zero, the NL method is applied to $\hat{x}_{0}$, $\hat{x}_{1}$ and $\hat{x}_{2}$ [7]. Therefore, by defining

[TABLE]

the Modified NL (MNL) method is obtained and can be represented as follows

[TABLE]

Apart from that, the error is unavoidable for very small values of $\sigma_{i}$ due to the finite precision hardware implementation. Fortunately, the latter divergence occurs only in a very small number of points of the NL estimation and is usually noticeable as high spikes.

Different simple techniques can be used in order to compensate for these issues. One approach could be using some of the linear combinations (weighted averages) of the existing estimations in substitution instead. This can lead to reducing the sign alternations of the errors. Another approach is applying the NL formula to transformed versions of the estimations in another specified domain. Assuming that the signals are bounded, simple methods such as Clipping, Substitution and Smoothing can be used in order to compensate for the undesirable spikes caused by the NL method. For each up-crossing of the foreknown levels (the maximum and the minimum) we either clip the estimated signal or substitute the best existing signal for the reconstructed signal. If the signal is too noisy, the Median Filter (MedFilt) can be used in order to smooth the NL signal. In this article, we use both MedFilt and Clipping.

One interesting feature of the NL formula is its symmetry. For a diverging algorithm we have

[TABLE]

For any selection of $\hat{x}_{1}$, $\hat{x}_{2}$ and $\hat{x}_{3}$, whether the error is decreasing or increasing, the NL formula leads to the same results since $\hat{x}_{3}$, $\hat{x}_{2}$ and $\hat{x}_{1}$ can be considered to be three successive converging estimations. For a constant rate of changes for three successive errors (either converging or diverging), using the NL method can lead to a perfect signal reconstruction. This property can be used in Watermarking and Steganography.

4 Simulation Results and Discussion

In this section we use the MNL method to accelerate some iterative algorithms. There are so many algorithms which solve the problems recursively. Even though the introduced NL method is capable of improving most of the existing iterative algorithms under some simple assumptions, we focus on recovering the signal from its with missing samples (special case of nonuniform samples). Missing samples occur at random indices with iid Bernoulli distribution with the parameter $p=\text{Loss-Rate (}LR$).

In order to illustrate the performance of image recovery algorithms, some Image Quality Assessment (IQA) methods are used. Peak Signal to Noise Ratio (PSNR), Inter-Patch and Intra-Patch based IQA (IIQA) [31], Structural Similarity (SSIM) [32] (Convex SIMilarity (CSIM) as its convex version [33]), Multi-Scale SSIM (MS-SSIM) [34], Edge Strength Similarity (ESSIM) [35] and Feature Similarity (FSIM) [36] are used as Full-Reference (FR) IQA methods while Blind/Referenceless Image Spatial Quality Evaluator (BRISQUE) [37, 38] and Naturalness Image Quality Evaluator (NIQE) [39] are used as No-Reference (NR) IQA methods.

The MNL formula was applied to the IM in order to reconstruct the 1D sampled signals (with the length $L=500$) with a specified Over Sampling Ratio (OSR). The LP signals were generated by filtering the white normal noise using the DFT Filter. In order to achieve fair results, the result of each experiment was averaged over 100 runs.

For a selection of parameters same as the one in [7] ($LR=33\%$ and $OSR=8$), the MNL can be stabilized by using Substitution and Clipping. Then, it can be applied to accelerate the IM even when $\lambda>1$, as shown in Fig.2.

By increasing and decreasing the LR and the OSR, respectively, the MNL starts to act unstably. This is because the chances for having a diverging case grow. Fortunately, simulation results show that the MNL estimated signal includes only a very few unstable points and therefore, can be stabilized using MedFilt and Substitution, as shown in Fig.3. Note that the MNL improves iterative algorithms in terms of convergence. Hence, what actually leads to lower SNR improvement in this experiment is the performance of the IM.

Increasing the relaxation parameter causes the IM to diverge. In order to avoid the latter problem, the MNL method can be used, as shown in Fig.4.

The MNL method can be easily generalized to improve iterative image recovery algorithms. It can be applied to the IM in order to reconstruct the image “Lenna” (with the size $512\times 512$) from its nonuniform samples, as shown in Fig.5, 6.

The MNL was used to stabilize image recovery using the IM algorithm, as shown in Fig.7.

Applying the MNL formula to the CA eventuates in the same results. Therefore, we focus on the stabilizing property of the MNL, as shown in Fig.8, 9.

There are two approaches to apply the MNL method to the SL0:

Applying the MNL formula to the last four estimations of the inner loop (MNL) 2. 2.

Applying the MNL formula to the last estimations of the main algorithm ($\text{MNL}_{2}$)

The signal is assumed to be sparse in the DFT domain. Non-zero components are independently generated at random indices with the probability $p_{{}_{NZ}}$ while zero components are assumed to be contaminated by zero-mean Gaussian noise (with the standard deviation $\sigma_{off}$). The algorithm is initialized using the Least Square Estimation (LSE) and $\sigma_{0}=2\times\text{max}(|\hat{x}_{0}[n]|)$). $\sigma_{k}$ decreases by $SDF$ with each iteration. The inner minimization loop uses the SD method with 3 iterations and the adaption rate $\mu_{0}=2$. Fig.10 and Fig.11 show the results with and without the presence of noise, respectively.

The MNL formula can be used to improve the performance of IMAT, as shown in Fig.12, 13. The latter figures show the performance curves of reconstructing the ”Baboon” (also known as the ”Mandrill”) image from its nonuniform samples.

Considering IMAT’s multiple parameters, it seems rather difficult to assess its sensitivity to parameter changes. To do so, the performance curves are depicted in terms of $\lambda$ and different values of $\alpha$. DCT was used as the function Tr(.) and the main algorithm was run in 5 iterations, as shown in Fig.14,15.

As it can be seen, the MNL can preserve the performance of the algorithm to a great extent. Also, the stability range of the algorithm is extended. A wide range of stability is important because of its changes for different images and sampling patterns since in each case there is an optimum value of $\lambda$ for which the algorithm is both stable and fast enough. Crossing that optimal point causes the algorithm to diverge. Hence, by extending the stability range the reliability of the main algorithm increases. A highly detailed image of an eye (with the size $1600\times 1030$) is reconstructed from its nonuniform samples as an example of sparse image recovery, as shown in Fig.16.

IMATI’s sensitivity to parameter changes can be studied, as shown in Fig.17,18.

As shown in Fig. 19, the MNL can be used to accelerate the ADMM algorithm for solving LASSO problems family with $\alpha,\rho,\lambda$ and $K$ as the over-relaxation, augmented Lagrangian, Lagrangian parameters and group size, respectively.

It must be mentioned that in the case of $\frac{m}{n}<0.5$, our simulations show that the MNL is not able to improve the ADMM; in fact, based on the Convergence Analysis, corresponding cases to $e_{1}\times e_{3}<0$ randomly occur and the NL diverges.

The convergence of the IRLS can be increased by MNL, as shown in Fig.20.

5 Conclusion

In this paper, a non-linear acceleration method (the NL) and its modification (the MNL) are introduced in order to accelerate and improve iterative algorithms. Besides, a complete analysis on convergence is given. It is stated that the proposed method can improve a wide variety of linearly convergent algorithms including optimization methods, band-limited signal recovery methods and sparse recovery algorithms. The proposed method is also capable of stabilizing diverging algorithms. It is shown that the MNL method improves iterative sparse recovery algorithms such as the IRLS, ADMM, SL0 and IMAT. Simulation results show that the performance of this method in terms of various quality assessments is noticeably better. By stabilizing and accelerating the CA method, it is shown that the MNL method can even be used to improve an acceleration method.

Acknowledgement

We wish to thank Nasim Bagheri for her assistance with the text editing.

Appendix A Convergence Analysis

At first, assume that $|\alpha_{1}|=|\alpha_{2}|$ (linear convergence):

Cases $1.1$ and $1.2$, if $\alpha_{1}=\alpha_{2}=\pm\alpha$ ($\alpha>0$) then

[TABLE]

which is equal to zero as expected.

Case 1.3, if $\alpha_{1}=-\alpha_{2}=\alpha>0$ then

[TABLE]

which means that for a convergent sequence ($0<\alpha<1$), the NL estimation diverges; in other words, we have $|\frac{e_{{}_{NL}}}{e_{3}}|>1$.

Case 1.4, if $-\alpha_{1}=\alpha_{2}=\alpha>0$ then

[TABLE]

which also results in the divergence of the NL method in the case of a converging sequence ($0<\alpha<1$).

In a more general case, it can be assumed that $|\alpha_{1}|\neq|\alpha_{2}|$. Hence, $e_{1}$ and $e_{3}$ should be rewritten in terms of $e_{2}$. In a converging algorithm, subsequent estimations satisfy

[TABLE]

We can assume that after each iteration, the algorithm becomes less capable of reducing the errors which results in

[TABLE]

The latter alongside the convergence of the algorithm leads to obtaining

[TABLE]

Hence, we have the following, which represent the case of sub-linear convergence:

Case 2.1, if all of the $e_{i}$’s have the same sign, in other words, if $e_{1}\times e_{2}>0$ and $e_{3}\times e_{2}>0$, then

[TABLE]

where the equality holds for $\delta=\frac{1}{\alpha}-1$; this implies that the NL method does not improve the existing estimation, as expected from (5) and (6).

Case 2.2, if $e_{1}\times e_{2}<0$ and $e_{3}\times e_{2}<0$ then

[TABLE]

Case 2.3, if $e_{1}\times e_{2}>0,e_{3}\times e_{2}<0$ and $0<\delta<\frac{1}{\alpha}-1$ then

[TABLE]

Case 2.4, if $e_{1}\times e_{2}<0$ and $e_{3}\times e_{2}>0$ then

[TABLE]

where $0<\delta<\frac{1}{\alpha}-1$.

In the case of super-linear convergent algorithms, it can be assumed that

[TABLE]

which results in the same procedure as for the case of sub-linear convergence except for the sign of $\delta$.

Case 3.1, if all of the $e_{i}$’s have the same sign, then

[TABLE]

where $\delta_{0}=\frac{\alpha^{2}-\alpha+1-\sqrt{2\alpha^{2}-2\alpha+1}}{\alpha^{2}}$. Hence, the NL estimation diverges for $\delta_{0}<\delta<1$. It can be seen that for a cubically convergent sequence ($\delta=1-\alpha$), the NL estimation diverges since $\delta_{0}<1-\alpha$.

Case 3.2, if $e_{1}\times e_{2}<0$ and $e_{3}\times e_{2}<0$, then

[TABLE]

where $\delta_{0}=\frac{\alpha^{2}+\alpha+1-\sqrt{2\alpha^{2}+2\alpha+1}}{\alpha^{2}}$. In this case, sequences with quadratic convergence can be accelerated if $\delta_{0}>1-\alpha$. In other words, these sequences can be accelerated if $0.36110<\alpha$.

Case 3.3, if $e_{1}\times e_{2}>0,e_{3}\times e_{2}<0$ and $0<\delta<1$ then

[TABLE]

Case 3.4, if $e_{1}\times e_{2}<0$ and $e_{3}\times e_{2}>0$ then

[TABLE]