A Nonlinear Acceleration Method for Iterative Algorithms

A Nonlinear Acceleration Metho d for Iterativ e Algorithms Mahdi Shamsi a , Mahmoud Ghandi , F arokh Marv asti a a Multime dia and Signal pr o c essing L ab (MSL), A dvanc e d Communic ations R ese arch Institute (ACRI), EE Dep artment of Sharif University of T echnolo gy, T ehran, I.R.Ir an Abstract Iterativ e metho ds hav e led to b etter understanding and solving problems suc h as missing sampling, deconv olution, inv erse systems, impulsiv e and Salt and P epp er noise remo v al problems. Ho wev er, the challenges suc h as the speed of con vergence and or the accuracy of the answer still remain. In order to improv e the existing iterative algorithms, a non-linear metho d is discussed in this pap er. The mentioned metho d is analyzed from differen t aspects, including its conv er- gence and its abilit y to accelerate recursive algorithms.W e sho w that this metho d is capable of impro ving Iterativ e Metho d (IM) as a non-uniform sampling re- construction algorithm and some iterative sparse recov ery algorithms suc h as Iterativ e Reweigh ted Least Squares (IRLS), Iterative Metho d with Adaptive Thresholding (IMA T), Smo othed ` 0 (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 b oth accelerating and stabilizing the well-kno wn Chebyshev Acceleration (CA) metho d. F urther- more, the prop osed algorithm can extend the stability range by reducing the sensitivit y of iterative algorithms to the changes of adaptation rate. Keywor ds: Non-Linear Acceleration; Iterativ e Metho ds; Sparse Reco very; Acceleration Metho ds; IMA T; LASSO. 1. In tro duction Ideally , the solution of a problem is obtained either in a closed form or by an analytical approac h. Ho wev er, this persp ective is not applicable to most cases Pr eprint submitte d to Journal of Elsevier Signal Pr oc essing Septemb er 23, 2024 since a closed form do es not necessarily exist and even if it do es, it might be to o complicated to ac hieve. One wa y to approac h suc h problems is using iterativ e algorithms. Despite being applicable to many cases, these algorithms ha ve their o wn disadv an tages. F or instance, they b ecome more complex as the num b er of iterations increases. Besides, their conv ergence/stabilit y should b e considered as w ell. In order to accelerate iterative algorithms, man y different metho ds hav e b een prop osed. Polynomial acceleration techniques are used to iteratively solve a large set of linear equations [1, 2]. The Chebyshev Algorithm (CA) for exam- ple, is a p olynomial acceleration metho d that has b een introduced to sp eed up the conv ergence rate of frame algorithms. Conjugate Gradient (CG) metho ds are amongst the most useful algorithms for solving optimization problems and can b e simply adapted to accelerate nonlinear iterative metho ds such as CG- Iterativ e Hard Thresholding (IHT) [1, 3, 4]. Accelerating metho ds are mostly prop osed in order to increase the conv ergence rate of iterative algorithms based on their target signals. As a result, each accelerating metho d is only capable of accelerating a limited n umber of iterative algorithms. In this pap er, a Non-Linear (NL) acceleration metho d is used to increase the conv ergence rate of any iterative algorithms. The prop osed method is also capable of stabilizing some diverging algorithms. Previously , a similar idea but with a different p oin t of view was used to accelerate an Iterativ e Metho d (IM) [5, 6] for non-uniform missing samples reco v ery problem regarding 1-D Lo w-Pass (LP) signals[7]. Before that, Aitk en used this metho d to accelerate the rate of con verging sequences only [8]. The NL metho d is capable of increasing the conv ergence rate of optimization algorithms. Iterative metho ds are widely used in gradient-based optimization algorithms suc h as AdaGrad for high dimensional sparse gradients [9], RMSprop for non-stationary and real-time scenarios [10] and AdaMax for a com bined case of online and sparse gradien t-based optimization problems [11]. Least Absolute Shrink age and Selection Op erator (LASSO) ha v e b een popu- larized as a regularized Least-Squares Estimation (LSE) problem [12, 13] which 2 induces some sparsity to the LS solution. Group-Lasso was introduced to allo w predefined groups of co v ariates to b e selected in to or out of a mo del together [14]. LSQR was prop osed to solve sparse linear equations and sparse least squares; it can b e used to improv e LASSO solving algorithms in the case of ill-conditioned measuremen t matrices [15]. There are lots of algorithms for solving the Lasso problem such as Alternating Direction Metho d of Multipliers (ADMM) whic h can iterativ ely solve LASSO problems family[16, 17]. In this pap er, after stabilizing the NL metho d, we extend it to accelerating image reco very algorithms as well as sparse recov ery metho ds. W e then study its in teresting capability of stabilizing diverging algorithms. In a nutshell, the present study consists of (1) improving previous w orks in order to accelerate and stabilize the IM with higher v alues of relaxation parameter, (2) accelerating iterativ e image recov ery algorithms and (3) applying the NL metho d in order to accelerate iterativ e sparse recov ery algorithms. The con vergence of the proposed metho d is analyzed in three different categories of sub-linear, linear and sup er-linear con v erging sequences based on the sign c hanges of the errors in three successive estimations. These statements are confirmed b y simulations of v arious iterativ e algorithms. This paper is organized as follo ws: In Section 2 w e briefly review some signal reco very algorithms as well as the CA. The NL algorithm is studied in Section 3. In Section 4 the simulation results are rep orted. Finally , in Section 5, we will conclude the pap er. 2. Preliminaries In this section, iterativ e algorithms are considered as a broad group of prob- lem solving approac hes 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 op erator it is desired to find G-1 to comp ensate for its distortion. The error op erator could b e defined as E , I − G 3 where I is the identit y operator. Hence we can write G − 1 = I G = I I − E = I + ∞ X n =1 E n ; k E k < 1 ⇒ G − 1 = I + K X n =1 E n where G − 1 is the k th order estimation of G − 1 . It is clear that G − 1 = I + E ( G − 1 ) . The con vergence rate of the IM can b e controlled by defining a relaxation pa- rameter suc h as λ in G − 1 = λI λG = λI I − E λ ; E λ , I − λG , k E λ k < 1 ⇒ G − 1 = I + K X n =1 E λ n whic h can b e recursiv ely implemented by the equation b elow: x k = λ ( x 0 − G ( x k − 1 )) + x k − 1 (1) where x k is the k th estimated signal. It has b een prov ed that the IM leads to the pseudo-inv erse solution and that the conv ergence (in the sense of stability and sp eed) can b e controlled by tuning the relaxation parameter ( λ ) [18]. The IM is suitable for reconstructing band-limited signals and by choosing a proper G ,it can b e used as a non-uniform missing sample recov ery algorithm [19]. Most signals are not band-limited. Ho wev er, they can b e sparse in some other domains. Sparse reco very is a broad problem in the literature of signal reco very . Assuming that a given signal is sparse in a sp ecific domain, it can b e p erfectly reconstructed ev en with a limited num b er of observ ations. The main problem in sparse reco very is the minimization of an ` 0 semi-norm minimiza- tion ( P 0 problem) [20]. Because of non-con vexit y , the P 0 problem is usually substituted with an ` 1 norm minimization problem ( P 1 problem). , an ` 1 norm minimization ( P 1 problem) is usually substituted for the P 0 . It has been sho wn 4 that under some conditions regarding the signal sparsity num b er and the ob- serv ation matrix, the solution of P0 can b e obtained by solving P 1 [20, 21]. W e ha ve P 0 : min s k s k 0 ; As = b , P 1 : min s k s k 1 ; As = b where A and b are the fat observ ation matrix and the observ ed signal, resp ec- tiv ely . The method of Iteratively Reweigh ted Least Squares (IRLS) is used to it- erativ ely solve approximated P 0 with a weigh ted LS problem [22, 23]. Another approac h to sparse recov ery is appro ximating the ` 0 semi-norm b y a smo oth function. Smo othed ` 0 (SL0) method is an iterative sparse recov ery method whic h can b e used to appro ximate the ` 0 semi-norm with a smooth function [24] suc h as f σ ( s ) = N − N X n =1 F σ ( s [ n ]) ; F σ ( s [ n ]) = exp( − | s [ n ] | 2 2 σ 2 ) where s is a sparse signal with the length N and s [ n ] is its n th comp onen t. It can b e seen that f 0 ( s ) = k s k 0 ; b y this approximation, the problem can be reduced to an ordinary optimization problem. As a result, its minimum can simply b e found using simple minimization metho ds suc h as Steep est Descent (SD) method. It should b e noted that by assigning a very small v alue to σ the algorithm is trapp ed in a lo cal minimum [25]. In order to escap e the lo cal minim um, the algorithm is run for a moderately larger σ and after some iterations, the estimated signal is used for initializing the next run of the algorithm with a smaller v alue of σ (it can be reduced b y a decreasing factor such as Sigma Decreasing F actor ( S D F )). This pro cess lasts until the algorithm conv erges. In order to satisfy the observ ation constraints, after eac h gradient step, a pro jection step is required, as sho wn in Alg.1. Another approach to solve sparse recov ery problems is mo difying inv erse al- gorithms. In order to use the IM for sparse recov ery it first needs to b e prop erly mo dified. Iterativ e algorithms such as the IHT which guarantees sub optimal signal recov ery with robustness against additive noise [26, 27]- and Iterative 5 Data: A, b, σ 0 , M , K , r , µ 0 Result: s s = A + b ; for m=0:M-1 do σ ← r m × σ 0 ; for k=0:K-1 do G : s ← s − µ 0 ∇ f σ ( s ) P : s ← s − A + ( As − b ) end end Algorithm 1: The SL0 algorithm. Metho d with Adaptive Thresholding (IMA T) [28], use thresholding (in a sp ec- ified transform domain) as an approach to sparse recov ery . IMA T can lead to faster and b etter reconstruction p erformance compared to the IHT. Besides, the IHT needs prior knowledge regarding the sparsit y num ber while it is not necessary for IMA T [29, 30]. IMA T can b e interpreted as a mo dification for the IM. This can b e realized b y using a threshold function after each iteration of the IM, as opp osed to low pass filtering (2) x k = T [ λ ( x 0 − G ( x k − 1 )) + x k − 1 ] (2) T ( x ) =                    X = T r ( x ) Y =      X ; T k ≤ | X | 0 ; else return I T r ( Y ) where T ( . ) is the thresholding function, T r ( . ) and I T r ( . ) are resp ectively a transformation op erator and its inv erse. T r ( . ) needs to b e prop erly chosen in order to transform the signal to its sparsity domain. One common approach is to reduce the exp onential threshold function in each iteration using the equation: T k = T 0 e − αk where T 0 and α are h yp er-parameters used to con trol the threshold v alues in each iteration. 6 In non-uniform sampling, G is defined as the sampling operator. Appropriate mo dification of G can lead to impro ving the performance of IMA T. A mo dified v ersion of IMA T is IMA T with Interpolation (IMA TI) in which IMA T is im- pro ved b y using a sampling op erator follo wed b y an interpolator or a smo othing function as the distortion op erator, as shown in Fig.1. Figure 1: Blo c k Diagram of IMA TI algorithm with distortion op erator G and relaxation parameter of λ . IMA T and IMA TI algorithms can also b e simply used for image reconstruc- tion and also 3-D signal reco very [29]. In order to increase the speed of frame algorithms, the CA metho d can b e used. It is represen ted by the following equations x 1 = 2 A + B x 0 , λ k = (1 − ρ 2 4 λ k − 1 ) − 1 ; for k > 1 : x k = ( x 1 + x k − 1 − 2 A + B G ( x k − 1 ) − x k − 2 ) λ k + x k − 2 where A, B > 0 are the frame b ounds which can con trol the conv ergence of the algorithm. Hence, inappropriate selection of these parameters can result in div ergence. 3. Theory and Calculation In this section, the imp ortance of the conv ergence rate and the stability of iterativ e algorithms are discussed. Ev en though these t wo sub jects are generally inconsisten t, the NL metho d and its modification are in tro duced in order to both sp eed up and stabilize the iterative algorithm. Assuming ˆ x k [ n ] is the k th estimation of the desired signal at time index n , the corresp onding recov ery error is giv en by e k [ n ] = ˆ x k [ n ] − x [ n ]. Generally , 7 e k [ n ] can b e written as a prop ortion of e k − 1 [ n ], e k [ n ] = α k − 1 [ n ] × e k − 1 [ n ] ; e k − 1 [ n ] 6 = 0 where α k [ n ] is the co efficient of prop ortionalit y . T o b e concise and in order not to lose generalit y , w e consider three typical estimated signals, suc h as ˆ x 1 , ˆ x 2 and ˆ x 3 , for a sp ecific time index. By assigning the same v alue to the first t wo successiv e α k ’s (i.e., α 1 = α 2 = α ), we can write ˆ x 3 − x = α ( ˆ x 2 − x ) = α 2 ( ˆ x 1 − x ) . (3) By computing x from (3) the follo wing Non-Linear (NL) formula is obtained: x N L = ˆ x 3 × ˆ x 1 − ˆ x 2 2 ˆ x 3 + ˆ x 1 − 2 ˆ x 2 . Considering ˆ x i = x + e i for i = 1 , 2 , 3, it can b e deduced from the NL form ula that x N L = x + e N L ; e N L , e 3 × e 1 − e 2 2 e 3 + e 1 − 2 e 2 . Due to the sign alternations of e i ’s, it can b e shown that in some cases | e N L | can b e larger than | e 3 | -whic h 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 con vergence of the NL metho d, it should b e noted that there exist 2 3 = 8 p ossible cases for the sign alternations of e i ’s. Assuming that sig n ( e 2 ) = 1, this num b er can b e reduced to 2 2 = 4. F urthermore, there are t wo cases based on the relative changes of | α 1 | and | α 2 | . The first case is linear con vergence (for | α 1 | = | α 2 | ) while we consider the second case to serve as the t wo cases of sub-linear and sup er-linear con vergence (for | α 1 | 6 = | α 2 | ). In order to mak e sure that the estimation do es not div erge when divided by zero, the NL metho d is applied to ˆ x 0 , ˆ x 1 and ˆ x 2 [7]. Therefore, b y defining σ 1 , ˆ x 3 + ˆ x 1 − 2 ˆ x 2 , σ 0 , ˆ x 2 + ˆ x 0 − 2 ˆ x 1 , the Mo dified NL (MNL) metho d is obtained and can b e represented as follows x M N L =      ˆ x 3 × ˆ x 1 − ˆ x 2 2 ˆ x 3 + ˆ x 1 − 2 ˆ x 2 ; | σ 0 | ≤ | σ 1 | ˆ x 2 × ˆ x 0 − ˆ x 2 1 ˆ x 2 + ˆ x 0 − 2 ˆ x 1 ; | σ 0 | > | σ 1 | . (4) 8 Apart from that, the error is unav oidable for v ery small v alues of σ i due to the finite precision hardware implemen tation. F ortunately , the latter div ergence o ccurs only in a v ery small num b er of p oin ts of the NL estimation and is usually noticeable as high spik es. Differen t simple techniques can b e used in order to compensate for these issues. One approach could be using some of the linear com binations (weigh ted a verages) of the existing estimations in substitution instead. This can lead to reducing the sign alternations of the errors. Another approac h is applying the NL formula to transformed versions of the estimations in another sp ecified do- main. Assuming that the signals are b ounded, simple methods suc h as Clipping, Substitution and Smo othing can b e used in order to comp ensate for the unde- sirable spikes caused by the NL metho d. F or each up-crossing of the foreknown lev els (the maxim um and the minimum) we either clip the estimated signal or substitute the b est existing signal for the reconstructed signal. If the signal is to o noisy , the Median Filter (MedFilt) can b e used in order to smooth the NL signal. In this article, we use b oth MedFilt and Clipping. One interesting feature of the NL formula is its symmetry . F or a div erging algorithm w e hav e | e 1 | < | e 2 | < | e 3 | . F or any selection of ˆ x 1 , ˆ x 2 and ˆ x 3 , whether the error is decreas ing or increasing, the NL form ula leads to the same results since ˆ x 3 , ˆ x 2 and ˆ x 1 can b e considered to be three successive con v erging estimations. F or a constant rate of c hanges for three successiv e errors (either conv erging or diverging), using the NL metho d can lead to a p erfect signal reconstruction. This prop ert y can b e used in W a- termarking and Steganograph y . 4. Sim ulation Results and Discussion In this section we use the MNL metho d to accelerate some iterativ e algo- rithms. There are so many algorithms whic h solve the problems recursively . 9 Ev en though the in tro duced NL method is capable of improving most of the existing iterative algorithms under some simple assumptions, w e focus on re- co vering the signal from its with missing samples (sp ecial case of nonuniform samples). Missing samples o ccur at random indices with iid Bernoulli distribu- tion with the parameter p = Loss-Rate ( LR ). In order to illustrate the performance of image recov ery algorithms, some Image Quality Assessment (IQA) metho ds are used. Peak Signal to Noise Ra- tio (PSNR), Inter-P atch and In tra-Patc h based IQA (I IQA) [31], Structural Similarit y (SSIM) [32] (Conv ex SIMilarity (CSIM) as its conv ex v ersion [33]), Multi-Scale SSIM (MS-SSIM) [34], Edge Strength Similarity (ESSIM) [35] and F eature Similarity (FSIM) [36] are used as F ull-Reference (FR) IQA metho ds while Blind/Referenceless Image Spatial Qualit y Ev aluator (BRISQUE) [37, 38] and Naturalness Image Quality Ev aluator (NIQE) [39] are used as No-Reference (NR) IQA metho ds. 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 achiev e fair results, the result of eac h exp erimen t w as av eraged o ver 100 runs. F or a selection of parameters same as the one in [7] ( LR = 33% and O S R = 8), the MNL can b e stabilized by using Substitution and Clipping. Then, it can b e applied to accelerate the IM even when λ > 1, as sho wn in Fig.2. (a) The MNL+Clipping. (b) The MNL+Substitution. Figure 2: SNR curves of the IM and the MNL, O S R = 8, LR = 1 3 , λ = 2 . 2. 10 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. F ortunately , simulation results show that the MNL estimated signal includes only a very few unstable p oin ts and therefore, can be stabilized using MedFilt and Substitution, as sho wn in Fig.3. Note that the MNL improv es iterative algorithms in terms of conv ergence. Hence, what actually leads to lo wer SNR impro vemen t in this experiment is the performance of the IM. (a) The MNL+MedFilt. (b) The MNL+Substitution. Figure 3: SNR curves of the IM and the MNL, O S R = 4, LR = 1 2 , λ = 2. Increasing the relaxation parameter causes the IM to diverge. In order to av oid the latter problem, the MNL metho d can b e used, as shown in Fig.4. (a) The MNL+Clipping. (b) The MNL+Substitution. Figure 4: SNR curves of the IM and the MNL, O S R = 8, LR = 1 3 , λ = 2 . 2. The MNL metho d can b e easily generalized to improv e iterative image recov ery algorithms. It can be applied to the IM in order to reconstruct the image Lenna (with the size 512 × 512) from its non uniform samples, as shown in Fig.5, 6. 11 Figure 5: PSNR and IIQA curves of the IM and the MNL(+Clipping), λ = 2, O S R = 4, LR = 1 3 , image: Lenna. Figure 6: SSIM and MS-SSIM curves of the IM and the MNL(+Clipping), λ = 2, O S R = 4, LR = 1 3 , image: Lenna. The MNL w as used to stabilize image reco v ery using the IM algorithm, as shown in Fig.7. Applying the MNL formula to the CA ev en tuates in the same results. Therefore, w e fo cus on the stabilizing prop ert y of the MNL, as sho wn in Fig.8, 9. There are t wo approaches to apply the MNL metho d to the SL0: I Applying the MNL formula to the last four estimations of the inner lo op (MNL) I I Applying the MNL form ula to the last estimations of the main algorithm (MNL 2 ) The signal is assumed to be sparse in the DFT domain. Non-zero components are indep endently generated at random indices with the probability p N Z while zero comp onen ts are assumed to b e contaminated by zero-mean Gaussian noise (with the standard deviation σ of f ). The algorithm is initialized using the Least 12 Figure 7: PSNR, IIQA, SSIM, and MS-SSIM curves of the IM and the MNL(+Clipping), λ = 3 . 5, O S R = 4, LR = 2 3 , image: Lenna. Square Estimation (LSE) and σ 0 = 2 × max( | ˆ x 0 [ n ] | )). σ k decreases by S D F with each iteration. The inner minimization lo op uses the SD metho d with 3 iterations and the adaption rate µ 0 = 2. Fig.10 and Fig.11 sho w the results with and without the presence of noise, resp ectively . The MNL formula can b e used to improv e the p erformance of IMA T, as shown in Fig.12, 13. The latter figures sho w the p erformance curv es of reconstructing the ”Bab oon” (also kno wn as the ”Mandrill”) image from its non uniform samples. Considering IMA T’s multiple parameters, it seems rather difficult to assess its sensitivit y to parameter changes. T o do so, the p erformance curv es are depicted in terms of λ and different v alues of α . DCT was used as the function T r(.) and the main algorithm w as 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 imp ortan t b ecause of its changes for different images and sampling patterns since in eac h case there is an optimum v alue of λ for whic h the 13 Figure 8: PSNR and I IQA curves of the IM and the MNL(+Clipping), A = 0 . 25, B = 0 . 6, λ 0 = 3 . 5, LR = 1 2 , image: Cameraman. Figure 9: SSIM and MS-SSIM curves of the IM and the MNL(+Clipping), A = 0 . 25, B = 0 . 6, λ 0 = 3 . 5, O S R = 2, LR = 1 2 , image: Cameraman. algorithm is b oth stable and fast enough. Crossing that optimal p oint causes the algorithm to diverge. Hence, by extending the stability range the reliabilit y of the main algorithm increases. A highly detailed image of an ey e (with the size 1600 × 1030) is reconstructed from its nonuniform samples as an example of sparse image reco very , as shown in Fig.16. IMA TI’s sensitivity to parameter c hanges can be studied, as shown in Fig.17,18. As shown in Fig. 19, the MNL can b e used to accelerate the ADMM algorithm for solving LASSO problems family with α , ρ, λ and K as the o ver-relaxation, augmen ted Lagrangian, Lagrangian parameters and group size, resp ectively . It m ust b e mentioned that in the case of m n < 0 . 5, our simulations show that the MNL is not able to improv e the ADMM; in fact, based on the Con vergence Analysis, corresp onding cases to e 1 × e 3 < 0 randomly o ccur and the NL div erges. The con vergence of the IRLS can b e increased by MNL, as shown in Fig.20. 14 (a) LR = 30%. (b) LR = 50%. Figure 10: SNR curves of the SL0 recovery metho d and the MNL, 1 D signal with L = 1000, S D F = 0 . 5, p N Z = 0 . 05, av eraged ov er 50 runs. 5. Conclusion In this pap er, a non-linear acceleration metho d (the NL) and its mo difi- cation (the MNL) are introduced in order to accelerate and impro ve iterative algorithms. Besides, a complete analysis on conv ergence is giv en. It is stated that the prop osed method can improv e a wide v ariety of linearly con v ergent algorithms including optimization metho ds, band-limited signal recov ery meth- o ds and sparse reco v ery algorithms. The proposed method is also capable of stabilizing diverging algorithms. It is shown that the MNL method impro v es iterativ e sparse recov ery algorithms such as the IRLS, ADMM, SL0 and IMA T. Sim ulation results show that the performance of this metho d in terms of v arious qualit y assessments is noticeably b etter. By stabilizing and accelerating the CA metho d, it is shown that the MNL metho d can even b e used to improv e an acceleration metho d. 15 (a) S D F = 0 . 5. (b) S D F = 0 . 7. Figure 11: SNR curves of the SL0 recovery metho d and the MNL, 1 D signal with L = 1000, LR = 30%, p N Z = 0 . 1, av eraged ov er 50 runs. Figure 12: PSNR and I IQA curves of IMA T and the MNL(+MedFilt), λ = 2, T 0 = 300, α = 1, LR = 30%, image: Babo on. Figure 13: SSIM and MS-SSIM curves of IMA T and the MNL(+MedFilt), λ = 2, T 0 = 300, α = 1, LR = 30%, image: Babo on. 16 (a) LR = 30%. (b) LR = 50%. Figure 14: PSNR curves of the IMA T and the MNL (+MedFilt), image: Pirate. 17 (a) LR = 30%. (b) LR = 50%. Figure 15: SSIM curves of the IMA T and the MNL (+MedFilt), image: Pirate. 18 (a) Main image, sampled image, and reconstructed images using IMA TI and the MNL (b) PSNR, I IQA, NIQE, and BRISQUE curves. (c) SSIM, MS-SSIM, ESSIM, and FSIM curves. Figure 16: Image recov ery using IMA TI (+Gaussian interpretor) and the MNL (+Clipping), σ = 2, LR = 30%, T 0 = 1000, α = 1, λ = 3 . 5. 19 (a) LR = 30%. (b) LR = 50%. Figure 17: PSNR curves of IMA TI and the MNL (+Clipping), image: Mandril. 20 (a) LR = 30%. (b) LR = 50%. Figure 18: SSIM curves of IMA TI and the MNL (+Clipping), image: Mandril. 21 (a) α = 0 . 3 , ρ = 0 . 8. (b) α = 0 . 3 , ρ = 0 . 5. (c) α = 0 . 3 , ρ = 0 . 5. Figure 19: SNR curves of ADMM for solving LASSO, LASSO-LSQR and group LASSO, num ber of simulations = 100. 22 Figure 20: SNR curves of the IRLS, num ber of simulations = 100. 23 Ac kno wledgement W e wish to thank Nasim Bagheri for her assistance with the text editing. App endix A. Conv ergence Analysis A t first, assume that | α 1 | = | α 2 | (linear con vergence): Cases 1 . 1 and 1 . 2, if α 1 = α 2 = ± α ( α > 0) then e N L = ( ± α − 1 )( ± α ) × e 2 2 − e 2 2 ( ± α − 1 + ( ± α ) − 2) e 2 = (1 − 1) ± α − 1 + ( ± α ) − 2 e 2 whic h is equal to zero as exp ected. Case 1.3, if α 1 = − α 2 = α > 0 then e N L = ( α − 1 )( − α ) − 1 α − 1 − α − 2 e 2 = − 2 α − 1 − α − 2 e 2 ⇒ | e N L e 3 | = 2 | α 2 + 2 α − 1 | = 2 | ( α + 1) 2 − 2 | whic h means that for a con v ergent sequence (0 < α < 1), the NL estimation div erges; in other words, we hav e | e N L e 3 | > 1. Case 1.4, if − α 1 = α 2 = α > 0 then e N L = ( − α − 1 )( α ) − 1 − α − 1 + α − 2 e 2 = − 2 − α − 1 + α − 2 e 2 ⇒ | e N L e 3 | = 2 | α 2 − 2 α − 1 | = 2 | ( α − 1) 2 − 2 | whic h also results in the divergence of the NL method in the case of a conv erging sequence (0 < α < 1). In a more general case, it can b e assumed that | α 1 | 6 = | α 2 | . Hence, e 1 and e 3 should b e rewritten in terms of e 2 . In a conv erging algorithm, subsequen t estimations satisfy | e 1 | > | e 2 | > | e 3 | . W e can ass ume that after each iteration, the algorithm b ecomes less capable of reducing the errors whic h results in | α 1 | = α , | α 2 | = (1 + δ ) α ; δ > 0 . (A.1) 24 The latter alongside the con vergence of the algorithm leads to obtaining | α 1 | , | α 2 | < 1 ⇒ α < 1 , δ < 1 α − 1 . (A.2) Hence, w e hav e the follo wing, whic h represent the case of sub-linear conv ergence: Case 2.1, if all of the e i ’s hav e the same sign, in other words, if e 1 × e 2 > 0 and e 3 × e 2 > 0, then e N L = α − 1 × (1 + δ ) α − 1 α − 1 + (1 + δ ) α − 2 e 2 = δ α ( α − 1) 2 + δ α 2 ⇒ | e N L e 3 | = δ (1 + δ )(( α − 1) 2 + δ α 2 ) ≤ 1 where the equalit y holds for δ = 1 α − 1; this implies that the NL metho d does not impro ve the existing estimation, as exp ected from (A.1) and (A.2). Case 2.2, if e 1 × e 2 < 0 and e 3 × e 2 < 0 then e N L = − α − 1 × − (1 + δ ) α − 1 − α − 1 − (1 + δ ) α − 2 e 2 = δ α − ( α + 1) 2 − δ α 2 ⇒ | e N L e 3 | = δ (1 + δ )(( α + 1) 2 + δ α 2 ) < 1 ; ∀ δ > 0 . Case 2.3, if e 1 × e 2 > 0 , e 3 × e 2 < 0 and 0 < δ < 1 α − 1 then e N L = − α − 1 × (1 + δ ) α − 1 α − 1 − (1 + δ ) α − 2 e 2 = − ( δ + 2) α − (( α + 1) 2 − 2 + δ α 2 ) ⇒ | e N L e 3 | = ( δ + 2) α (1 + δ )(( α + 1) 2 − 2 + δ α 2 ) > 1 . Case 2.4, if e 1 × e 2 < 0 and e 3 × e 2 > 0 then e N L = α − 1 × − (1 + δ ) α − 1 − α − 1 + (1 + δ ) α − 2 e 2 = − ( δ + 2) α ( α − 1) 2 − 2 + δ α 2 ⇒ | e N L e 3 | = ( δ + 2) α (1 + δ )(( α − 1) 2 − 2 + δ α 2 ) > 1 where 0 < δ < 1 α − 1. In the case of sup er-linear conv ergent algorithms, it can b e assumed that | α 1 | = α , | α 2 | = (1 − δ ) α ; 0 < δ < 1 (A.3) 25 whic h results in the same pro cedure as for the case of sub-linear conv ergence except for the sign of δ . Case 3.1, if all of the e i ’s ha ve the same sign, then | e N L e 3 | = δ (1 − δ )(( α − 1) 2 − δ α 2 ) < 1 ; for δ < δ 0 (A.4) where δ 0 = α 2 − α +1 − √ 2 α 2 − 2 α +1 α 2 . Hence, the NL estimation diverges for δ 0 < δ < 1. It can be seen that for a cubically con vergen t sequence ( δ = 1 − α ), the NL estimation div erges since δ 0 < 1 − α . Case 3.2, if e 1 × e 2 < 0 and e 3 × e 2 < 0, then | e N L e 3 | = δ (1 − δ )(( α + 1) 2 − δ α 2 ) < 1 ; δ < δ 0 (A.5) where δ 0 = α 2 + α +1 − √ 2 α 2 +2 α +1 α 2 . In this case, sequences with quadratic conv er- gence can b e accelerated if δ 0 > 1 − α . In other words, these sequences can b e accelerated if 0 . 36110 < α . Case 3.3, if e 1 × e 2 > 0 , e 3 × e 2 < 0 and 0 < δ < 1 then | e N L e 3 | = (2 − δ ) α | (1 − δ )(( α + 1) 2 − 2 − δ α 2 ) | > 1 . Case 3.4, if e 1 × e 2 < 0 and e 3 × e 2 > 0 then | e N L e 3 | = ( − δ + 2) α | (1 − δ )(( α − 1) 2 − 2 − δ α 2 ) | > 1 ; 0 < δ < 1 . References References [1] K. Gro c henig, Acceleration of the frame algorithm, IEEE T ransactions on Signal Pro cessing 41 (12) (1993) 3331–3340. [2] Y. Saad, Cheb yshev acceleration tec hniques for solving nonsymmetric eigen v alue problems, Mathematics of Computation 42 (166) (1984) 567– 588. 26 [3] J. No cedal, S. J. W righ t, Conjugate gradi ent metho ds, Numerical optimiza- tion (2006) 101–134. [4] J. D. Blanchard, J. T anner, K. W ei, Cgiht: conjugate gradient iterative hard thresholding for compressed sensing and matrix completion, Informa- tion and Inference: A Journal of the IMA 4 (4) (2015) 289–327. [5] H. G. F eich tinger, K. Gr¨ oc henig, Theory and practice of irregular sampling. [6] F. A. Marv asti, An iterativ e metho d to comp ensate for the interpolation distortion, IEEE T ransactions on Acoustics, Sp eech, and Signal Pro cessing 37 (10) (1989) 1617–1621. [7] M. Ghandi, M. M. J. Y ekta, F. Marv asti, Some nonlinear/adaptiv e metho ds for fast recov ery of the missing samples of signals, Signal Pro cessing 88 (3) (2008) 624–638. [8] A. C. Aitken, Xxv.on b ernoulli’s numerical solution of algebraic equations, Pro ceedings of the Roy al So ciet y of Edin burgh 46 (1927) 289–305. [9] J. Duchi, E. Hazan, Y. Singer, Adaptiv e subgradient metho ds for online learning and sto chastic optimization, Journal of Machine Learning Re- searc h 12 (Jul) (2011) 2121–2159. [10] T. Tieleman, G. Hin ton, Lecture 6.5-rmsprop: Divide the gradien t b y a running av erage of its recent magnitude, COURSERA: Neural netw orks for mac hine learning 4 (2) (2012) 26–31. [11] D. P . Kingma, J. Ba, Adam: A method for sto c hastic optimization, arXiv preprin t [12] R. Tibshirani, Regression shrink age and selection via the lasso, Journal of the Roy al Statistical Society: Series B (Metho dological) 58 (1) (1996) 267–288. 27 [13] T. Hastie, R. Tibshirani, J. F riedman, R. Tibshirani, S. P ow ers, J. Qian, K. Jung, A. Sch uler, N. Shah, Q. Zhao, et al., The elements of statistical learning: Data mining, inference, and prediction. [14] M. Y uan, Y. Lin, Model selection and estimation in regression with grouped v ariables, Journal of the Ro yal Statistical So ciety: Series B (Statistical Metho dology) 68 (1) (2006) 49–67. [15] C. C. P aige, M. A. Saunders, Lsqr: An algorithm for sparse linear equations and sparse least squares, A CM T ransactions on Mathematical Soft ware (TOMS) 8 (1) (1982) 43–71. [16] B. W ahlb erg, S. Boyd, M. Annergren, Y. W ang, An admm algorithm for a class of total v ariation regularized estimation problems, IF A C Pro ceedings V olumes 45 (16) (2012) 83–88. [17] S. Bo yd, N. P arikh, E. Ch u, B. Peleato, J. Ec kstein, et al., Distributed optimization and statistical learning via the alternating direction metho d of multipliers, F oundations and T rends R  in Machine learning 3 (1) (2011) 1–122. [18] A. Amini, F. Marv asti, Reconstruction of multiband signals from non- in vertible uniform and p eriodic non uniform samples using an iterative metho d, in: Pro ceedings of SampT a2007, Thessaloniki, Greece, 2007. [19] F. Marv asti, Nonuniform sampling: theory and practice, Springer Science & Business Media, 2012. [20] D. L. Donoho, M. Elad, V. N. T emlyak ov, Stable recov ery of sparse ov er- complete representations in the presence of noise, IEEE T ransactions on information theory 52 (1) (2006) 6–18. [21] D. L. Donoho, F or most large underdetermined systems of linear equations the minimal ` 1-norm solution is also the sparsest solution, Comm unica- tions on Pure and Applied Mathematics: A Journal Issued b y the Courant Institute of Mathematical Sciences 59 (6) (2006) 797–829. 28 [22] R. Chartrand, W. Yin, Iterativ ely reweigh ted algorithms for compressive sensing, in: 2008 IEEE In ternational Conference on Acoustics, Sp eec h and Signal Pro cessing, IEEE, 2008, pp. 3869–3872. [23] I. Daub ec hies, R. DeV ore, M. F ornasier, C. S. G¨ unt¨ urk, Iterativ ely rew eighted least squares minimization for sparse reco v ery , Comm unications on Pure and Applied Mathematics: A Journal Issued by the Courant In- stitute of Mathematical Sciences 63 (1) (2010) 1–38. [24] H. Mohimani, M. Babaie-Zadeh, C. Jutten, A fast approac h for ov ercom- plete sparse decomp osition based on smoothed ` 0 norm, IEEE T ransactions on Signal Pro cessing 57 (1) (2009) 289–301. [25] A. Eftekhari, M. Babaie-Zadeh, C. Jutten, H. A. Moghaddam, Robust-sl0 for stable sparse represen tation in noisy settings, in: Acoustics, Sp eec h and Signal Pro cessing, 2009. ICASSP 2009. IEEE International Conference on, IEEE, 2009, pp. 3433–3436. [26] T. Blumensath, M. E. Davies, Iterative thresholding for sparse approxima- tions, Journal of F ourier analysis and Applications 14 (5-6) (2008) 629–654. [27] T. Blumensath, M. E. Da vies, Iterativ e hard thresholding for compressed sensing, Applied and computational harmonic analysis 27 (3) (2009) 265– 274. [28] F. A. Marv asti, A unified approac h to zero-crossings and nonuniform sam- pling of single and multi-dimensional signals and systems, Chicago, Ill.: F arohk A. Marv asti, 1987. [29] M. Azghani, F. Marv asti, Iterative methods for random sampling and com- pressed sensing recov ery , in: 10th In ternational Conference on Sampling Theory and Applications (SAMPT A), 2013. [30] M. Azghani, P . Kosmas, F. Marv asti, Microw av e medical imaging based on sparsit y and an iterative metho d with adaptiv e thresholding, IEEE trans- actions on medical imaging 34 (2) (2015) 357–365. 29 [31] F. Zhou, Z. Lu, C. W ang, W. Sun, S.-T. Xia, Q. Liao, Image quality as- sessmen t based on inter-patc h and intra-patc h similarit y , PloS one 10 (3) (2015) e0116312. [32] Z. W ang, A. C. Bovik, H. R. Sheikh, E. P . Simoncelli, Image quality as- sessmen t: from error visibility to structural similarity , IEEE transactions on image pro cessing 13 (4) (2004) 600–612. [33] A. Jav aheri, H. Zayy ani, M. A. Figueiredo, F. Marv asti, Robust sparse reco very in impulsiv e noise via con tinuous mixed norm, IEEE Signal Pro- cessing Letters 25 (8) (2018) 1146–1150. [34] Z. W ang, E. Simoncelli, A. Bovik, et al., Multi-scale structural similarity for image quality assessment, in: ASILOMAR CONFERENCE ON SIGNALS SYSTEMS AND COMPUTERS, V ol. 2, IEEE; 1998, 2003, pp. 1398–1402. [35] X. Zhang, X. F eng, W. W ang, W. Xue, Edge strength similarit y for image qualit y assessment, IEEE Signal pro cessing letters 20 (4) (2013) 319–322. [36] L. Zhang, L. Zhang, X. Mou, D. Zhang, et al., Fsim: a feature similarity index for image qualit y assessmen t, IEEE transactions on Image Pro cessing 20 (8) (2011) 2378–2386. [37] A. Mittal, A. K. Mo orth y , A. C. Bovik, Blind/referenceless image spatial qualit y ev aluator, in: Signals, Systems and Computers (ASILOMAR), 2011 Conference Record of the F ort y Fifth Asilomar Conference on, IEEE, 2011, pp. 723–727. [38] A. Mittal, A. K. Mo orthy , A. C. Bovik, No-reference image quality as- sessmen t in the spatial domain, IEEE T ransactions on Image Pro cessing 21 (12) (2012) 4695–4708. [39] A. Mittal, R. Soundarara jan, A. C. Bovik, Making a” completely blind” image quality analyzer., IEEE Signal Proce ss. Lett. 20 (3) (2013) 209–212. 30