Conjugate gradient acceleration of iteratively re-weighted least squares methods

Introduction

Iteratively Re-weighted Least Squares (IRLS) is a method for solving minimization problems by transforming them into a sequence of easier quadratic problems which are then solved with efficient tools of numerical linear algebra. Contrary to classical Newton methods smoothness of the objective function is not required in general. We refer to the recent paper OBD14 for an updated and rather general view about these methods.

In the context of constructive approximation, an IRLS algorithm appeared for the first time in the doctoral thesis of Lawson in 1961 LW61 in the form of an algorithm for solving uniform approximation problems. It computes a sequence of polynomials that minimize a sequence of weighted $L_{\tau}$–norms. This iterative algorithm is now well-known in classical approximation theory as Lawson’s algorithm. In Cline72 it is proved that this algorithm essentially obeys a linear convergence rate.

In recent years, there has been an explosion of papers on applications and variations on the theme of IRLS, especially in the engineering community of signal processing, and it is by now almost impossible to give a complete account of the developments. (Presently Scholar Google reports more than 3180 papers since 2010 containing the phrase “Iteratively Re-weighted Least Squares” and more than 100 with it in the title since 1970, half of which appeared after 2003.)

2 Contribution of this paper

Besides analyzing the effect of CG in an IRLS for problems of the type (1), we further extend it in Section 4 to a class of problems of the type

for $0<\tau\leqslant 1$, used for sparse recovery in signal processing. In the work xulaiyin ; sergei ; davo15 a convergence analysis of IRLS towards the solution of (2) has been carried out with two limitations:

In xulaiyin the authors do not consider the use of an iterative algorithm to solve the appearing system of linear equations and they do not show the behavior of the algorithm when the measurements $y$ are given with additional noise;

Also in sergei ; davo15 a precise analysis of convergence is missing when iterative methods are used to solve the intermediate sequence of systems of linear equations. Also the non-convex case of $\tau<1$ is not specifically addressed.

Regarding these gaps, we contribute in this work by

giving a proper analysis of the convergence when inaccurate CG solutions are used;

extending the results of convergence in sergei ; davo15 to the case of $0<\tau<1$ by combining our analysis with findings in RamlauZarzer12 ; Zarzer09 ;

performing numerical tests which evaluate possible speedups via the CG method, also taking problems into consideration where measurements may be affected by noise.

Our work on CG accelerated IRLS for (2) does not analytically address rates of convergence because this turned out to be a very technical task.

We illustrate the theoretical results of this paper described above by several numerical experiments. We first show that our versions of IRLS yield significant improvements in terms of computational time and may outperform state of the art first order methods such as Iterative Hard Thresholding (IHT) Blumensath09 and Fast Iterative Soft-Thresholding Algorithm (FISTA) beck09 , especially in high dimensional problems ($N\geqslant 10^{5}$). These results are somehow both surprising and counterintuitive as it is well-known that first order methods should be preferred in higher dimension. However, they can be easily explained by observing that in certain regimes preconditioning in the conjugate gradient method (as we show at the end of Subsection 5.3) turns out to be extremely efficient. This is perhaps not a completely new discovery, as benefits of preconditioning in IRLS have been reported already in minimization problems involving total variation terms VO98 . The second significant outcome of our experiments is that CG-IRLS not only is faster than state of the art first order methods, but also shows higher recovery rates, i.e., requires less measurements for successful sparse recovery. This will be demonstrated with corresponding phase transition diagrams of empirical success rates (Figure 4).

3 Outline of the paper

The paper is organized as follows: In Section 2, we introduce definitions and notation and give a short review on the CG method. Although this brief introduction on CG retraces very well-known facts of the numerical linear algebra literature, it is necessary for us for the sake of a consistent presentation also in terms of notation. We hope that this small detour will help readers to access more easily the technical parts of the paper. In Section 3, we present the IRLS method tailored to problems of the type (1) and its modification including CG for the solution of the quadratic optimizations. We present a detailed analysis of the convergence and rate of convergence. The approach is further extended to problems of type (2) in Section 4, where we also analyze the convergence of the method. We conclude with numerical experiments in Section 5 showing that the modifications to IRLS inspired by our theoretical results make the algorithm extremely efficient, also compared to state of the art first order methods, especially in high dimension.

Definitions, Notation, and Conjugate Gradient method

In this section, we introduce the main terms and notation used in this paper. In addition to this, we shortly review the basics around the Conjugate Gradient method. For a more detailed introduction to conjugate gradient methods, we refer to respective text books, e.g., nowr2006 ; sarisa00 . In order to simplify cross-reading, we use the same notation as in dadefogu10 .

where $\lambda_{\max}(\cdot)$ denotes the largest eigenvalue of a square matrix (compare Definition 5).

for all sets $T\subseteq\{1,\ldots,N\}$ with $\#T\leq{K}$ and all $\eta\in\ker\Phi\backslash\{0\}$. We say in short that $\Phi$ has the $({K},\gamma_{K})$-NSP.

We give an important consequence of the NSP codade09 ; fora13 , (dadefogu10, , Lemma 7.6).

It is well-known that the NSP for $0<\tau\leqslant 1$ can be shown via the restricted isometry property Chartrand08 ; fora13 , but also direct proofs of the NSP are available for certain random matrices giving often better constants and working under weaker assumptions chgulepa12 ; dilera15 ; fora13 ; kara13 ; leme15 . In particular, Gaussian random matrices satisfy the NSP of order $K$ with high probability if $m\geqslant CK\log(K/N)$. Structured random matrices including random partial Fourier and discrete cosine matrices, and partial random circulant matrices – both important in applications – satisfy the RIP and hence, the NSP with high probability provided that $m\geqslant CK\log^{4}(N)$ cata06 ; fora13 ; krmera14 ; ra10 ; ruve08 . Note that for these types of structured matrices, fast matrix vector multiplication routines are available.

We denote with $\Lambda(A)$ the set of eigenvalues of a square matrix A. Respectively, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ are the smallest and largest eigenvalues. We define by $\sigma_{\min}(A)$ and $\sigma_{\max}(A)$ the smallest and largest singular value of a rectangular matrix $A$.

The following theorem establishes the convergence and the convergence rate of CG.

Let the matrix $A$ be Hermitian and positive definite. The Algorithm CG converges to the solution of the system $Ax=y$ after at most $N$ steps. Moreover, the error $x^{i}-x$ is such that

where $\kappa_{A}=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}$ is the condition number of the matrix $A$ and $\sigma_{\max}(A)$ (resp. $\sigma_{\min}(A)$) is the largest (resp. smallest) singular value of $A$.

Theorem 2.1 is slightly modified with respect to the formulation in sarisa00 . There, the matrix $A$ is considered to be symmetric instead of being Hermitian. However, in the complex case, the proof can be performed similarly by replacing the transpose by the conjugate transpose.

2 Modified conjugate gradient method (MCG)

In Section 3, we are interested in a vector which solves the weighted least-squares problem

where $D\mathrel{\mathop{:}}=\operatorname*{diag}{[w_{i}^{-1}]}_{i=1}^{N}$. Hence, in order to determine $\hat{x}$, we first solve the system

and then we compute $\hat{x}=D\Phi^{*}\theta$. Notice that the system (7) has the general form

for all $i\geqslant 0$, where $c_{TT^{*}}{=\frac{\sqrt{\kappa(TT^{*})}-1}{\sqrt{\kappa(TT^{*})}+1}=\frac{\sigma_{\max}(T)-\sigma_{\min}(T)}{\sigma_{\max}(T)+\sigma_{\min}(T)}}$ is defined as in Theorem 2.1, and $\bar{x}^{0}=T^{*}\theta^{0}$ is the initial vector. Moreover, by setting $D\mathrel{\mathop{:}}=\operatorname*{diag}{[w_{i}^{-1}]}_{i=1}^{N}$, and $\hat{x}^{i}=D^{\frac{1}{2}}\bar{x}^{i}$ as well as $\hat{x}=D^{\frac{1}{2}}\bar{x}$, we obtain

for $\theta$ as given in (8). By the identity

In this section, we start with a detailed introduction of the IRLS algorithm and its modified version that uses CG for the solution of the successive quadratic optimization problems. Afterwards, we present two results providing the convergence and the rate of convergence of the modified algorithm. As crucial feature, we give bounds on the accuracies of the (inexact) CG solutions of the intermediate least squares problems which ensure convergence of the overall IRLS methods. In particular, these tolerances must depend on the current iteration and should tend to zero with increasing iteration count. In fact, without this condition, one may observe divergence of the method. The proofs of the theorems are developed into several lemmas.

The following functional turns out to be a crucial tool for the analysis of the IRLS algorithm and its modified variant.

The convergence of IRLS is by now well-established and we refer to dadefogu10 and (fora13, , Section 15.3) for details, which we in part extend in our analysis in Section 3.3.

and thus, since $\Phi$ has full rank and $\Phi D\Phi^{*}$ is invertible, we conclude

As a consequence, we see that at step 2 of Algorithm IRLS the minimizer of the least squares problem is explicitly given by the equation

where we introduced the $N\times N$ diagonal matrix

Furthermore, the new weight vector in step 4 of Algorithm IRLS is explicitly given by

Taking into consideration that $w_{j}>0$, this formula can be derived from the first order optimality condition ${\partial}\mathcal{J}_{\tau}(x^{n+1},w,{\varepsilon}^{n+1})/{\partial}w=0$.

2 The algorithm CG-IRLS

In contrast to Algorithm IRLS, the value $\beta$ in step 3 is introduced to obtain flexibility in tuning the performance of the algorithm. While we prove in Theorem 3.1 convergence for any positive value of $\beta$, Theorem 3.1(iii) below guarantees instance optimality only for $\beta<\left(\frac{1-\gamma}{1+\gamma}\frac{K+1-k}{N}\right)^{\frac{1}{\tau}}$ in the case that $\lim\limits_{n\rightarrow\infty}{\varepsilon}^{n}\neq 0$. Nevertheless in practice, choices of $\beta$ which do not necessarily fulfill this condition may work very well. Section 5, investigates good choices of $\beta$ numerically.

The vector $\hat{x}=\hat{x}^{n+1}$ is not known a priori;

The computation of the condition number $c_{TT^{*}}$ is possible, but it requires the computation of eigenvalues with additional computational cost which we prefer to avoid.

We use $\theta^{n+1,i}=(\Phi D_{n}\Phi^{*})^{-1}(y-\rho^{n+1,i})$ from step 5 of MCG to obtain

The last inequality above results from $\lambda_{\min}\left(\Phi D_{n}\Phi^{*}\right)=\sigma_{\min}^{2}\left(\Phi D_{n}^{\frac{1}{2}}\right)$ and

In inequality (15), the computation of $\sigma_{\min}\left(\Phi\right)$ and $\|\Phi\|$ is necessary. The computation of these constants might be demanding, but has to be performed only once before the algorithm starts. Furthermore, in practice it is sufficient to compute approximations of these values and therefore these operations are not critical for the computation time of the algorithm.

3 Convergence results

After introducing Algorithm CG-IRLS, we state below the two main results of this section. Theorem 3.1 shows the convergence of the algorithm to a limit point that obeys certain error guarantees with respect to the solution of (1). Below $K$ denotes the index used in the ${\varepsilon}$-update rule, i.e., step 3) of Algorithm CG-IRLS.

Let $0<\tau\leqslant 1$. Assume $K$ is such that $\Phi$ satisfies the Null Space Property (5) of order $K$, with $\gamma<1$. If $\textnormal{tol}_{n+1}$ in Algorithm CG-IRLS is chosen such that

If ${\varepsilon}>0$, then for each $\bar{x}\in\mathcal{Z}_{\tau}(y)\neq\emptyset$, we have $\left\langle\bar{x},\eta\right\rangle_{\hat{w}(\bar{x},{\varepsilon},\tau)}=0$ for all $\eta\in\mathcal{N}_{\Phi}$, where $\hat{w}(\bar{x},{\varepsilon},\tau)=\left[\left||\bar{x}_{i}|^{2}+{\varepsilon}^{2}\right|^{-\frac{2-\tau}{2}}\right]_{i=1}^{N}$. Moreover, in the case of $\tau=1$, $\bar{x}$ is the single element of $\mathcal{Z}_{\tau}(y)$ and $\bar{x}=x^{{\varepsilon},1}\mathrel{\mathop{:}}=\operatorname*{arg\,min\,}\limits_{x\in\mathcal{F}_{\Phi}(y)}\sum\limits_{j=1}^{N}|x_{j}^{2}+{\varepsilon}^{2}|^{\frac{1}{2}}$ (compare (42)).

Denote by $\mathcal{X}_{{\varepsilon},\tau}(y)$ the set of global minimizers of $f_{{\varepsilon},\tau}(x)\mathrel{\mathop{:}}=\sum\limits_{j=1}^{N}|x_{j}^{2}+{\varepsilon}^{2}|^{\frac{\tau}{2}}$ on $\mathcal{F}_{\Phi}(y)$. If ${\varepsilon}>0$ and $\bar{x}\in\mathcal{Z}_{\tau}(y)\cap\mathcal{X}_{{\varepsilon},\tau}(y)$, then for each $x\in\mathcal{F}_{\Phi}(y)$ and any $\beta<\left(\frac{1-\gamma}{1+\gamma}\frac{K+1-k}{N}\right)^{\frac{1}{\tau}}$, we have

Knowing that the algorithm converges and leads to an adequate solution, one is also interested in how fast one approaches this solution. Theorem 3.2 states that a linear rate of convergence can be established in the case of $\tau=1$. In the case of $0<\tau<1$ this rate is even asymptotically super-linear.

Assume $\Phi$ satisfies the NSP of order $K$ with constant $\gamma$ such that $0<\gamma<1-\frac{2}{K+2}$, and that $\mathcal{F}_{\Phi}(y)$ contains a $k$-sparse vector $x^{*}$. Define $\Lambda\mathrel{\mathop{:}}=\operatorname*{supp}(x^{*})$. Suppose that $k<K-\frac{2\gamma}{1-\gamma}$ and $0<\nu<1$ are such that

If $a_{n+1}$ and $tol_{n+1}$ are chosen as in Theorem 3.1 with the additional bound

Note that the second bound in (23), which implies (25), is only of theoretical nature. Since the value of $E_{n}$ is unknown it cannot be computed in an implementation. However, heuristic choices of $\textnormal{tol}_{n+1}$ may fulfill this bound. Thus, in practice one can only guarantee the “asymptotic” (super-)linear convergence (24).

In the remainder of this section we aim to prove both results by means of some technical lemmas which are reported in Section 3.3.1 and Section 3.3.2.

One important issue in the investigation of the dynamics of Algorithm CG-IRLS is the relationship between the weighted norm of an iterate and the weighted norm of its predecessor. In the following lemma, we present some helpful estimates.

hold for all $n\geqslant 1$, where $W_{n}\mathrel{\mathop{:}}=\left\|D_{n}^{-\frac{1}{2}}D_{n-1}^{\frac{1}{2}}\right\|$.

where the last inequality is due to (26).

The functional $\mathcal{J}_{\tau}(x,w,{\varepsilon})$ obeys the following monotonicity property.

The first inequality follows from the minimization property of ${w}^{n+1}$. The second inequality follows from ${\varepsilon}^{{n+1}}\leq{\varepsilon}^{{n}}$.

The following lemma describes how the difference of the functional, evaluated in the exact and the approximated solution can be controlled by a positive scalar $a_{n+1}$ and an appropriately chosen tolerance $\textnormal{tol}_{n+1}$.

where $\bar{W}_{n+1}$ was defined in (18). By choosing $\textnormal{tol}_{n+1}$ as in (16), we obtain

where we have used the Cauchy-Schwarz inequality in the first inequality, (26) and (27) in the fifth inequality, (32) in the third inequality, the definition of $c_{n}$ in (17), and the Assumption (16) on $\textnormal{tol}_{n+1}$ in the last inequality.

Since $1\leqslant\bar{W}_{n+1}$, we obtain (30) by

with the same arguments as above. Lemma 5 yields

where the first inequality follows from (29), the second and third by (28), and the last by (30).

Identity (33) follows by insertion of the definition of ${w}^{n+1}$ in step 4 of Algorithm CG-IRLS.

By the minimizing property of $\hat{x}^{n+1}$ and the fact that $\hat{x}^{n}\in\mathcal{F}_{\Phi}(y)$, we have

and thus, together with (31), it follows that

Inequality (34) then follows from (29) and

Consequently, the bound (35) follows from

Inequality (36) is a direct consequence of (35).

Notice that (34) states the boundedness of the iterates. The lower bound (35) on the weights ${w}^{n}$ will become useful in the proof of Lemma 8.

where C is the constant of Lemma 7 and $\hat{x}^{n}=\operatorname*{arg\,min\,}\limits_{x\in\mathcal{F}_{\Phi}(y)}\mathcal{J}_{\tau}\left(x,{w}^{n-1},{\varepsilon}^{n-1}\right)$. As a consequence we have

Here we used the fact that $\hat{x}^{n}-\hat{x}^{n+1}\in\mathcal{N}_{\Phi}$ and therefore, $\left\langle\hat{x}^{n+1},\hat{x}^{n}-\hat{x}^{n+1}\right\rangle=0$ and in the last step we applied the bound (36). Summing these inequalities over $n\geq 1$, we arrive at

Letting $N\rightarrow\infty$ yields the desired result.

The following lemma will play a major role in our proof of convergence since it shows that not only (38) holds but that also the difference between successive iterates becomes arbitrarily small.

By (36) of Lemma 7 and the condition (16) on $\textnormal{tol}_{n}$, we have

Together with Lemma 8 we can prove our statement:

where the first and last term vanish because of (40) and the other term due to (38).

In this section, we introduce an auxiliary functional which is useful for the proof of convergence. From the monotonicity of ${\varepsilon}_{n}$, we know that ${\varepsilon}=\lim\limits_{n\rightarrow\infty}{\varepsilon}_{n}$ exists and is nonnegative. We introduce the functional

In the case of $0<\tau<1$, we denote by $\mathcal{X}_{{\varepsilon},\tau}(y)$ the set of global minimizers of $f_{{\varepsilon},\tau}$ on $\mathcal{F}_{\Phi}(y)$. For both cases, the minimizers are characterized by the following lemma.

Let ${\varepsilon}>0$ and $x\in\mathcal{F}_{\Phi}(y)$. If $x=x^{{\varepsilon},1}$ or $x\in\mathcal{X}_{{\varepsilon},\tau}(y)$, then $\left\langle x,\eta\right\rangle_{\hat{w}(x,{\varepsilon},\tau)}=0$ for all $\eta\in\mathcal{N}_{\Phi}$, where $\hat{w}(x,{\varepsilon},\tau)=\left[\left||x_{i}|^{2}+{\varepsilon}^{2}\right|^{-\frac{2-\tau}{2}}\right]_{i=1}^{N}$. In the case of $\tau=1$ also the converse is true.

The proof is an adaptation of (dadefogu10, , Lemma 5.2, Section 7) and is presented for the sake of completeness in Appendix A.

3.3 Proof of convergence

By the results of the previous section, we are able now to prove the convergence of Algorithm CG-IRLS. The proof is inspired by the ones of (dadefogu10, , Theorem 5.3, Theorem 7.7), see also (fora13, , Chapter 15.3), which we adapted to our case.

To prove (iii), assume that $\bar{x}\in\mathcal{Z}_{\tau}(y)\cap\mathcal{X}_{{\varepsilon},\tau}(y)$, and follow the proof of (dadefogu10, , Theorem 5.3, and 7.7) to conclude.

3.4 Proof of rate of convergence

We apply the characterization (12) with $w={w}^{n}$, $\hat{x}=\hat{x}^{n+1}=x^{*}+\hat{\eta}^{n+1}$, and $\eta=\hat{x}^{n+1}-x^{*}=\hat{\eta}^{n+1}$, which gives

Rearranging the terms and using the fact that $x^{*}$ is supported on $\Lambda$, we obtain

By assumption there exists $n_{0}$ such that $E_{n_{0}}\leqslant R^{*}$. We prove (24), and $E_{n}\leqslant R^{*}\Rightarrow E_{n+1}\leqslant R^{*}$ to obtain the validity for all $n\geqslant n_{0}$. Assuming $E_{n}\leqslant R^{*}$, we have for all $j\in\Lambda$,

Hence, (43) combined with (44) and the NSP leads to

Combining (dadefogu10, , Proposition 7.4) with the above estimate yields

Note that this is also valid if $\hat{\eta}_{\Lambda^{c}}^{n+1}=0$ since then the left-hand side is zero and the right-hand side non-negative. We furthermore obtain

In addition to this, we know by (dadefogu10, , Lemma 4.1, 7.5), that

where we used the triangle inequality in the first inequality, (36) in the third inequality, and $C$ is the constant from Lemma 7.

Equation (25) then follows by condition (23). By means of (20), we obtain

and therefore the linear convergence for $\tau=1$, and the super-linear convergence for $\tau<1$ as soon as $n\geqslant n_{0}$.

Lai, Xu, and Yin in xulaiyin and Daubechies and Voronin in sergei ; davo15 showed independently that computing the optimizer of the problem (49) can be approached by an alternating minimization of the functional $J_{\tau,\lambda}$ with respect to $x$, $w$, and $\varepsilon$. The difference between these two works is the definition of the update rule for $\varepsilon$. Here, we chose the rule in step 4 of Algorithm 5 proposed by Daubechies and Voronin because it allows us to show that the algorithm converges to a minimizer of (49) for $\tau=1$ and to critical points of (49) for $\tau<1$ (more precise statements will be given below). However, we were not able to prove similar statements for the rule of Lai, Xu, and Yin. It only allows to show the convergence of the algorithm to a critical point of the smoothed functional

We approach the first step of the algorithm by computing a critical point of $J_{\tau,\lambda}(\cdot,w,\varepsilon)$ via the first order optimality condition

We denote the solution of this system by $x^{n+1}$. The new weight $w^{n+1}$ is obtained in step 3 and can be expressed componentwise by

Similarly to the previous section we propose the combination of Algorithm 5 with the CG method. CG is used to calculate an approximation of the solution of the linear system (52) in line 3 of the algorithm. After including the CG method, the modified algorithm which we shall consider is Algorithm CG-IRLS-$\lambda$.

Notice that the matrix $\Phi^{*}\Phi$ is positive semi-definite and $\lambda\tau D_{n}^{-1}=\lambda\tau\operatorname*{diag}\left[w_{j}^{n}\right]_{j=1}^{N}$ is positive definite. Therefore, $A_{n}$ is positive definite and invertible, and furthermore

The second factor of (55) is estimated by

where we used (54) in the inequality. Thus, we obtain

In the remainder of this section, we clarify how to choose the tolerance $\textnormal{tol}_{n+1}$, and establish a convergence result of the algorithm. In the case of $\tau=1$, the problem (49) is the minimization of the well-known LASSO functional. It is convex, and the optimality conditions can be stated in terms of subdifferential inclusions. We are able to show that at least a subsequence of the algorithm is converging to a solution of (49). If $0<\tau<1$, the problem is non-convex and non-smooth. Necessary first order optimality conditions for a global minimizer of this functional were derived in (BrediesLorenz14, , Proposition 3.14), and (ItoKunisch14, , Theorem 2.2). In our case, we are able to show that the non-zero components of the limits of the algorithm fulfill the respective conditions. However, as soon as the algorithm is producing zeros in some components of the limit, so far, we were not able to verify the conditions mentioned above. On this account, we pursue a different strategy, which originates from Zarzer09 . We do not directly show that the algorithm computes a solution of problem (49). Instead we show that a subsequence of the algorithm is at least computing a point $x^{\dagger}$, whose transformation $\breve{x}^{\dagger}=\mathcal{N}_{\upsilon/\tau}^{-1}(x^{\dagger})$ is a critical point of the new functional

is a continuous bijective mapping and $1<\upsilon\leqslant 2$. It was shown in Zarzer09 ; RamlauZarzer12 that assuming $\breve{x}^{\dagger}$ is a global minimizer of $\breve{F}_{\upsilon,\lambda}(x)$ implies that $x^{\dagger}$ is a global minimizer of $F_{\tau,\lambda}$, i.e., a solution of problem (49). Furthermore, it was also shown that this result can be partially extended to local minimizers. We comment on this issue in Remark 7. These considerations allow us to state the main convergence result.

As we argued in Remark 3, a possible relaxation of the tolerance bound (16) is allowed to further boost the convergence, the same applies to the bound (59).

In the case $0<\tau<1$, the theorem includes the possibility that there may exist several converging subsequences with different limits. Potentially only one of these limits may have the nice property that its transformation is a critical point. In the proof of the theorem, which follows further below, an appropriate subsequence is constructed. Actually this construction leads to the following hint, how to practically choose the subsequence: Take a converging subsequence $x_{n_{l}}$ for which the $n_{l}$ satisfy equation (85).

It will be important below that a minimizer $x^{\sharp}$ of $F_{1,\lambda}(x)$ is characterized by the conditions

Note that in the (less important) case $x^{\lambda}=0$, our theorem does not give a conclusion about $x^{\lambda}$ being a minimizer of $F_{1,\lambda}(x)$.

The result of Theorem 4.1 for $0<\tau<1$ can be partially extended towards local minimizers. For the sake of completeness we sketch the argument from RamlauZarzer12 . Assume that $\breve{x}^{\lambda}$ is a local minimizer. Then there is a neighborhood $U_{\epsilon}(\breve{x}^{\lambda})$ with $\epsilon>0$ such that for all $x^{\prime}\in U_{\epsilon}(\breve{x}^{\lambda})$:

By continuity of $\mathcal{N}_{\upsilon/\tau}$ there exists an $\hat{\epsilon}>0$ such that the neighborhood $U_{\hat{\epsilon}}(x^{\lambda})\subset\mathcal{N}_{\upsilon/\tau}(U_{\epsilon}(\breve{x}^{\lambda}))$. Thus, for all $x\in U_{\hat{\epsilon}}(x^{\lambda})$, we have $x^{\prime}=\mathcal{N}^{-1}_{\upsilon/\tau}(x)\in U_{\epsilon}(\breve{x}^{\lambda})$, and obtain

For the proof of Theorem 4.1, we proceed similarly to Section 3, by first presenting a sequence of auxiliary lemmas on properties of the functional $J_{\tau,\lambda}$ and the dynamics of Algorithm CG-IRLS-$\lambda$.

Since the functional is composed of positive summands, its definition and (66) imply

since $\hat{x}^{n+1}$ is the exact minimizer. From (69) we obtain

Since (69) holds in addition to (65) and (66), we conclude, also for the exact solution $\hat{x}^{n+1}$, the bound

Additionally using (71), we are able to estimate the second summand of (70) by

where we used the Cauchy-Schwarz inequality in the second inequality, the triangle inequality in the third inequality, and the bounds in (67) and (71) in the last inequality.

The following pivotal result of this section allows us to control the difference between the exact and approximate solution of the linear system in line 3 of Algorithm CG-IRLS-$\lambda$.

For a given positive number $a_{n+1}$ and a choice of the accuracy $\textnormal{tol}_{n+1}$ satisfying (59), the functional $J_{\tau,\lambda}$ fulfills the two monotonicity properties

where $C_{w^{n}}$ was defined in (60), we can estimate

where we used (73) in the second inequality, Cauchy-Schwarz in the third inequality, and (68), (67), and (72) in the sixth inequality. Thus we obtain (74). To show (75), we use (68) in the second to last inequality, condition (59) in the last inequality and the fact that $\hat{x}^{n+1}=\underset{x}{\operatorname*{arg\,min\,}}J_{\tau,\lambda}(x,w^{n},\varepsilon^{n})$ (and thus fulfilling (51)) in the second identity below:

Besides Lemma 12 there are two more helpful properties of the functional. First, the identity

where the estimate (68) is used in the first inequality.

2 Proof of convergence

We use the properties of $J$, which we derived in the previous subsection. First, we show (82):

We used (81) in the first inequality, (75) in the second inequality, (63) and (64) in the third inequality, (74) in the fourth inequality and a telescoping sum in the identity. Letting $M\rightarrow\infty$ we obtain

Second, we show (83). From line 1 and 3 of (76) we know that

Since the second summand is positive, we conclude

and thus taking limits on both sides we get

The following lemma provides a lower bound for the $\varepsilon^{n}$, which is used to show a contradiction in the proof of Theorem 4.1. Recall that \phi\in\big{(}0,\frac{1}{4-\tau}\big{)} is the parameter appearing in the update rule for $\varepsilon$ in step 4 of both the algorithms CG-IRLS-$\lambda$ and IRLS-$\lambda$.

Following exactly the steps of the proof of (sergei, , Lemma 4.5.6.) yields the assertion. Observe that all of these steps are also valid for $0<\tau<1$, although in (sergei, , Lemma 4.5.6) the author restricted it to the case $\tau\geqslant 1$.

The observation in the previous proof that $(\varepsilon^{n})$ converges to will be again important below.

We are now prepared for the proof of Theorem (4.1).

Consider the case $\tau=1$ and $x^{\lambda}\neq 0$. We first show that

It follows from equation (51) and the boundedness of the residual (71) that the sequence $(\hat{x}^{n_{k}+1}w_{j}^{n_{k}})_{n_{k}}$ is bounded, i.e.,

Therefore, there exists a converging subsequence, for simplicity again denoted by $(\hat{x}^{n_{k}+1}w_{j}^{n_{k}})_{n_{k}}$. To show the identity in (86), we estimate

In order to show condition (62) for $j$ such that $x^{\lambda}_{j}=0$, we follow the main idea in the proof of Lemma 4.5.9. in sergei . Assume

which is a contradiction, and thus the assumption (87) is false. By means of this result and again a continuity argument, we show condition (62) by

Consider the case $0<\tau<1$. The transformation $\mathcal{N}_{\zeta}(x)$ defined in (58) is continuous and bijective. Thus, $\breve{x}^{\lambda}:=\mathcal{N}_{\upsilon/\tau}^{-1}(x^{\lambda})$ is well-defined, and $x^{\lambda}_{j}=0$ if and only if $\breve{x}^{\lambda}_{j}=0$. At a critical point of the differentiable functional $\breve{F}_{\tau,\lambda}$, its first derivative has to vanish which is equivalent to the conditions

We replace $x^{\lambda}=\mathcal{N}_{\upsilon/\tau}(\breve{x}^{\lambda})$ and obtain

We multiply this identity by $\frac{\upsilon}{\tau}|x_{j}|^{\frac{\upsilon-\tau}{\tau}}$ and obtain (92).

If $\breve{x}^{\lambda}$ is also a global minimizer of $\breve{F}_{\upsilon,\lambda}$, then $x^{\lambda}$ is a global minimizer of $F_{\tau,\lambda}$. This is due the equivalence of the two problems which was shown in (RamlauZarzer12, , Proposition 2.4) based on the continuity and bijectivity of the mapping $\mathcal{N}_{\upsilon/\tau}$ (Zarzer09, , Proposition 3.4).

Numerical Results

We illustrate the theoretical results of this paper by several numerical experiments. We first show that our modified versions of IRLS yield significant improvements in terms of computational time and often outperform the state of the art methods Iterative Hard Thresholding (IHT) Blumensath09 and Fast Iterative Soft-Thresholding Algorithm (FISTA) beck09 .

Before going into the detailed presentation of the numerical tests, we raise two plain numerical disclaimers concerning the numerical stability of CG-IRLS and CG-IRLS-$\lambda$:

The first issue concerns IRLS methods in general: The case where $\varepsilon^{n}\rightarrow 0$, i.e., $x_{j}^{n}\rightarrow 0$, for some $j\in\{1,\ldots,N\}$ and $n\rightarrow\infty$, is very likely since our goal is the computation of sparse vectors. In this case $w^{n}_{j}$ will for some $n$ become too large to be properly represented by a computer. Thus, in practice, we have to provide a lower bound for $\varepsilon$ by some $\varepsilon^{\textrm{min}}>0$. Imposing such a limit has the theoretical disadvantage that in general the algorithms are only calculating an approximation of the respective problems (1) and (49). Therefore, to obtain a “sufficiently good” approximation, one has to choose $\varepsilon^{\textrm{min}}$ sufficiently small. This raises yet another numerical issue: If we choose, e.g., $\varepsilon^{\textrm{min}}=1\text{\sc{e}-}8$ and assume that also $x_{j}^{n}\ll 1$, then $w^{n}_{j}$ is of the order $1\text{\sc{e}+}8$. Compared to the entries of the matrix $\Phi$, which are of the order $1$, any multiplication or addition by such a value will cause serious numerical errors. In this context we cannot expect that the IRLS method reaches high accuracy, and saturation effects of the error are likely to occur before machine precision.

In the following, we start with a description of the general test settings, which will be common for both Algorithms CG-IRLS and CG-IRLS-$\lambda$. Afterwards we independently analyze the speed of both methods and compare them with state of the art algorithms, namely IHT and FISTA. We respectively start with a single trial, followed by a speed-test on a variety of problems. We will also compare the performance of both CG-IRLS and CG-IRLS-$\lambda$ for the noiseless case which leads to surprising results.

All tests are performed with MATLAB version R2014a. For the sake of faster tests (in some cases experiments run for several days) and simplicity, we restrict ourselves to experiments with models defined by real numbers although everything can be similarly done over the complex field. To exploit the advantage of fast matrix-vector multiplications and to allow high dimensional tests, we use randomly sampled partial discrete cosine transformation matrices $\Phi$. We perform tests in three different dimensional settings (later we will extend them to higher dimension) and choose different values $N$ of the dimension of the signal, the amount $m$ of measurements, the respective sparsity $k$ of the synthesized solutions, and the index $K$ in Algorithm (CG-)IRLS:

For each of these settings, we draw at random a set of 100 synthetic problems on which a speed-test is performed. For each synthetic problem the support $\Lambda$ is determined by the first $k$ entries of a random permutation of the numbers $1,\ldots,N$. Then we draw the sparse vector $x^{*}$ at random with entries $x^{*}_{i}\sim\mathcal{N}(0,1)$ for $i\in\Lambda$ and $x^{*}_{\Lambda^{c}}=0$, and a randomly row sampled normalized discrete cosine matrix $\Phi$, where the full non-normalized discrete cosine matrix is given by

In practice, we set the MSNR first and choose the noise level $\sigma=\frac{\sqrt{k}}{\text{MSNR}\sqrt{m}}$. If $\text{MSNR}=\infty$, the problem is noiseless, i.e., $e=0$.

2 Algorithm CG-IRLS

To get an immediate impression about the general behavior of CG-IRLS, we compare its performance in terms of accuracy and speed to IRLS, where the intermediate linear systems are solved exactly via Gaussian elimination (i.e., by the standard MATLAB backslash operator). We choose IHT as a first order state of the art benchmark, to get a fair comparison with another method which can exploit fast matrix-vector multiplications.

For this test, we set the parameter $\beta$ in the $\varepsilon$-update rule to 2. We comment on the choice of this particular parameter in a dedicated paragraph below.

As we have shown by a single trial in the previous paragraph, CG-IRLS as it is presented in Section 3.2 is not able to outperform IHT. Therefore, we introduce the following practical modifications to the algorithm:

We introduce the parameter maxiter_cg, which defines the maximal number of inner CG iterations. Thus, the inner loop of the algorithm stops as soon as maxiter_cg iterations were performed, even if the theoretical tolerance $\textnormal{tol}_{n}$ is not reached yet.

The left plot of Figure 1 reveals that in the beginning CG-IRLS reduces the error more slowly than IHT, and it gets faster after it reached a certain ball around the solution. Therefore, we use IHT as a warm up for CG-IRLS, in the sense that we apply a number start_iht of IHT iterations to compute a proper starting vector for CG-IRLS.

We call CG-IRLSm the algorithm with modifications (i) and (ii), and IHT+CG-IRLSm the algorithm with modifications (i), (ii), and (iii). We set $\texttt{maxiter\_cg}=\lfloor m/12\rfloor$, $\texttt{start\_iht}=150$, and we set $\beta$ to 0.5. If these algorithms are executed on the same trial as in the previous paragraph, we obtain the result which is shown on the right plot in Figure 1. For this trial, the modified algorithms show a significantly reduced computational time with respect to the unmodified version and they now converge faster than IHT. However, the introduction of the practical modifications (i)–(iii) does not necessarily comply anymore with the assumptions of Theorem 3.1. Therefore, we do not have rigorous convergence and recovery guarantees anymore and recovery might potentially fail more often. In the next paragraph, we empirically investigate the failure rate and explore the performance of the different methods on a sufficiently large test set.

In order to investigate the influence of the tolerance $\textnormal{tol}_{n}$ and the number of (inner) iterations of the MCG procedure performed along the IRLS iterations, we plot both quantities in Figure 2 for CG-IRLS, CG-IRLSm, and IHT+CG-IRLSm. Obviously the tolerance quickly decreases in any method. For IHT+CG-IRLSm also the number of MCG iterations decreases (until the method becomes unstable), while for the other two methods the number of MCG iterations first increases and then decreases again until the methods also become unstable. In CG-IRLSm, the number of MCG iterations is bounded by $\texttt{maxiter\_cg}=\lfloor m/12\rfloor$. This more economical behavior only slightly influences the approximation of the MCG solutions and leads to reduced computational time.

Another natural modification to CG-IRLS consists in the introduction of a preconditioner to compensate for the deterioriation of the condition number of $\Phi D_{n}\Phi^{*}$ as soon as $\varepsilon^{n}$ becomes too small (when $w^{n}$ becomes very large). The matrix $\Phi\Phi^{*}$ is very well conditioned, while the matrix $\Phi D_{n}\Phi^{*}$ “sandwiching” $D_{n}$ becomes more ill-conditioned as $n$ gets larger, and, unfortunately, it is hard to identify additional “sandwiching” preconditioners $P_{n}$ such that the matrix $P_{n}\Phi D_{n}\Phi^{*}P_{n}^{*}$ is suitably well-conditioned. In the numerical experiments standard preconditioners failed to yield any significant improvement in terms of convergence speed. Hence, we refrained from introducing further preconditioners. Instead, as we will show at the end of Subsection 5.3, a standard (Jacobi) preconditioning of the matrix

where the source of singularity is added to the product $\Phi^{*}\Phi$, leads to a dramatic improvement of computational speed.

In each setting we check for each trial which methods succeeds or fails. If all methods succeed, we compare the computational time, determine the fastest method, and count the computational time of each method for the respective mean computational time. The results are shown in Figure 3. By analyzing the diagrams, we are able to distill the following observations:

Especially in Setting A and B, CG-IRLSm and IHT+CG-IRLSm are better or comparable to IHT in terms of mean computational time and provide in most cases the fastest method. CG-IRLS performs much worse. The failure rate of all the methods is negligible here.

The gap in the computational time between all methods becomes larger when $N$ is larger.

With increasing dimension of the problem, the advantage of using the modified CG-IRLS methods subsides, in particular in Setting C.

In the literature Chartrand07 ; Chartrand08 ; ChartrandYin08 ; dadefogu10 superlinear convergence is reported for $\tau<1$, and perhaps one of the most surprising outcomes is that the best results for all CG-IRLS methods are instead obtained for $\tau=1$. This can probably be explained by observing that superlinear convergence kicks in only in a rather small ball around the solution and hence does not necessarily improve the actual computation time!

Not only the computational performance, but also the failure rate of the CG-IRLS based methods increases with decreasing $\tau$. However, as expected, CG-IRLS succeeds in the convex case of $\tau=1$. The failure of CG-IRLS for $\tau<1$ can probably be attributed to non-convexity.

We conclude that CG-IRLSm and IHT+CG-IRLSm perform well for $\tau=1$ and for the problem dimension $N$ within the range of 1000 – 10000. They are even able to outperform IHT. However, by extrapolation of the numerical results IHT is expected to be faster for $N>10000$. (This is in compliance with the general folklore that first order methods should be preferred for higher dimension. However, as we will see in Subsection 5.3, a proper preconditioning of CG-IRLS-$\lambda$ will win over IHT for dimensions $N\geqslant 10^{5}$!) As soon as $N<1000$, direct methods such as Gaussian elimination are faster than CG, and thus, one should use standard IRLS with $\tau<1$.

The numerical tests in the previous paragraph were preceded by a careful and systematic investigation of the tuning of the parameters $\beta$, maxiter_cg, and start_iht. While we fixed start_iht to 100, 150, and 200 for Setting A, B, and C respectively to produce a good starting value, we tried $\beta\in\{1/N,0.01,0.1,0.5,0.75,1,1.5,2,5,10\}$, and $\texttt{maxiter\_cg}\in\{\lfloor m/8\rfloor,\lfloor m/12\rfloor,\lfloor m/16\rfloor\}$ for each setting. The results of this parameter sensitivity study can be summarized as follows:

The best computational time is obtained for $\beta\sim 1$. In particular the computational time is not depending substantially on $\beta$ in this order of magnitude. More precisely, for CG-IRLS the choice of $\beta=0.5$ and for (IHT+)CG-IRLSm the choice of $\beta=2$ works best.

The choice of maxiter_cg very much determines the tradeoff between failure and speed of the method. The value $\lfloor m/12\rfloor$ seems to be the best compromise. For a smaller value the failure rate becomes too high, for a larger value the method is too slow.

Besides the empirical analysis of the speed of convergence, we also investigate the robustness of CG-IRLS with respect to the achievable sparsity level for exact recovery of $x^{*}$. Therefore, we fix $N=2000$ and we compute a phase transition diagram for IHT and CG-IRLS on a regular Cartesian $50\times 40$ grid, where one axis represents $m/N$ and the other represents $k/m$. For each grid point we plot the empirical success recovery rate, which is numerically realized by running both algorithms on 20 random trials. CG-IRLS or IHT is successful if it is able to compute a solution with a relative error of less than $1\text{\sc{e}-}4$ within 20 or 500 (outer) iterations respectively. Since we aim at simulating a setting in which the sparsity $k$ is not known exactly, we set the parameter $K=1.1\cdot k$ for both IHT and CG-IRLS. The interpolated plot is shown in Figure 4. It turns out that CG-IRLS has a significantly higher success recovery rate than IHT for less sparse solutions.

3 Algorithm CG-IRLS-λ𝜆\lambda

For the noisy setting we set $\text{MSNR}=10$. According to birits09 ; candes2009 , we choose $\lambda=c\sigma\sqrt{m\log N}$ as a near-optimal regularization parameter, where we empirically determine $c=0.48$. Since we work with relatively large values of $\lambda$ in the regularized problem (49), we cannot use the synthesized sparse solution $x^{*}$ as a reference for the convergence analysis. Instead, we need another reliable method to compute the minimizer of the functional. In the convex case of $\tau=1$, this is performed by the well-known and fast algorithm FISTA beck09 , which shall also serve as a benchmark for the speed analysis. In the non-convex case of $\tau<1$, there is no method which guarantees the computation of the global minimizer, thus, we have to omit a detailed speed-test in this case. However, we describe the behavior of Algorithm CG-IRLS-$\lambda$ for $\tau$ changing.

If the problem is noiseless, i.e., $e=0$, the solution $x^{\lambda}$ of (49) converges to the solution of (1) for $\lambda\rightarrow 0$. Thus, we choose $\lambda=m\cdot 1\text{\sc{e}-}8$, and assume the synthesized sparse solution $x^{*}$ as a good proxy for the minimizer and a reference for the convergence analysis. (Actually, this can also be seen the other way around, i.e., we use the minimizer $x^{\lambda}$ of the regularized functional to compute a good approximation to $x^{*}$.) It turns out that for $\lambda\approx 0$, as we comment below in more detail, FISTA is basically of no use.

As in the previous subsection, we first show that the CG-method within IRLS-$\lambda$ leads to significant improvements in terms of the computational speed. Therefore we choose a noisy trial of Setting B, and compare the computational time of the methods IRLS-$\lambda$, CG-IRLS-$\lambda$, and FISTA. The result is presented on the left plot of Figure 5. We observe, that CG-IRLS-$\lambda$ computes the first iterations in much less time than IRLS-$\lambda$, but due to bad conditioning of the inner CG problems it performs much worse afterwards. Furthermore, as may be expected, the algorithm is not suitable to compute a highly accurate solution. For the computation of a solution with a relative error in the order of $1\text{\sc{e}-}3$, CG-IRLS-$\lambda$ outperforms FISTA. FISTA is able to compute highly accurate solutions, but a solution with a relative error of $1\text{\sc{e}-}3$ should be sufficient in most applications because the goal in general is not to compute the minimizer of the Lagrangian functional but an approximation of the sparse signal.

To further decrease the computational time of CG-IRLS-$\lambda$, we propose the following modifications:

To overcome the bad conditioning in the CG loop, we precondition the matrix $A_{n}=\Phi^{*}\Phi+\operatorname*{diag}\left[\lambda\tau w_{j}^{n}\right]_{j=1}^{N}$ by means of the Jacobi preconditioner, i.e., we pre-multiply the linear system by the inverse of its diagonal, $\left(\operatorname*{diag}A_{n}\right)^{-1}$, which is a very efficient operation in practice.

We introduce the parameter maxiter_cg which defines the maximal number of inner CG iterations and is set to the value $\texttt{maxiter\_cg}=4$ in the following.

The algorithm with modification (i) is called PCG-IRLS-$\lambda$, and the one with modification (i) and (ii) PCGm-IRLS-$\lambda$. We run these algorithms on the same trial of Setting B as in the previous paragraph. The respective result is shown on the right plot of Figure 5. This time, preconditioning effectively yields a strong decrease of computational time, especially in the final iterations where $A_{n}$ is badly conditioned. Furthermore, modification (ii) importantly increases the performance of the proposed algorithm also in the initial iterations. However, again we have to take into consideration that we may violate the assumptions of Theorem 4.1 so that convergence is not guaranteed anymore and failure rates might potentially increase. In the two paragraphs below that are entitled Empirical test on computational time and failure rate with noisy/noiseless data, we present simulations on noisy and noiseless data, which give a more precise picture of the speed and failure rate of the previously introduced methods in comparison to FISTA and IHT.

We investigate the influence of the tolerance $\textnormal{tol}_{n}$ and the number of (inner) iterations of the CG procedure performed along the IRLS iterations in Figure 6 for CG-IRLS-$\lambda$, PCG-IRLS-$\lambda$, and PCGm-IRLS-$\lambda$. We see that the methods do not differ much in terms of the tolerance. In particular CG-IRLS-$\lambda$ and PCG-IRLS-$\lambda$ have nearly the same sequence of $\textnormal{tol}_{n}$, however, due to the bad conditioning, the number of inner CG iterations tremendously increases with growing number of outer IRLS iterations in CG-IRLS-$\lambda$. In contrast, the number of inner CG iterations in PCG-IRLS-$\lambda$ stays very low. A bound on the CG iterations of $\texttt{maxiter\_cg}=4$ does only very slightly change the behavior of $\textnormal{tol}_{n}$ in PCGm-IRLS-$\lambda$ and leads to a further advantage in the computational time, as can be seen in Figure 5.

In the previous paragraph, we observed that the CG-IRLS-$\lambda$ methods are only computing efficiently solutions with a low relative error. Thus we now focus on this setting and compare the three methods PCG-IRLS-$\lambda$, PCGm-IRLS-$\lambda$, and FISTA with respect to their computational time and failure rate in recovering solutions with a relative error of $1\text{\sc{e}-}1$, $1\text{\sc{e}-}2$, and $1\text{\sc{e}-}3$. We only consider the convex case $\tau=1$. Similarly to the procedure in Section 5.2, we run these algorithms on 100 trials for each setting with the respectively chosen values of $\lambda$. In Figure 7 the upper bar plot shows the result for the mean computational time and the lower stacked bar plot shows how often a method was the fastest one. We do not present a plot of the failure rate since none of the methods failed at all. By means of the plots, we demonstrate that both PCG-IRLS-$\lambda$, and PCGm-IRLS-$\lambda$ are faster than FISTA, while PCGm-IRLS-$\lambda$ always performs best.

In the noiseless case, we compare the computational time of FISTA and PCGm-IRLS-$\lambda$ to IHT and IHT+CG-IRLSm. We set $\texttt{maxiter\_cg}=40$ for PCGm-IRLS-$\lambda$. In a first test, we run these algorithms on one trial of Setting A, and C respectively, and plot the results in Figure 8.

As already mentioned, FISTA is not suitable for small values of $\lambda$ on the order of $m\cdot 1\text{\sc{e}-}8$ and converges then extremely slowly, but PCGm-IRLS-$\lambda$ can compete with the remaining methods. IHT+CG-IRLSm is in some settings able to outperform IHT, at least when a high accuracy is needed. PCGm-IRLS-$\lambda$ is always at least as fast as IHT with increasing relative performance gain for increasing dimensions. This observation suggests the conjecture that PCGm-IRLS-$\lambda$ provides the fastest method also in rather high dimensional problems. To validate this hypothesis numerically, we introduce two new high dimensional settings (to reach higher dimensionalities and retaining low computation times for the extensive tests it is again very beneficial to use the real cosine transform as a model for $\Phi$):

We run the most promising algorithms IHT and PCGm-IRLS-$\lambda$ on a trial of the large scale settings D and E. The result, which is plotted in Figure 9, shows that PCGm-IRLS-$\lambda$ is able to outperform IHT in these settings unless one requires an extremely low relative error ($\leqslant 1\text{\sc{e}-}8$), because of the error saturation effect. We confirm this outcome in a test on 100 trials for Setting D and E and present the result in Figure 10.

In the last experiment of this paper, we are interested in the influence of the parameter $\tau$. Of course, changing $\tau$ also means modifying the problem resulting in a different minimizer. Due to non-convexity also spurious local minimizers may appear. Therefore, we do not compare the speed of the method to FISTA. In Figure 11, we show the performance of Algorithm PCGm-IRLS-$\lambda$ for a single trial of Setting C and the parameters $\tau\in\{1,0.9,0.8,0.7\}$ for the noisy and noiseless setting. As reference for the error analysis, we choose the sparse synthetic solution $x^{*}$, which is actually not the minimizer here.

In both the noisy and noiseless setting, using a parameter $\tau<1$ improves the computational time of the algorithm. In the noiseless case, $\tau=0.9$ seems to be a good choice, smaller values do not improve the performance. In contrast, in the noisy setting the computational time decreases with decreasing $\tau$.

Appendix A Proof of Lemma 10

“$\Rightarrow$”(in the case $0<\tau\leqslant 1$) Let $x=x^{{\varepsilon},1}$ or $x\in\mathcal{X}_{{\varepsilon},\tau}(y)$, and $\eta\in\mathcal{N}_{\Phi}$ arbitrary. Consider the function

Now $G_{{\varepsilon},\tau}(0)=0$ and from the minimization property of $f_{{\varepsilon},\tau}(x)$, $G_{{\varepsilon},\tau}(t)\geq 0$. Therefore,

“$\Leftarrow$”(only in the case $\tau=1$) Now let $x\in\mathcal{F}_{\Phi}(y)$ and $\left\langle x,\eta\right\rangle_{\hat{w}(x,{\varepsilon},1)}=0$ for all $\eta\in\mathcal{N}_{\Phi}$. We want to show that $x$ is the minimizer of $f_{{\varepsilon},1}$ in $\mathcal{F}_{\Phi}(y)$. Consider the convex univariate function $g(u)\mathrel{\mathop{:}}=[u^{2}+{\varepsilon}^{2}]^{1/2}$. For any point $u_{0}$ we have from convexity that

because the right-hand-side is the linear function which is tangent to $g$ at $u_{0}$. It follows, that for every point $v\in\mathcal{F}_{\Phi}(y)$ we have

where we have used the orthogonality condition and the fact that $(v-x)\in\mathcal{N}_{\Phi}$. Since $v$ was chosen arbitrarily, $x=x^{{\varepsilon},1}$ as claimed.

References