Iteratively Reweighted Least Squares for Basis Pursuit with Global Linear Convergence Rate
Year: 2022

Abstract:
The recovery of sparse data is at the core of many applications in machine learning and signal processing. While such problems can be tackled using `1-regularization as in the LASSO estimator and in the Basis Pursuit approach, specialized algorithms are typically required to solve the corresponding high-dimensional non-smooth optimization for large instances. Iteratively Reweighted Least Squares (IRLS) is a widely used algorithm for this purpose due to its excellent numerical performance. However, while existing theory is able to guarantee convergence of this algorithm to the minimizer, it does not provide a global convergence rate. In this paper, we prove that a variant of IRLS converges with a global linear rate to a sparse solution, i.e., with a linear error decrease occurring immediately from any initialization, if the measurements fulfill the usual null space property assumption. We support our theory by numerical experiments showing that our linear rate captures the correct dimension dependence. We anticipate that our theoretical findings will lead to new insights for many other use cases of the IRLS algorithm, such as in low-rank matrix recovery.

heading: 1 Introduction
text: The field of sparse recovery deals with the problem of recovering an (approximately) sparse vector x from only few linear measurements, presented by an underdetermined system of linear equations of the form y = Ax. One approach to solve this problem is to consider the `0-minimization under linear constraints, which is NP-hard in general [24, 54]. For computational reasons, instead of `0-minimization, it is common practice to consider its convex relaxation
min x∈RN
||x||1 subject to Ax = y, (P1)
∗kuemmerle@jhu.edu †verdun@ma.tum.de ‡dominik.stoeger@ku.de §The name order of the authors is alphabetical.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
where A ∈ Rm×N , b ∈ Rm are given, which is referred to as `1(-norm) minimization [20, 27, 64] or basis pursuit [18, 19] in the literature.
Unlike `0-minimization, the optimization program (P1) is computationally tractable in general, and a close relationship of their minimizers has been recognized and well-studied in the theory of compressive sensing [12, 13, 26, 32]. In statistics and machine learning, an unconstrained variant of (P1), often called LASSO, amounts to the most well-studied tractable estimator for variable selection in high-dimensional inference [37, 48, 65]. `1-minimization has many other applications and it was even called the modern least squares [14].
The tractability of (P1) becomes evident from the fact that it can be reformulated as a linear program [66]. However, as many problems of interest in applications are high-dimensional and therefore challenging for standard linear programming methods, many specialized solvers for (P1) have been proposed, such as the Homotopy Method [28], primal-dual methods [15, 52], Alternating Direction Method of Multipliers [9], Bregman iterative regularization [71] and Semismooth Newton Augmented Lagragian Methods [46] and Iteratively Reweighted Least Squares (IRLS), the latter of which is in the focus of this paper.
Iteratively Reweighted Least Squares corresponds to a family of algorithms that minimizes nonsmooth objective functions by solving a sequence of quadratic problems, with its idea going back to a method proposed by Weiszfeld for the Fermat-Weber problem [7, 69]. A variety of different problems such as robust regression in statistics [39, 53], total variation regularization in image processing [2, 33, 55], joint learning of neural networks [72], robust subspace recovery [45] and the recovery of low-rank matrices [31, 42, 43, 51] can be solved efficiently by IRLS in practice, as it relies on simple linear algebra to solve the linear systems arising from the quadratic problems at each iteration, without the need of a careful initialization or intricate parameter tuning. On the other hand, the analysis of IRLS methods is typically challenging: General convergence results are often weak, and stronger convergence results are only available in particular cases; see Section 2.3 for more details.
IRLS for sparse recovery. In the sparse recovery context, the first variants of IRLS were introduced in [34, 59] for the `p-quasinorm minimization problem (Pp) with 0 < p ≤ 1 that is similar to (P1), but with ‖x‖p instead of ‖x‖1 as an objective. In [16], modifications of the method of [34, 59] using specific smoothing parameter update rules were observed to exhibit excellent numerical performance for solving (Pp), retrieving the underlying sparse vector when most of the methods fail. A useful fact is that IRLS is one of the few methods (ADMM being the other one [9]) that provides a framework to solve both constrained and unconstrained formulations of `p-minimization problems.
A major step forward in the theoretical understanding of IRLS was achieved in [23], where the authors showed that a variant of IRLS for (P1) converges globally to an `1-minimizer if the measurement operator A fulfills the null space property of sufficient order, which essentially ensures that an `1-minimizer is actually sparse. However, since this proof relies on the existence of a convergent subsequence, their proof does not reveal any rate for global convergence. The analysis of [23] provides, furthermore, a locally linear convergence rate, but this local linear rate has the drawback that it only applies if the support of the true signal has been discovered, which is arguably the difficult part of `0-minimization—cf. Proposition 3.1 below and Section 4.1.
While several extensions and modifications of the IRLS algorithm in [23] have been proposed (see, e.g., [3, 30]), this following fundamental algorithmic question has remained unanswered:
What is the global convergence rate of the IRLS algorithm for `1-minimization?
Our contribution. We resolve this question, formally stated in [62], and present a new IRLS algorithm that converges linearly to a sparse ground truth, starting from any initialization, as stated in Theorem 3.2. Our algorithm returns a feasible solution with δ-accuracy, i.e., ‖x∗ − xk‖1 ≤ δ, where x∗ is the underlying s-sparse vector, in k = O(N √ (logN)/m log(1/δ)) iterations. Analogous to [23], it is assumed that the measurement matrix A satisfies the so-called null space property [21]. We also provide a similar result for approximately sparse vectors. Our proof relies on a novel quantification of the descent of a carefully chosen objective function in the direction of the ground truth. Additionally, we support the theoretical claims by numerical simulations indicating that we capture the correct dimension dependence. We believe that the new analysis techniques in this paper are of independent interest and will pave the way for establishing global convergence rates for other variants of IRLS such as in low-rank matrix recovery [31].
Notation. We denote the cardinality of a set I by |I| and the support of a vector x ∈ RN , i.e., the index set of its nonzero entries, by supp(x) = {j ∈ [N ] : xj 6= 0}. We call a vector s-sparse if at most s of its entries are nonzero. We denote by xI the restriction of x onto the coordinates indexed by I , and use the notation Ic := [N ] \ I to denote the complement of a set I . Furthermore, σs(x)`1 denotes the `1-error of the best s-term approximation of a vector x ∈ RN , i.e., σs(x)`1 = inf{‖x− z‖1 : z ∈ RN is s-sparse}.

heading: 2 IRLS for sparse recovery
text: We now present a simple derivation of the Iteratively Reweighted Least Squares (IRLS) algorithm for `1-minimization which is studied in this paper. IRLS algorithms can be interpreted as a variant of a Majorize-Minimize (MM) algorithm [63], as we will lay out in the following. It mitigates the non-smoothness of the ‖ · ‖1-norm by using the smoothed objective function Jε : RN → R, which is defined, for a given ε > 0, by
Jε(x) := N∑ i=1 jε(xi) with jε(x) := { |x|, if |x| > ε, 1 2 ( x2 ε + ε ) , if |x| ≤ ε. (1)
The function Jε can be considered as a scaled Huber loss function which is widely used in robust regression analysis [40, 50]. Moreover, the function Jε is continuously differentiable and fulfills |x| ≤ jε(x) ≤ |x| + ε for each x ∈ R. Instead of minimizing the function Jε directly, the idea of IRLS is to minimize instead a suitable chosen quadratic function Qε(·, x), which majorizes Jε such that Qε (z, x) ≥ Jε (z) for all z ∈ RN . This function is furthermore chosen such that Qε (x, x) = Jε (x) holds, which implies that min
z∈Rn Qε (z, x) ≤ Jε (x). The latter inequality implies
that by minimizing Qε(·, x), IRLS actually achieves an improvement in the value of Jε as well. More specifically, Qε (·, x) is defined by
Qε(z, x) := Jε(x) + 〈∇Jε(x), z − x〉+ 1
2 〈(z − x),diag(wε(x))(z − x)〉
= Jε(x) + 1
2 〈z,diag(wε(x))z〉 −
1 2 〈x, diag(wε(x))x〉,
(2)
where ∇Jε(x) = ({ xi |xi| , if |xi| > ε xi ε , if |xi| ≤ ε )N i=1 is the gradient of Jε at x and the weight vector wε (x) ∈ RN is a vector of weights such that wε(x)i := [max(|xi|, ε)]−1 for i ∈ [N ]. The following lemma shows that Qε (·, ·) has indeed the above-mentioned properties. We refer to the supplementary material for a proof.
Lemma 2.1. Let ε > 0, let Jε : RN → R be defined as in (1) and Qε : RN ×RN → R as defined in (2). Then, for any z, x ∈ RN , the following affirmations hold:
i. diag(wε(x))x = ∇Jε(x), ii. Qε(x, x) = Jε(x), iii. Qε(z, x) ≥ Jε(z).
As can be seen from the equality in (2), minimizing Qε(·, x) corresponds to a minimization of a (re-)weighted least squares objective 〈·,diag(wε(x))·〉, which lends its name to the method. Note that unlike a classical MM approach, however, IRLS comes with an update step of the smoothing parameter ε at each iteration. We provide an outline of the method in Algorithm 1.
The weighted least squares update (3) can be computed such that xk+1 = W−1k A ∗(AW−1k A ∗)−1(y) with Wk = diag (wk), with the solution of the (m×m) linear system (AW−1k A∗)z = y as a main computational step. This linear system is positive definite and suitable for the use of iterative solvers. In [30], an analysis of how accurately the linear system of a similar IRLS method needs to be solved to ensure overall convergence. We note that for small εk, the Sherman-Woodbury formula [70] can be used so that the calculation of xk+1 boils down to solving a smaller linear system that is well-conditioned, c.f. the supplementary material for details. This numerically advantageous property is not shared by the methods of [3, 23, 30], as our smoothing update (4) is slightly different from the ones proposed in these papers. We refer to Section 2.2 for a discussion.
The update step of the smoothing parameter ε (4) for the IRLS algorithm under consideration requires an a priori estimate of the sparsity of the ground truth of the signal, a piece of information that is
Algorithm 1 Iteratively Reweighted Least Squares for `1-minimization
Input: Measurement matrix A ∈ Rm×N , data vector y ∈ Rm, initial weight vector w0 ∈ RN (default: w0 = (1, 1, . . . , 1)). Set ε0 =∞. for k = 0, 1, 2, . . . do
xk+1 := arg min z∈RN 〈z,diag (wk) z〉 subject to Az = y, (3)
εk+1 := min ( εk,
σs(x k+1)`1 N
) , (4)
(wk+1)i := 1 max ( |xk+1i |, εk+1 ) for each i ∈ [N ], (5) end for return Sequence (xk)k≥1.
also needed by most of the methods for sparse reconstruction. In practice, an overestimation of s is not a problem for similar numerical results if the overestimation remains within small multiples of the sparsity of the signal. We note, however, that there are also versions of IRLS which do not require a-priori knowledge of s, e.g. [30,68], as the update rule for the smoothing parameter is chosen differently. An interesting future research direction is to extend the analysis presented here to IRLS with such a smoothing parameter update.
A consequence of Lemma 2.1, step (4), the fact that ε 7→ Jε(z) is monotonously non-decreasing, and that k 7→ εk is non-increasing is that k 7→ Jε(z) is non-increasing in k. This implies that the iterates xk, xk+1 of Algorithm 1 fulfill
Jεk+1(xk+1) ≤ Jεk(xk+1) ≤ Qεk(xk+1, xk) ≤ Qεk(xk, xk) = Jεk(xk). (6) This shows in particular that the sequence { Jεk ( xk )}∞ k=0
is non-increasing. For this reason, it can be shown that each accumulation point of the sequence of iterates (xk)k≥0 is a (first-order) stationary point of the smoothed `1-objective Jε(·) subject to the measurement constraint imposed by A and y, where ε = limk→∞ εk (see [23, Theorem 5.3]).

heading: 2.1 Null space property
text: As in [23], the analysis we present is based on the assumption that the measurement matrix A satisfies the so-called null space property [21, 36], which is a key concept in the compressed sensing literature (see, e.g., [32, Chapter 4] for an overview). Definition 2.2. A matrix A ∈ Rm×N is said to satisfy the `1-null space property (`1-NSP) of order s ∈ N with constant 0 < ρs < 1 if for any set S ⊂ [N ] of cardinality |S| ≤ s, it holds that ‖vS‖1 ≤ ρs‖vSc‖1, for all v ∈ ker(A).
In [32, Chapter 4], the property of Definition 2.2 was called stable null space property. The importance of the null space property is due to the fact that it gives a necessary and sufficient criterion for the success of basis pursuit for sparse recovery, as the following theorem shows. Theorem 2.3 ([32, Theorem 4.5]). Given a matrix A ∈ Rm×N , every vector x ∈ RN such that ||x||0 ≤ s is the unique solution of (P1) with Ax = y if and only if A satisfies the null space property of order s for some 0 < ρs < 1.
The `1-NSP is implied by the restricted isometry property (see, e.g., [11]), which is fulfilled by a large class of random matrices with high probability. For example, this includes matrices with (sub-)gaussian entries and random partial Fourier matrices [6, 60].

heading: 2.2 Existing theory
text: A predecessor of IRLS for the sparse recovery problem (P1), and more generally, for `p-quasinorm minimization with 0 < p ≤ 1, is the FOCal Underdetermined System Solver (FOCUSS) as proposed by Gorodnitsky, Rao and Kreutz-Delgado [34, 59]. Asymptotic convergence of FOCUSS to a stationary point from any initialization was claimed in [59], but the proof was not entirely accurate, as
pointed out by [17]. One limitation of FOCUSS is that, unlike in IRLS as presented in Algorithm 1, no smoothing parameter ε is used, which leads to ill-conditioned linear systems.
To mitigate this, [16] proposed an IRLS method that uses smoothing parameters ε (such as used in Qε defined above) that are updated iteratively. It was observed that this leads to a better condition number for the linear systems to be solved in each step of IRLS and, furthermore, that this smoothing strategy has the advantage of finding sparser vectors if the weights of IRLS are chosen to minimize a non-convex `p-quasinorm for p < 1.
Further progress for IRLS designed to minimize an `1-norm was achieved in the seminal paper [23]. In [23], it was shown that if the measurement operator fulfills a suitable `1-null space property as in Definition 2.2, an IRLS method with iteratively updated smoothing converges to an s-sparse solution, coinciding with the `1-minimizer, if there exists one that is compatible with the measurements. This method uses not exactly the update rule of (4), but rather updates the smoothing parameter such that εk+1 = min(εk, R(xk+1)s+1/N), where R(xk+1)s+1 is the (s+ 1)st-largest element of the set {|xk+1j |, j ∈ [N ]}. Furthermore, a local linear convergence rate of IRLS was established [23, Theorem 6.1] under same conditions.
However, the analysis of [23] has its limitations: First, there is a gap in the assumption of their convergence results between the sparsity s of a vector to be recovered and the order ŝ of the NSP of the measurement operator. Recently, this gap was circumvented in [3] with an IRLS algorithm that uses a smoothing update rule based on an `1-norm, namely, εk+1 = min(εk, η(1− ρs)σs(x
k+1)`1/N), where η ∈ (0, 1), and ρs is the NSP constant of the order s of the NSP fulfilled by the measurement matrix A—this rule is quite similar to the rule (4) that we use in Algorithm 1. In particular, [3, Theorem III.6] establishes convergence with local linear rate similar to [23] without the gap mentioned above. The main limitation, however, of the theory of [23] (which is shared by [3]) is that the linear convergence rate only holds locally, i.e., in a situation where the support of the sparse vector has already been identified, see also Section 3 and Section 4.1 for a discussion.
We finally mention three relevant papers for the theoretical understanding of IRLS. [5] established the correspondence between the IRLS algorithms and the Expectation-Maximization algorithm for constrained maximum likelihood estimation under a Gaussian scale mixture distribution. By doing so, they established similar results as those from [23], i.e., the global convergence of IRLS with local linear convergence rate (as can be seen from their equation (38), which similar to (7) below) but by using different techniques based on such correspondence. [62] explores the relationship of IRLS for `1-minimization and a slime mold dynamics, interpreting both as an instance of the same meta-algorithm. Without requiring any connection between sparse recovery and `1-minimization, [29] shows that an IRLS-like algorithm for (P1), requires O(N1/3 log(1/δ)/δ2/3 + log(N)/δ2) iterations to obtain a multiplicative error of 1 + δ on the minimizer ||x||1. Unlike our result Theorem 3.2, this corresponds not to a linear, but to a sublinear convergence rate.

heading: 2.3 Related work
text: As mentioned in the introduction, IRLS has a long history and has appeared under different names within different communities, e.g., similar algorithms are usually called half-quadratic algorithms in image processing [2, 41] and the Kac̆anov method in numerical PDEs [25]. Probably the most common usage of IRLS has been in robust regression [35, 39], c.f. [10] for a survey that also covers applications in approximation theory. For p-norm regression, [1] proposed a version of IRLS for which convergence results for p ∈ [2,∞) were established, solving a problem that was open for over thirty years. Also, for robust regression, by using an `1-objective on the residual, [53] showed recently global convergence of IRLS with a linear rate, with high probability for sub-Gaussian data. We note that our proof strategy is different from the one of [53] due to a structural difference of (P1) from robust regression.
In [57], the authors provide a general framework for formulating IRLS algorithms for the optimization of a quite general class of non-convex and non-smooth functions, however, without updated smoothing. They use techniques developed in [4] to show convergence of the sequence of iterates to a critical point under the Kurdyka-Łojasiewicz property [8]. However, no results about convergence rates were presented.
For the sparse recovery problem, the topic discussed in this paper, the references [30, 44, 68] analyzed IRLS for an unconstrained version of (P1), which is usually a preferable formulation if
the measurements are corrupted by noise. Additionally, the work [30] addressed the question of how to solve the successive quadratic optimization problems. The authors developed a theory that shows, under the NSP, how accurately the quadratic subproblems need to be solved via the conjugate gradient method in order to preserve the convergence results established in [23].
Finally, for the related problems of low-rank matrix recovery and completion, IRLS strategies have emerged as one of the most successful methods in terms of data-efficiency and scalability [31, 42, 43, 51].
While we were writing this paper, the manuscript [58] appeared providing new insights about IRLS. It describes a surprisingly simple reparametrization of the IRLS formulation for `p-minimization (with p ∈ (2/3, 1)) that leads to a smooth bilevel optimization problem without any spurious minima, i.e., the stationary points of this new formulation are either global minima or strict saddles. It is an interesting future direction to explore the connection between this new approach and our global convergence theory.

heading: 3 IRLS for Basis Pursuit with Global Linear Rate
text: As discussed in Section 2.2, the main theoretical advancements for IRLS for the sparse recovery problem were achieved in the work [23]. Proposition 3.1. [23, Theorem 6.1] Assume that A ∈ Rm×N satisfies the NSP of order ŝ > s with constant ρŝ such that 0 < ρŝ < 1− 2ŝ+2 and ŝ > s+ 2ρŝ 1−ρŝ hold. Let x∗ ∈ R
N be an s-sparse vector and set y = Ax∗. Assume that there exists an integer k0 ≥ 1 and a positive number ξ > 0 such that
ξ := ‖xk0 − x∗‖1
mini∈S |(x∗)i| < 1. (7)
Then the iterates {xk0 , xk0+1, xk0+2, . . .} of the IRLS method in [23] converge linearly to x∗, i.e., for all k ≥ k0, the kth iteration of IRLS satisfies
‖xk+1 − x∗‖1 ≤ ρŝ(1 + ρŝ)
1− ξ
( 1 +
1
ŝ− 1− s
) ‖xk − x∗‖1. (8)
The main contribution of this paper is that we overcome a local assumption such as (7) and show that IRLS as defined by Algorithm 1 exhibits a global linear convergence rate, i.e., there is a linear convergence rate starting from any initialization, as early as in the first iteration.
Exactly sparse case. Our first main result, Theorem 3.2, deals with the scenario that the ground truth vector x∗ is exactly s-sparse. Our second result, presented in the supplementary material, generalizes the first one to the more realistic situation of approximately sparse vectors.
Theorem 3.2. Consider the problem of recovering an unknown s-sparse vector x∗ ∈ RN from known measurements of the form y = Ax∗. Assume that the measurement matrix A ∈ Rm×N fulfills the `1-NSP of order s with constant ρs < 1/2. Let the IRLS iterates { xk } k
and {εk}k be defined by the IRLS algorithm (3) and (4) with initialization x0. Then, for all k ∈ N, it holds that
Jεk(xk)− ‖x∗‖1 ≤ (
1− c ρ1N
)k ( Jε0(x0)− ‖x∗‖1 ) (9)
as well as
‖xk − x∗‖1 ≤ 9 (
1− c ρ1N
)k ‖x0 − x∗‖1. (10)
Here c = 1/768 is an absolute constant and ρ1 < 1/2 denotes the `1-NSP constant of order 1.
Inequality (9) says that the difference Jεk(xk)−‖x∗‖1 converges linearly with a uniform upper bound of 1− cρ1N on the linear convergence factor. As our proof, which is detailed in the supplementary material, shows, this implies inequality (10), which implies that also ‖x∗ − xk‖1 exhibits linear convergence in the number of iterations k. In particular, this means that for some error tolerance δ > 0, we obtain ‖x∗ − xk‖1 ≤ δ after O ( ρ1N log ( ‖x∗−x0‖1
δ
)) iterations.
Remark 3.3. Note that it follows directly from Definition 2.2 that the constant ρ1 of the `1-NSP of order 1 satisfies ρ1 ≤ ρs ≤ 1, which implies that δ-accuracy is obtained after O ( N log ( ‖x∗−x0‖1
δ )) iterations. This bound can be improved in many scenarios where one can obtain more explicit bounds on ρ1, for example, when A is a Gaussian matrix. Namely, inspecting [32, p. 142 and Thm. 9.2], we observe in this scenario that ρ1 . √ (logN)/m with high probability. Hence, in this scenario, at
most O ( N √
logN m log
( ‖x∗−x0‖1
δ
)) iterations are needed to achieve δ-accuracy.
The key idea in our proof is to use fact that the quadratic functional Qεk(·, xk) approximates the `1-norm in a neighborhood of the current iterate xk. For this reason, we also expect that for t > 0 sufficiently small, we have that Qεk(x k + tvk, xk) < Qεk(x k, xk) if vk = x∗ − xk is the vector between xk and the ground truth x∗. Then by choosing t properly, we can guarantee a sufficient decrease of the functional Jεk ( xk ) in each iteration.
In Section 4, we conduct experiments that indeed verify the linear convergence in (9) and (10). Moreover, we study numerically whether one can observe a dependence of the convergence rate on the problem parameters N , s and m. We construct a worst-case example which indicates that the convergence rate indeed may depend on the dimension N in a way as described by (9). In a certain sense, this indicates that there are two convergence phases, a global one, where a dimensiondependent constant cannot be avoided and a local convergence phase, where a local convergence result such as described in Proposition 3.1 kicks in.
Finally, let us mention that we have undertaken no efforts to optimize the constant c = 1/768 in Theorem 3.2. Nevertheless, we note that the constant c can be replaced by the sharper constant cρs as defined in Proposition B.3.

heading: 4 Numerical experiments
text: In this section, we support our theory with numerical experiments. First, we examine whether IRLS indeed exhibits two distinct convergence phases, a “global” one, as described in this paper, and a local one, as described in [3, 23], corresponding to different linear convergence rate factors. Second, we explore to which extent the dimension dependence in the convergence rates (9) and (10) indicated by Theorem 3.2 is necessary, or if we rather can expect a dimension-free linear convergence rate factor. All experiments are conducted on an iMac computer with a 4 GHz Quad-Core Intel Core i7 CPU, using MATLAB R2020b.

heading: 4.1 Local and global convergence phase
text: We first note that the local convergence result of [23, Theorem 6.1] depends on the locality condition ξ(k) := ‖x
k−x∗‖1 mini∈S |(x∗)i| < 1, cf. (7). Under this condition (and an appropriate null space condition),
Proposition 3.1 stated above implies that ‖xk+1 − x∗‖1 ≤ µ‖xk − x∗‖1 with an absolute constant µ < 1 which, in particular, does not depend on the dimension N,m, and s. This corresponds to a locally linear rate for IRLS. A very similar condition to (7) is required by the comparable and more recent local convergence statement [3, Theorem III.6, inequality (III.14)] for the IRLS variant considered in [3].
However, a closer look at the locality condition (7) reveal that its basin of attraction is very restrictive: This condition means that the support identification problem underlying the sparse recovery has already been solved, as can be seen from the following proposition, whose proof we provide in the supplementary material.
Proposition 4.1. Let xk, x∗ ∈ RN , let S ⊂ [N ] be the support set of x∗ of size |S| = s. If (7) holds, i.e., if ‖xk − x∗‖1 < mini∈S |(x∗)i|, then the set Sk ⊂ [N ] of the s largest coordinates of xk coincides with S.
We now explore the behavior of the IRLS algorithm for `1-minimization, Algorithm 1, and the sharpness of Proposition 4.1 in experiments that build on those of [23, Section 8.1]. For this purprose, forN = 8000, we sample independently a 200-sparse vector x∗ ∈ RN with random support S ⊂ [N ], s = 200 = |S|, chosen uniformly at random such that (x∗)S is chosen according the Haar measure
on the sphere of a 200-dimensional unit `2-ball, and a measurement matrix A ∈ Rm×N with i.i.d. Gaussian entries such that Aij ∼ N (0, 1/m), while setting m = b2s log(N/s)c. Such a matrix is known to fulfill with high probability the `1-null space property of order s with constant ρs < 1 [32, Theorem 9.29].
In Figure 1, we track the decay of the `1-error ‖xk − x∗‖1 of the iterates xk returned by Algorithm 1 via the values of ζ(k) := ‖xk − x∗‖1/mini∈S |(x∗)i|, depicted in red, and the behavior of the factor µ`1(k) := ‖xk − x∗‖1/‖xk−1 − x∗‖1, depicted in blue. We observe that the condition (7) for local convergence with the fast, dimension-less linear rate (8) is satisfied after k = 33 iterations, as indicated by the vertical dashed red line.
In the first few iterations, ζ(k) is larger than 1 by several orders of magnitudes, suggesting that the local convergence rate results of [3, 23] do not apply until the later stages of the simulation: In fact, we observe that the support S of x∗ is already perfectly identified via the s largest coordinates of xk as soon as k ≥ 18. For iterations 18 ≤ k ≤ 50, the linear rate µ`1(k) remains very stably around ≈ 0.7, after which an accelerated linear rate can be observed.5 Before k = 18, for this example, the rate µ(k) hovers around 0.7 with slight variations. For all iterations k, µ(k) is smaller than 1, in line with the global linear convergence rate implied by Theorem 3.2.
Repeating a similar experiment for a larger ambient space dimension N = 16000 and a smaller measurement-to-sparsity ratio such that m = b1.75s log(N/s)c results in a qualitatively similar situation, as seen in Figure 2(a): In Figure 2(a), we add also a plot of the linear convergence factor µ(k) := Jεk (x
k)−‖x∗‖1 Jεk−1 (xk−1)−‖x∗‖1 that tracks the behavior of the linear convergences in the smoothed
`1-norm objective J , cf. (17). In addition to what have been observed in Figure 1, we see that µ(k) and µ`1(k) exhibit a very similar behavior for this example.
Hence, these experiments indicate that we can distinguish two phases. In the first, global phase linear convergence already sets in, but the instantaneous linear convergence rate has not yet stabilized. In the second one, when the support identification problem has been solved, the instantaneous linear convergence stabilizes.
Remark 4.2. There are other methods in the literature, such as proximal algorithms, for which convergence results with a two phase behaviour were already established. For example, [47] showed that a forward-backward method applied to the Lasso problem exhibits local linear convergence, and that after a finite number of iterations, the region of fast convergence is reached. In particular, [47, Proposition 3.6(ii)] provides a bound on this number of iterations, which scales proportionally with ||x∗−x0||22. On the other hand, (3.2) for IRLS provides a bound on the number of iterations until the fast linear convergence rate is reached that scales proportionally with log(||x∗ − x0||2), but also proportionally with the dimension N . Moreover, most of these results require stronger assumptions than the NSP, such as the restricted isometry property or a restricted strong convexity/smoothness property.
5The latter phenomenon cannot be observed for the IRLS algorithm of [23] as it uses a slightly different objective function than Algorithm 1.
Linear convergence factors µ`1(k) := ‖xk−x∗‖1 ‖xk−1−x∗‖1 (in blue) and µ(k) := Jεk (x k)−‖x∗‖1 Jεk−1 (xk−1)−‖x∗‖1 (in green), filled circles if Sk = S (perfect support identification), and error parameter ζ(k) := ‖xk − x∗‖1/mini∈S |(x∗)i| (in red), horizontal and vertical red lines as in Figure 1.

heading: 4.2 Global convergence rate and its dimension dependence
text: In this section, we explore to which extent the dependence on N in the convergence rates (10) and (14) is necessary or if we can rather expect a dimension-free linear convergence rate factor. To this end, we run a variation of IRLS that initializes the weight vector w0 ∈ RN not uniformly as in Algorithm 1, but based on an adversary initialization, here denoted by zadv. More specifically, we first compute a minimizer
zadv ∈ arg min z∈RSc :AScz=y ‖z‖1
of the `1-minimization problem restricted to the off-support coordinates of x∗ indexed by Sc and set then x0 ∈ RN such that x0Sc := zadv and x0S = 0. Based on this initialization x0, we compute ε0 := σs(x 0)`1 N and set the first weight vector such that for all i ∈ [N ],
(w0)i := 1
max (|x0i |, ε0) , (11)
before proceeding with the IRLS steps (3), (4) and (5) until convergence.
We observe in Figure 2(b) that this initialization, which is adversary as it sets very large initial weights on the coordinates of S that correspond to the true support of x∗, eventually results in the same behavior of Algorithm 1 as for the standard initialization by uniform weights, identifying the true support at iteration k = 39 compared to k = 30. However, in the first few iterations, we see that the instantaneous linear convergence factor µ(k) is close to 1 with µ(1) = 0.980, decreasing only slowly before stabilizing around 0.79 after around k = 30.
While this is just one example, this already indicates that in general, a linear rate such as (8), i.e., without dependence on the dimension N (which has been proven locally in [23, Theorem 6.1] and [3, Theorem III.6]) might not hold in general.
In our next experiment, we further investigate numerically the dimension dependence of the worstcase linear convergence factor µ(k) := Jεk (x k)−‖x∗‖1
Jεk−1 (xk−1)−‖x∗‖1 , which is upper bounded by the result of
Theorem 3.2. We saw that in the experiment using the adversary initialization mentioned above and depicted in Figure 2(b), the maximal value was attained in the first iteration, i.e., for µ(1), as the effect of the adversary initialization is most eminent for k = 1.
We now run IRLS starting from the adversary initialization for different ambient dimensions N = 125 · 2`/2 for ` = 0, 1, . . . , 14. For each of the values of N , we sample vectors x∗ ∈ RN of sparsity s = 40 from the same random model as above, and scale the number of i.i.d. Gaussian measurements with m = b2s log(N/s)c. We average the resulting values for µ(1) across 500 independent realizations of the experiment.
Jε0 (x0)−‖x∗‖1 is the linear convergence
factor, for IRLS initialized from adversary initialization.
In Figure 3, we see that dependence on N of linear convergence factor µ(1) that is observed for this experiment is quite well described by the upper bound (9) provided by our main result Theorem 3.2, as 11−µ(1) scales almost linearly with N . As a footnote in Section 3 indicates, the constant ρ1 of the
null space property of order 1 scales with √
logN m , and therefore a precise dependence on all the
parameters such as m and s might be more complicated than what can be observed in this experiment.
Nevertheless, we interpret Figure 3 as strong evidence that the linear convergence rate factor of Proposition B.3 is tight in its dependence on N , and that a dimension-less factor µ cannot be expected in general.
Conclusion
In this paper, we solved an open problem in the algorithmic theory for sparse recovery. In particular, we established a new variant of the IRLS algorithm for Basis Pursuit or `1-minimization for which we show a global linear convergence rate under a suitable and sharp assumption, namely, the null-space property. Moreover, we have corroborated our theory with numerical experiments that, first, discussed the difference between the local and global convergence phase and, second, that elucidated the optimality of the dimension dependence of convergence rate given by our main theorem.
We think that the results in this paper give rise to a number of interesting research directions for follow-up work. While the numerical experiments in Section 4 substantiate the hypothesis that the dependence of the convergence rate on N and ρ1 in our theory is not an artifact of our proof, we also observed in this section that for a generic initialization no such dependence can be observed. In view of this, it is interesting to investigate whether a dimension-independent global convergence rate is possible, for example via a smoothed analysis [22, 61]. Furthermore, there are currently no convergence rates available for IRLS optimizing a nuclear norm-type objective, which is of great interest for low-rank matrix recovery [31, 43, 51], and we expect that our analysis may be generalizable to this setting as well.
Finally, it was observed that sparse vectors can be recovered from even fewer measurements via the optimization of a non-convex `p-quasinorm (with 0 < p < 1), and that IRLS exhibits excellent performance in this case [16, 23]. While a thorough understanding has remained elusive so far for this non-convex case, we consider our results as a first step towards a global convergence theory for IRLS for the optimization of `p-quasinorms or similar non-convex surrogate objectives.
Acknowledgments and Disclosure of Funding
The authors want to thank Massimo Fornasier for inspiring discussions and the anonymous reviewers for their detailed comments. C.K. is grateful for support by the NSF under the grant NSF-IIS-1837991. C.M.V. is supported by the German Federal Ministry of Education and Research (BMBF) in the context of the Collaborative Research Project SparseMRI3D+ and by the German Science Foundation (DFG) in the context of the project KR 4512/1-1. D.S. was supported by the Air Force Office of Scientific Research under award #FA9550-18-1-00788.

Review summaries:
In this paper, the authors study the convergence rate of the Iterative Reweighted Least Square to solve Basis Pursuit. After recalling the existing results and describing the IRLS algorithm with a specific smoothing parameter updates, the authors state their main result showing the global convergence of IRLS under the $\ell_1$-null space property. Finally, some  numerical experiments are presented to corroborate the theoretical results, showing a two stage convergence, one global described by this paper and one local, once the support is identified, which has been studied previously.

The paper studies the Basis Pursuit (BP) optimization problem (minimal l1 norm solution to a linear system). Iterative Reweighted Least Squares (IRLS) is a majorization minimization algorithm to solve this problem. The paper shows that it converges *globally* and linearly, under a nullspace property assumption (with parameter <1/2). Existing results only prove local linear convergence rates, first showing that IRLS identifies the support of the solution, after which the rate becomes linear. The key is a new update rule for the smoothing parameter used in the majorization of |.| by a scaled Huber function. Experiments highlight the two regimes : global convergence shown in this paper, followed by a linear convergence with a better rate once the iterate enter a basin of attraction of the solution.  

The authors propose an iteratively reweighted least squares (IRLS) algorithm for the basis pursuit (BP) problem. The proposed algorithm differs from existing IRLS variants for basis pursuit mainly in terms of the iterative update of the smoothing parameter $\varepsilon_k$. Provided that the measurement matrix $A$ satisfies the null space property (NSP) of order $s$ with constant $\rho_s = \frac12$, the authors prove that their method converges globally at a linear rate and can recover any $s$-sparse vector $x^*$ from an arbitrary initial point $x_0$. The main contribution of the paper is the establishment of a global linear convergence rate for an IRLS algorithm for BP.