Global Linear and Local Superlinear Convergence of IRLS for Non-Smooth Robust Regression
Authors: Liangzu Peng
Christian Kümmerle
Rene Vidal
Venue: NeurIPS 2022 Accept (2022)
Keywords: Convergence Rate Analysis
Non-Smooth Optimization
Robust Regression
Outliers
Iteratively Reweighted Least-Squares
Sparsity

Abstract: We advance both the theory and practice of robust $\ell_p$-quasinorm regression for $p \in (0,1]$ by using novel variants of iteratively reweighted least-squares (IRLS) to solve the underlying non-smooth problem. In the convex case, $p=1$, we prove that this IRLS variant converges globally at a linear rate under a mild, deterministic condition on the feature matrix called the stable range space property. In the non-convex case, $p\in(0,1)$, we prove that under a similar condition, IRLS converges locally to the global minimizer at a superlinear rate of order $2-p$; the rate becomes quadratic as $p\to 0$. We showcase the proposed methods in three applications: real phase retrieval, regression without correspondences, and robust face restoration. The results show that (1) IRLS can handle a larger number of outliers than other methods, (2) it is faster than competing methods at the same level of accuracy, (3) it restores a sparsely corrupted face image with satisfactory visual quality.

Global Linear and Local Superlinear Convergence of IRLS for Non-Smooth Robust Regression
Liangzu Peng
Christian Kümmerle
Mathematical Institute for Data Science Johns Hopkins University 1peng258 jhhu. edu
Department of Computer Science University of North Carolina at Charlotte
kuemme 1 eluncc. edu
René Vidal
Mathematical Institute for Data Science Johns Hopkins University rvidalg edu
Abstract
We advance both the theory and practice robust f, Ep-quasinorm regression for p (0, 1] by using novel variants of iteratively reweighted least-squares (IRLS) solve the underlying non-smooth problem. In the convex case, p = 1, we prove that this IRLS variant converges globally linear rate under mild, deterministic condition on the feature matrix called the stable range space property. the non- convex case, p E (0, 1). prove that under similar condition, IRLS converges locally to the global minimizer at a superlinear rate of order 2 - p: the rate becomes quadratic as p - O. We showcase the proposed methods in three applications: real phase retrieval, regression without correspondences, and robust face restoration. The results show that (1) IRLS can handle larger number of outliers than other methods, (2) it is faster than competing methods at the same level of accuracy, (3) it restores a sparsely corrupted face image with satisfactory visual quality. https: / http//gttuuuttt thub. com 1 www.gnnnngaaaaae.aae.aaaeiauauaa TRLS-MA NeSS-NNurrPP2000
Introduction
Given feature matrix A E Rexn with m >> and response vector y R.m. the problem
min #ER"
Aa -yll yll,
of fp-regression has a variety of different applications depending on the choice of ||0/| =
(E, (u,P)) 3/P, While p = 2 corresponds to standard linear regression, the choice of p > 2 arises naturally in semi-supervised learning on graphs [1, 2, 3], and a lot of activity has been dedicated recently to the computational complexity analysis for the case p > [4, 6,
In this paper, we assume there is a coefficient vector x E IR" such that the residual r= Ax k-spaase, in which case choice p E (0, of interest. Indeed, for the convex and non non-smooth case p = 1, (1) is known as least absolute deviation, which dates back to the time of Boscovich around in the middle of the 18th century [7, 8). Since then, it has been known intuitively that (1)
(p = 1) is robust to large but few measurement errors (quoting [9]), i.e., that (1) is robust to outliers.
The algorithmic and theoretical understanding of (1) for p 1 has been of long-standing interest to statisticians, which has led to vast literature spanning from the classic 1964 paper of Huber [10] recent contributions [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. That being said, for p there is
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

Algorithm 1: IRLS for Ey Regressies (IRLSp)
Input: A [aI...., aml E R y= WV, Umll R." E 10;
2 Weight initialization (0)
t [w] (0) ;...; wm ] E R // Instead, one can initialize some vector x: (1)
For & + 0, l....
x: (1+1) argmin
Σ ii (a. x - yi)
(2)
1=1
Update ((++1) suitably based on (1+1), A, and y: // See Section 2.2 for details.
Wi (++1) t- max {lag ((+1) vel. ((+1)) p-2 Vi = 1, , m
(3)
close, but somewhat underexplored connection between (1) and the basis pursuit problem [21], and more specifically, the theoretical [22. 23, 24] and algorithmic aspects [25, 26, 27] of compressed sensing [28] (see also Section 2.1). For p E (0. 1). problem (1) is non-smooth and non-convex, and generally less well understood. To the best of our knowledge, the only paper that considered (1) with p E (0, 1) is [29], where the authors presented a condition (based on restricted isometry constants) that guarantees exact recovery of x from (1). Other related works come either from the compressed sensing literature, where an E (and noisy) version of basis pursuit has been considered as aparsity-promottin formulation [30, 31, 32, 26, 33, 34, 35, 36, 37J, from the literature of matrix recovery/completio, where the Schatten-p norm has come into play as a non-convex surrogate for rank minimization [38, 39, 40, 41, 42, 43]. A key message from these works is that lp Schatten-p coomation-thhorette properties (e.g., requires fewer
minimization with p E (0, 1) offers better samples for exact recovery) than minimization with 1.
Here, we study an iteratively reweighted least-squares method (IRLS,, to solve (1) with p E (0. 1]. As listed in Algorithm 1, I alternates between solving weighted least-squares problem (2) and updating the weights (3); see Section 2.2 for more elaboration on IRLS. The simplicity of this idea (with its impressive performance) justifies its popularity in many machine learning [44, 45, 46, 47] and computer vision [48, 49, 50, 51, 52] applications. Based on recent advances on IRLS [26. 27], we make several contributions for understanding (1) and IRLSp. We state our contributions next.
The Stable Range Space Property (Section 3.1). We put forward the use of the (stable) range space property (RSP) for studying the robust regression problem (1). The stable RSP was proposed by [53] in different context to analyze compressed sensing algorithm for solving weighted basis pursuit problem at each iteration. In analogy to the nullspace property in compressed sensing [24, 28], we show that the RSP is necessary and sufficient condition for guaranteeing that the fp-regression problem in (1) admits x² as its unique solution (Proposition 2). Moreover, we show that if A E mm has i.idd (O, 1) entries and if m is large enough, then the (stable) RSP holds with high probability (Proposition 3). This justifies its use as the core assumption in our analysis.
Global Linear Convergence (Section 3.2). We prove in Theorem that, under a stable RSP assumption and with e(t) suitably updated, the IRLS, Algorithm with p = 1 converges linearly to the ground-trum truth x² from any initial weight (0))) (or equivalently from any initial point 2e()). Note that while (accelerated) first-order methods (e.g.. (sub-)gradient descent and proximal algorithms) can also solve the convex and non-smooth problem (1) with p = 1, they exhibit, in general, (global) sub-linear rates at best [54, 55]. To our knowledge, [17] is the only paper that claims global linear convergence of different IRLS variant for (1) with p = 1; we compare our results with the ones of [17] in Section 3.2. On the other hand, Theorem is inspired by the IRLS method of [27] for basis pursuit; based on [27], we suitably modify their proof strategy, and thus obtain a faster linear rate.
Local Superlinear Convergence (Section 3.3). We prove Theorem 2 that, under stable RSP assumption and with e(t) suitably updated, the IRLS, Algorithm 1 with p E (0, 1) converges to x superlinearly, provided that (1) falls into a certain neighborhood x=" To the best of our knowledge, no similar result exists for e--reressioon (1). While Theorem 2 is inspired by the IRLS algorithm of [26] for ba-basis pursuit, their IRIS method does not work well for small p (Section 2.2), and unlike
in our work, their radius of local convergence diminishes greatly p + 0 or + co.
2

---

Applications (Section 4). We illustrate the performance of IRLSp for real phase retrieval [56, 57], linear regression without correspondences [58. 59, 60, 61], and face restoration from sparsely corrupted measurements [62, 63, 64]. For real phase retrieval (Section 4. we show that 1 needs only m = 2n measurements to recover x: up to sign, with an additional assumption. Theoretically, even brute-force can fail to identify -xx for fewer than 2n- measurements [65]. Empirically, many methods, including Kaczmarz [66, 67, 57], PhaseLamp [68], truncated Wirtinger flow [69], and coordinate descent [70], fail with For linear regression without correspondences (Section 4.2), we show in Figures 2b-2c that, I is uniformly faster (20-100x) and more accurate than PDLP [71] (merged into Google or-tools) and the commercial solver Gurobi [72], both of which solve (1) with p = 1 as a linear program, and also than subgradient descent implemented by Beck & Guttmann-Bes [73]. For face restoration (Section 4.3), we present both quantitative and qualitative results on the Extended Yale dataset [74].
2nz Gaussian measurements (Figure 2a).
2 Background
2.1 Connection of Robust Regression and Compressed Sensing
The ep-regression problem (1) has a natural correspondence to the sparse recovery problem
min MER
|r|l,
s.t. Dr =
(4)
where D is (m- -n) n) x m matrix, 2 E IR=-" vector and the objective ||lle penalizes coefficient vectors with many non-zero coordinates. For
po = 1, (4) is also called basis pursuit [75]. Problems (1) and (4) are known to be related in the following sense (as implicitly stated in [21, 291):
Proposition 1. Suppose that the runge space of equal to the nullspace of D, and that z = Dy. If œ globally minimizes (I), then rı == Aœi globully minimizes (4), On the other hand, if r2 globally minimizes (4), then there exists some 22 with r22 = Ar y that globally minimizes (1).
Proposition sheds light on how we "transfer", with new insights, from the analysis of [26, 27] for (4) to results for (1). In what follows, we treat [26, 27] in the context of robust regression and highlight the contribution of our work relative to [26, 27] and other existing results, whenever possible.
2.2 Iteratively Reweighted Least-Squares: The Basics and New Insights
We first discuss two different updating rules for the smoothing parameter ((t) of Algorithm 1: the fixed rule and the dynamic rule. We then consider the weight updating strategy (3), Along the way we use synthetic experiments to illustrate ideas; see Appendix ?? for the experimental setup. In doing so, we intend to provide review of the state-of-the- on the variants of IRLS.
Fixed Smoothing Parameter. Instead of updating c(9) at each iteration, most works on IRLS use a fixed and small positive number, e.g., (0)) 0.001 for each [44, 76, 45, 47, 49, 51, 50, 77). The intuition is avoid division by the potentially very small residual |a, (t) ya-. But this common practice of fixing comes with at least three issues. First, IRLS with a fixed c > 0 converges only to "c-approximate" point, not exactly the global minimizer ([78, 79, 44, 76). Figure la). Second, even if setting e small (e.g.. = 10-15) could lead to an accurate enough solution for 1 (Figure la), it can fail for p (Figure lb). Finally, it makes obtaining a global linear convergence rate guarantee difficult: We are not aware of any theory about global linear rate of IRLS with fixed e.
Dynamic Smoothing Parameter. Researchers have used different insights to reach the consensus of dynamically updating e(t) [26, 17, 80. 27, 81]. The first insight is that convergence to global minimizers ensues if (() is suitably decreased to 0 [26). With / E (0, 1). one such decreasing rule is
((+1) + (9) if certain conditions are satisfied, or keep (1+1) +- (1) otherwise [17, 81].
(5)
Figure Ia depicts the performance of [17] with an arbitrary choice B = 0.5 and (0) = 1. Also shown in Figure la are the IRLS methods with update rules of [26] and [27]. which we discuss next.
The basis pursuit paper [26] also proposed dynamic updating rule for e In our context (1), it sets (0) t- 00, (++1) + min ((()) [m((1) a+1/m]]
(t+1) = Az (++1) - y E R." (6)
3

---

10-2-
10
Relative Error
10
IRLS (fa... 10-15 16y
10
27] <6.
117]
IRLS 16..1.
(o 7)
[26] (6)
[6]] Ma... 6)
[27] (fo.s: 7)
[27] (7)
[26] (fa. ()
IRLS (r-1E-5)
[27] (aa.-
A IRLS (c=1E-15)
IRLS (c=1E-10) IE-10)
[27] (fo 7)
10-16
16
16
10
20
30
40
50
200
400
600
100
300
500
700
900
Number of Iterations
(a) p
(b) < 1
(c) Sensitivity to a (k 200)
Figure 1: Figure la: The relative error ||x(() x² / //|l="
2l2 at cach iteration t of different IRLS variants for & regrreessio (k 200). Figure 1b: The relative error of IRLS, Algorithm I and that of [26] for re-resression (50 iterations). Figure Ic: The sensitivity of (6) and (7) to miissppeciiiattinn of (50 iterations). In Figure 1, we
1000, 10, and results are averaged over 20 trials.
where the hyper-ypaaee parameter' a is a non-negative integer, and [r(+++]]++] is the (a + 1)-th largest element of the residual (t-1)) in absolute values. Clearly, (6) creates a non-increasing sequence of smoothing parameters e()), and the decay rate of (1) is adaptive to the data (in this case to the residual ())+))) While at first glance it is not clear how fast e**) decays, [26] showed that the decay rate of (9) in (6) is locally linear for p = 1, a = under mild conditions, in accordance with the decay of the objective values.
The final dynamic update rule for the parameter (() is from [27]' and improves upon (6) [26] via:
( 20, c(++1) min ( ((+1))//1))) /m}
where o (1+1) min {[]nltt)) =||,
T is a-sparse
Observe that (7) computes the -n-mr (l+1) of the best a-term approximation of the residual p(++), which is in general larger than [r(++1)] of (6), while both update rules (6) and (7) rely on the thyper-paramed
a. Ideally, if k, and if IRLS with (6) or (7) converges to x=", then both [r(++1)]a++ 1++1 and (t++)) will approach 0 for large enough f. While update rules (6) and (7) perform similarly for é f-rreressioon (Figure la), (6) fails more easily for small p than (7) (Figure lb). Finally, note that overestimating the sparsity level & by setting a > k deteriorates the performance of (7) only slightly (Figure Ic).
Summary. In Section 2.2 delivered two take-away messages: (1) IRLS with fixed (9) finds some approximate solution, sometimes good enough, (2) dynamically updating e( per (5)-(7) leads to local [26] or global [17, 27] linear convergence guarantees. In view of Figure 1, we will next consider the IRLSp Algorithm 1 with update rule (7) for e and analyze its convergence rates.
3 Convergence Theory of IRLS for Robust Regression
In Section 3.1 we introduce the stable range space property and justify it as our core assumption. Under this assumption, we prove the global linear convergence (Theorem 1) and local superlinear convergence (Theorem 2) of the IRLS, Algorithm in Sections 3.2 and 3.3, respectively.
3.1 The Stable Range Space Property
Definition (Stable Range Space Property). A matrix A E m is said to satisfy the range space property (RSP) of order or the k-RSP for short, if the following holds for any vector d in the range
Inn [26] and [27]. is the sparsity level k but l might be unknown in practice so we treat it as a hyper- parameter. That being said, in the experiments we set = le by default, for simplicity and in light of Figure Ic.

---

space of A and any set S C (I....,m, m of cardinality at most k:
Eidl Σ [d:|
(8)
iES
The stable RSP of A is defined the same way as the RSP except that we now require
Eld S di
iES"
for some E (). We write "(kt (7 stable RSP" to emphasize the parameters i and
By definition, the stable RSP implies the RSP. Note that the (k, -ttbbbe RSP was proposed in [53] to analyze a reweighted ( f--minimization algorithm for compressed sensing. With the notation of Proposition 1, we for all d in the nullspace of D, the latter being the well-known nullspace property (NSP) [24, 28). In other words, the (stable) RSP of A is equivalent to the (stable) NSP of D.
that checking (8) for all d in the range space of A is equivalent to checking
In Sections 3.2 and 3.3, we will use the (kt -tabbelle RSP as an assumption for analysis. Arguably, this is a weak assumption to make the (stable) RSP is very close to sufficient and necessary condition for exact recovery of the coefficient vector via fp-robust regression (I):
Proposition (Exact Recovery RSP). Let p (0, 1). For all x² and y such that Ax" y is k-sparse x² is the unique solution to (1) if and only if A E R=xx" satisfies the k-RSP.
follows from [82, Theorem 5] and Proposition that checking whether a given matrix A satisfies the stable RSP is co-NP-complete [83]. On the other hand, show that random Gaussian matrices of size m x n with sufficiently large m satisfy the stable RSP with high probability, which further justifies the usage of the stable RSP an assumption in our analysis of IRLSp:
Proposition 3 (Gaussian + RSP). Suppose m 2 2k. Let A E Rm be a matrix with i.i.d. (0, 1) entries. Let 5 (0,1) and r7 E (0, 1) be fixed constants. if holds that
m +1 > 2k In(em/k)) (1.67 1,67+-71
(m n)²
181n((2550-1)
(10)
2k In (akkm//m/)))
then A satisfies the (ks 17) -ttblle RSP with probability at least 1
The assumption
2 2k of Proposition 3 is necessary, it is not hard to prove that, if m - n < 2k and Ax y is k-sparse, then the fg-minimization problem minzeR= || yloo has multiple solutions. With this assumption, (10) roughly becomes constant c if m n is large enough; thus, ignoring logarithmic and constant factors, condition (10) becomes m n 2 (k)), which nearly optimal, as the condition m 2 2k is necessary.
2 ck log(em/k) for some
Proposition 2 follows directly from [28, Theorem 4.9], so we omit its proof. The reader might also find Proposition 3 corresponds to [28, Corollary 9.34]. This corollary assumes D of (4) has i.i.d. N(0, 1) entries, from which is not immediate what the distribution of A is, thus it does not imply Proposition 3 directly; this is why provide complete proof for Proposition in Appendix ??.
3.2 Global Linear Convergence of IRLS
With the background provided Section 2, we are now ready to state the first main result:
Theorem 1 (Global Linear Convergence). Suppose A m2mxn obeys the (kt, n))-table RSP with ? E (0, 3/4) (cf. Definition 1). Lei Ax" y be k-sparse. If the smoothing parameter (!) is updated
as per (7) with a = k during the execution of the IRLS, Algorithm 1, then it holds for any 1 2 1 that
(1--5(1-- (3-44)) 294mm
Am (141) y
||
y11
3 3111 yll
(11)
meaning that the IRLS Algorithm converges linearly and globally in objective value.
Discussion. The condition 7 E (0,3/4) is better understood via Proposition 3, which asserts that this condition holds with high probability if A has ii.d. Gaussian entries and if m is large enough
5

---

(compared to n. k). Moreover, if the (k, -tabbe RSP holds, Proposition 2 implies that x* is a global minimizer of (1) with p = The strength of Theorem is in that guarantees linear convergence for any initialization, starting at the first iteration. However, the price to pay is that the convergence rate seems conservative with constant 1/294, and depends on the number m of samples: (11) predicts that 1 converges with accuracy 5 in (m log(1/8)) iterations and m can be large. We argue that such dependency is not an artifact of our analysis, because similar dependencies have been established for IRLS variants for fp-regression (1) with p > 2 [2, Theorem 3.1]. for interior point methods in convex optimization with m inequality constraints (184, Chapter 111, [55, Chapter 5), and for SGD in smooth & strongly convex optimization with sum of function components [85, Theorem 2.1 Moreover, an adversarial initialization for which the rate empirically depends on m was pointed out for IRLS for basis pursuit in [27]. That being said, we conjecture that the linear dependency
in (1I) can be improved to vm log m for some different choice of eO;; see [2] and [84, Sections 11.5 and 11 11.8.2] for why this conjecture makes sense. Furthermore, Figure la empirically depicts that, with the least-squares initialization, i.e., with w; (0) 1 for all TRLS) typically converges to x²' in 30 iterations. This hints to the possibility of alleviating the dependency on m. by analyzing the least-squares initialization, a challenging task which we leave to future work.
Comparison to [17]. To our knowledge, [17] is the only paper that claimed the global linear convergence of an IRLS variant with e(t) updated as per (5) and with 1. In particular, [17] sets c(t+)) + ()( whenever [|æ((++1] - 201 < 2,3c(t); see Figure 1a for the convergence of IRLS under this rule. The update rule e((+1) + (e(() of [17] has the advantage of being faster to compute than (7). In spite of this, and of their beautiful proof idea, their result has several disadvantages compared to Theorem 1. First, their proof involves sophisticated probabilistic arguments, and is tailored towards a feature matrix that satisfies strong concentration properties. Our analysis, on the other hand, uses the (k, -tt-ba RSP, which is close to a necessary and sufficient condition for the success of & é-rregression (ef. Proposition 2) and expected to hold for a much larger range of matrices than just (sub-) Sub-)Gaussiian matrices (see [86, 87] for related results in the context of basis pursuit). Second, some statements in their proof are inaccurate (e.g.. their Lemma does not hold for k < n, the last paragraph in their proving Lemma 8 is not rigorous). Third, their update rule incurs two hyper-parameters, (0) and 33 while Theorem 1 is based on the decreasing rule (7) that is adaptive to data, and the only dyper-parametee a is easier to set (Figure Ic).
Connection to [27]. A related global linear convergence was established for the IRLS method for basis pursuit in [27, Theorem 3.2]. Our proof is inspired by [27]. but also improves on its strategy. For example, their Theorem 3.2 involves two parameters for the stable mullspace property, while we only have a single RSP parameter 7), simplifying matters. Also, we obtain a constant of 1 1/294 in Theorem 1, which is better than the value of 1/768 for the respective constant in [27, Theorem 3.2].
It is important note that the proofs of [17, 27] and ours heavily rely certain rules for decreasing
the smoothing parameter e(t); all these proofs of global linear convergence for IRLS would break down if (()) were fixed. While the main theorems of both [17] and [27) are limited to the p = 1 case, we present our Theorem 2 for p E (0, 1) next.
3.3 Local Superlinear Convergence of IRLSp
The f, en-regression problem (1) with p E (0, 1) is more challenging than the case p = due to the lack of convexity. However, here we show that at least locally, IRLSp converges with a superlinear rate of order 2 p. Moreover, since it is valid to run IRLSp even with p = 0, we carefully design the proof such that the following result holds not only for p E (0, 1]. but also for p O, in which case IRLSp can be interpreted as an algorithm minimizing a oflooorlo objective (see Appendix ??).
Theorem 2 (Local Superlinear Convergence). Run IRLS, with p E [0, 1] and update e!!)) by (7) with = k. Assume A satisfies the (ks -7)-sable RSP (Definition 1). Let c E (0, 1) be sufficiently small constant such that 2ct pn(n + 1) < (1-c)2-P, -)- P. Let Ax - y be k-sparse with support St. Define
# (m + 1)(1 - c)p-2 2, min la x' yil " EES-
(12)
If the initialization (1) is in a neighhorhood of x=, in the sense that
Azx Ax ||, S c- min EE
aa x= Vols
(13)

---

then IRLSp achieves the following superlinear convergence rate of order 2 foror every 1:
|| AA' p ((IAx')))) - -x= 2
Axxl) (t) Ax'1l
(14)
In particular, for E [0, 1), the supertinear rate (14) implies
(x((++)- Ax" A1 - (uV(--) ||Axx()) Ax=)l1) (2-p)--1 Axx())
Ax=
(15)
Example 1 (Quadratic versus linear rates). We illustrate the practical benefit of a superlinear
(quadratic) convergence rate with the following simple calculation. Suppose that the incquality in (13)
is barely fulfilled such that in the case of p = 0, x1 satisfies μ. || Axx 1la S 0.9999999.
Then (15) implies that after 30 iterations, the residual error is far below numerical precision already
since ((1)) Ax ||2 S 0.99999992
||Ax(()) Ax"I| R 10-47, II || Ax ||1- On
the other hand, for 1, linear convergence factor of S 0.9999999 is only able to guarantee a decay of the order 1 < 0.999999930 xs 0.999997 after 30 iterations.
Discussion. Even though IRLS, with p E (0, 1) is tailored to Ey-regression, we measure the progress of I in l moormm e.g., we provide upper bounds on ||Az((+1) Ax [|2: this is because doing so avoids the use of Hölder's inequality, which allows us to give tighter results that improve on [26] (see the next paragraph). The second point that deserves discussion is the local convergence neighhorhood defined by the right-hand side of (13). Note that, we typically assume that every outlier sample (a, 84) (with i E 5") would result in a relatively large residual at x", say aa x - y/1 > O. Hence the minimization term of the right-hand side in (13) is well-behaved. Next, one might ask how to obtain such an initialization x! that satisfies (13). One might run IRLS1 to produce such an initialization; note though that the upper bound of (13) can not be computed, so one could not detect when to switch from IRLSI to SRLSp (p E [0. 1)). Or alternatively, one could run IRLSp with least-squares initialization and count on empirical global (super)linear convergence of TRLS,; the latter is what we did in the experiments and is what we recommend. A final and important remark is that f re-resression (I) with E (0, is in general NP-hard [28, Exercise 2.10]. but this does not contradict the theoretical local superlinear convergence of Theorem 2 and the empirical global convergence in Figures 1b and 2b, and this does not mean that IRLS, solves an NP-hard problem in polynomial time. The catch is that operate under the assumption that x leads to k-sparse residual and is a unique global minimizer of (1) (cf. Proposition 2). With this assumption and small enough k/ k/m, re-regression is tractable and can be solved via IRLS, to high accuracy (Figure lb).
Connection to [26]. [26, Theorem 7.9] proves the local superlinear convergence of IRLS with update rule (6) for b--basis pursuit (4), which motivates Theorem 2. However, [26, Theorem 7.9] requires c of (13) to satisfy = O(1mm3/p
which means that, as - 0 or m + oc, the eligible value of c quickly becomes vanishingly small, and thus the radius of local convergence (13) diminishes greatly, In contrast, in our condition 2c- 217-717(77 + 1) < (1 -
(1-c)2 c)2-P, c has no direct dependency on m and is well-behaved even if p + O. The above difference is due partly to our improved proof strategy, and partly to the different update rules of the smoothing parameter (f), Figure 1b showed that rule (6) of [26] does not work well for small p. in contrast to rule (7) that we use; a possible explanation for this phenomenon is that rule (6) decreases e(t) too fast, yielding a poor initialization for the next iteration of weighted least-squares (similarly the interior point method, cf. [84]). Finally, [26, Theorem 7.9] does not hold for p = 0, for which Theorem 2 still holds and suggests a local quadratic rate.
Applications and Experiments
We now explore the performance of IRLSp in three different applications, real phase retrieval (Section 4.1), linear regression without correspondences (Section 4.2), and face restoration (Section 4.3). Note that the first two applications are special examples of the recent algebraic-geomenn framework called homomorphic sensing [88], [89, Section 4), [90, Section 1.2]. In Section 4.4, we examine the behavior of IRLSp for different values of p and under noise.
4.1 Real Phase Retrieval
In real phase retrieval [56, 57], given measurement matrix A [ai,.... aml [T and y [y1..... Um]] where y |a, x"). and we need to find either x" -x". This problem
7

---

10
Relative Enror
10
Relative Error
Time (Seconds)
2
Kaczmarz (67.
10 Q
PlaseLamp [68]
TRI.S
Coordinate Def [70]
Truncated Wininger Flow (69)
Li Gurobi
IRLSp
0.1
18181
PDLP [71]
121So.1
Subgradiont Descent (78)
0.02
10
10 15
0.005
10
30
50
70
90
1000 3000 5000 7006 9000
1000
3000 5000 7000 9000
(16)
m
772
(8) 12
200
(b) n. = 50, k/m 0.5
(c) n = 50, k/m = 0.5
399
Figure 2: Figure 2a: Relative error min min xx" 2| [ + x: ||2//||==* of the methods that produce estimates è for real phase retrieval. Figures 2b-2c: Relative errors || x æ*112/ |122/|| ||2= l2 and running times for linear regression without correspondences. IRLS, run for most 50 iterations. 50 trials.
relative of the complex phase retrieval problem [65], where all data y and A, well as the ground-truth x", are complex-valued, and which has applications in X-ray crystallography [91].
Here we show that real phase retrieval can be solved via f.p regresssion (1), Consider the index sets
1+ = (i: a. x= > 0}: r = (i: a x < <0.
(16)
We might assume a xx * 0 for every i without loss of generality. Then we know that either Ax" -y y has its Eg norm equal to m - 1/- |. or A(-x") 30--- has its Eg norm equal to m - - With such sparsity patterns on these two residuals, relaxation of mo-minimization, to recover x or -x=, whichever corresponds to a sparser residual. Here, we essentially treat one of the two clusters defined by (16) inliers, and the other as outliers. This robust regression point of view on real phase retrieval seems to be known by experts [92, Section 1.2], but we have not found a paper that actually proposes to solve (1) for real phase retrieval.
minimize (I), which serves as a (non-convex)
Here we argue with Figure that, solving (1) (by IRLS0.11 for real phase retrieval can be bene- ficial. Indeed, theoretically. IRLS0.1 converges to +x* locally but superlinearly (as corollary of Theorem 2). Empirically, Figure 2a shows that IRLS succeeds in recovering +=² by using only m 2n - 1 samples, which is exactly the theoretical minimum for real phase retrieval [65, Propo- sition 2.5]. Figure 2a also shows that multiple other state-of-the-art methods" [57, 67, 68, 69, 70] fail to recover =x=t in this extreme situation, even though some of them have nearly optimal sample complexity, which requires roughly O(n) samples up to a logarithmic factor to succeed (caution: nearly optimal & optimal). As suggested by Proposition 3, for [+| | fixed, minimizing (1) with Gaussian matrix A enjoys m = (n)) sample complexity, up to also a logarithmic factor.
Despite these advantages, we also emphasize that solving (1) (via IRLSp.1) for real phase retrieval has an inherent limit: The performance depends on the number | of positive signs and the number m of samples; we might call |/+l/m the inlier rate or outlier rate. At n m 200/399, we see that IRLS0.1 succeeds | < 70, and if |+| 2 399 70 by symmetry, but it fails if [ and | get closer (Figure 2a). For IRLSo. to successfully handle the case |/+ [| we need more samples, empirically, say m 2 5n. It is an interesting future direction to design an algorithm that can handle the case of minimum samples (m 2n 1) and balanced data ( (-)).
4.2 Linear Regression without Correspondences
Compared to (real) phase retrieval, the problem of linear regression without correspondences is also important, but less developed subject; see [58, 94, 59, 60, 61, 95, 96, 97] for recent advances. In this problem are given y IR." and A m and we need solve the equations
y IIA
(17)
for some unknown permutation matrix II there exists solution (II", x²). Solving (17) is in general ha-hard for > 1 [94, 98]. However, if
R mmmx
xm and vector æ E R.", under the assumption that
'Reviewing these algorithus is beyond the scope of this paper. Some of them are designed for complex phase retrieval but is applicable to the real case. We used the implementations from PhasePack [93].
8

---

Relative Eiror
0.37
Least-Son Squares Least-Squares RLSp
Ground-Truth
Measurements
fi
Least-Squares"
Least-Squares
IRLSo.1
0.3
0.1
10%
30%
50%
70%
k/m
(a) Quantitative results
(b) Qualitative results (randomly selected images)
Figure 3: Recover a face image of the Extended Yale dataset [74] that corrupted by salt & pepper noise [63]. Figure 3a: relative error as function of the amount of random salt & pepper noise or sparsity level, averaged over trials, faces, and all individuals. Figure 3b: qualitative results.
we further assume that II" permutes most & of A, then know that the residual Ax" is k-sparse. This is an insight of [59], where it was proposed to estimate x' by solving (1) with = 1.
We show that IRLSp with p E (0, is more suitable than several baselines for solving this robust regression problem reduced from linear regression without correspondences. As baselines, we use the PDLP [71] and Gurobi [72] solvers to solve (1) with p 1 linear program, and also use subgradient descent method implemented the FOM toolbox of Beck & Guttmann Beck [73]. Figure 2b shows that TRLSSp.1 is the most accurate, and Figure 2c shows that I TRIS0.1 is 30 times faster than Gurobi and subgradient descent, and more than 100 times faster than PDLP. Interestingly. IRLS1 very competitive in terms accuracy (Figure 2b). Moreover, IRLS; is slower than IRLS as IRL converges faster and thus terminates earlier than IRIS1; we terminate IRLS, whenever the objective stops decreasing (up to a tolerance 10 10-15). Finally, see also [45, 48, 49, 50, 51, 52] where their IRLS variants outperform different set of baselines (for different problems).
4.3 Face Restoration from Sparsely Corrupted Measurements
Here we consider a simple face restoration experiment on the Extended Yale B dataset [74], downsam- pled as per [99]. This dataset contains the face images of 38 individuals under about R 60 different illuminations. Let f iff..... fi R is (vectorized) face image with low-rank (which is true under the Lambertian reflectance assumption [100]). thus each f. can be approximately represented
fu+1] be the matrix of faces of the same individual, where each
2, 016 pixels. We assume that F is approximately
linear combination of other faces of the individual, i.e., fi Faan for some x. E R2, where F; E mxxx is the same as F, except with the i-th column excluded.
In our experiments of Figure 3, we corrupt f. by random salt & pepper noise [63] and obtain the noisy face f. as measurement; so F: - f: is approximately sparse. We consider three methods to estimate n and thus to find estimate f. of fę: (1) Least-Squares that solves (1) with A F, y = f., and p = 2, (2) Least-Squares' that solves (1) with A F, y f:r and p = 2, used as golden baseline, and (3) TRLSS.11 which solves (1) with A F, and y f.. Figure 3a shows that TRISS..1 is more robust to salt & pepper noise, whose amount is k/m, than Least-Squares, and has lower error lff. fill2/11 fllz//lfillz (recall though that Least-Squares is statistically optimal for Gaussian noise). Figure 3b shows some images randomly chosen from the dataset; observe that Least-Squares tends "average" the illumination in all images F, while IRLS0.1 delivers faithful restoration.
4.4 The Choice of and Performance of IRLS, Under Noise
Here we outline two experiments that illustrate two different aspects IRLS, empirically.
In the first experiment, visualized in Figure 4a, we compare the decay of the relative errors of the iterates of IRLS, for different choices of p the case that residual is exactly k-sparse. We observe that TRLS exhibits a linear error decay the one hand and a superlinear decay for IRLS, with = 0.5 and p = that accelerates as p decreases towards confirming the rates predicted by Theorem and Theorem 2. (p 0.5, 0.1 are rather arbitrary choices in our experiments.)

---

IE-1
Relative Error
Relative Error
10-9
IRLS0s I8LS0
IE-2
10
IRLS1 TRt:So.s
IRLSua
10-16
5E-4
10
Number of Iterations
Number of Iterations
(a) The exact s-spas case
(b) The approximate k--parse case
Figure 4: Empirical convergence rates of I for different values of p (m 1000, 10, 200, averaged over 20 trials). Figure 4a: Ax* y is exactly k-sparse. Figure 4b: Ax" y is approximately k-sparse with 1% noise (see Section ?? for data generation).
In another experiment, visualized in Figure 4b, explore the robustness of IRLSp if the model assumption of k-sparse residuals is not satisfied, but only approximately. particular, the input vector R. of IRLSp is now such that its restriction on S with |SS k has i.i.d. A/0, 1) entries (corresponding to sparse corruptions), but its other coordinates are distributed as independent Gaussian variables with mean a: x and variance 0.01, where i (S*)² (corresponding to dense noise); see also the code provided in Appendix ??. We observe that also for this case of only approximately k-sparse residuals, IRLSp performs well, with faster and more accurate convergence for p 0.1 and = 0.5 compared to p = 1. We note that, strictly speaking, Theorem and Theorem 2 do not apply directly to this case of approximately k-sparse residuals, but it is not hard generalize to this case; see [27, Theorem A.1] for a similar result for IRLS applied to basis pursuit.
Discussion and Future Work
In this work, we established novel global linear and local superlinear rates of IRLSp for ep-regression (1) under the assumption of the stable range space property. Furthermore, we explored several applications for which IRLS, exhibits state-of-the-art results.
There are several directions that deserve further investigation. Theoretically, the Gaussian assumption of Proposition might be weakened; similar results would hold for sub-Gaussian distributions [86, 87]. Note also that, in experiments (Figures and 2), IRLSp works well even when m n - 2k. This leaves the question of whether the constant or logarithmic factors of condition (10) are too stringent for guarantecing the stable RSP to hold; can this be improved? Finally, it is tempting to carry out a global rate analysis of IRLS, with p (0, 1) for rerrgssssion; can one prove global (linear) rates under the stable RSP or other assumption?
Algorithmically, designing inexact solvers for the inner (weighted) least-squares problem (2) based on Krylov-subspace methods [101, 102] is expected to further accelerate our current IRLS implementa- tion (our current implementation is attached in Appendix ??). Also, our empirical experiments suggest that TRLS, (in particular p E (0, 1)) is worth being adopted and applied to many other outlier-robust estimation tasks beyond robust regression [45, 48, 49, 50, 51. 52, 103, 104, 105, 106, 107, 108].
Acknowledgments and Disclosure of Funding
This work was supported by grants NSF 1704458, NSF 1934979, NSF-IIS-1837991, ONR MURI 503405-78051, and the NSF-Simons Collaboration on the Mathematical Foundations of Deep Learning (NSF grant 1031.
References
[1] A. El Alaoui, X. Cheng, A. Ramdas. M. J. Wainwright, and M. I. Jordan. "Asymptotic behavior of Ey based Laplacian regularization in semi-supervised learning." in Conference Learning Theory.
10

---

pp. 879-906, 2016.
[2] D. Adil, R. Peng, and S. Sachdeva, "Fast, provably convergent IRLS algorithm for p-norm linear regression." Advances in Neural Information Processing Systems, vol. 32, 2019.
[3] D. Slepcev and M. Thorpe. "Analysis of p-Laplacian regularization in semisupervised learning," SIAM Journal Mathematical Analysis, vol. 51, no. pp. 2085-2120, 2019.
[4] Bubeck, M. B. Cohen, Y. T. Lee, and Y. Li, "An homotopy method for l, regression provably beyond self-concordance and in input-sparsity time," in ACM Symposium on Theory of Computing, pp. 1130-1137, 2018.
[5] D. Adil. Kyng, R. Peng, and Sachdeva, "Iterative refinement for l--norm regression." in ACM-SIAM Symposium on Discrete Algorithms, pp. 1405-1424, 2019.
[6] A. Jambulapati, Y. P. Liu, and A. Sidford, "Improved iteration complexities for overconstrained p-norm regression," " tech. rep., arXiv:21 (1ii::2111.0144::5 [cs.DS] 2021.
[7] J. Boscovich, "De litteraria expeditione per pontificiam ditionem, et synopsis amplioris operis, ac habentur plura ejus ex exemplaria etiam sensorum impessa," Bononiensi Scientiarum et Artum Instuto Atque Academia Commentarti, vol. 4, pp. 353-396, 1757.
[8] L. Plackett, "Studies in the History of Probability and Statisties. XXIX: The discovery of the method of least squares," Biometrika, vol. 59, pp.
08 1972.
[9] J. Wright and Y. Ma, Righ-Diimnsional Data Analysis with Low-Dinensional Models: Principles, Computation, and Applications. Cambridge University Press, 2020.
[10] P. J. Huber, "Robust estimation a location parameter," The Annals of Mathematical Statistics, vol. 35, no. 1, pp. 73-101, 1964.
[11] E. Candes, M. Rudelson. T. Tao, and R. Vershynin, "Error correction via linear programming." in IEEE Symposium on Foundations of Computer Science, pp. 668-681,22 2005.
[12] K. Bhatia, P. Jain, and P. Kar, "Robust regression via hard thresholding," Advances in Neural Information Processing Systems, vol. 28, 2015.
[13] K. Bhatia, P. Jain, Kamalaruban, and Kar, "Consistent robust regression." Advances in Neural Information Processing Systems, vol. 30, 2017.
[14] J. Yang. Y.-L. Chow, C. Ré, and M. W. Mahoney, "Weighted SGD for εp regression with randomized preconditioning." The Journal of Machine Learning Rexearch, vol. 18, no. I, pp. 7811-7853, 2017.
[15] D. Durfee, K. A. Lai, and S. Sawlani, "f, regression using lewis weights preconditioning and stochastic
gradient descent," in Conference On Learning Theory, pp. 1626 16666-166,2660 2018.
[16] A. S. Suggala, K. Bhatia. P. Ravikumar, and P. Jain, "Adaptive hard thresholding for near-optimal consistent robust regression," in Conference on Learning Theory, pp. 2892-2897, 2019.
[17) B. Mukhoty, G. Gopakumar, P. Jain, and P. Kar, "Globally-converger iteratively reweighted least squares for robust regression problems," in International Conference on Artificial Intelligence and Statistics,
pp. 313-322, 2019.
[18] S. Pesme and N. Flammarion, "Online robust regression via sgd on the f loss," Advances in Neural Information Processing Systems, vol. 33, pp. 2540-2552, 2020.
[19] X. Chen and M. Derezinski, "Query complexity of least absolute deviation regression via robust uniform convergence," Conference on Learning Theory, pp. 1144-1179, 2021.
[20] A. Parulekar, A. Parulekar, and F. Price, "11 regression with Lewis weights subsampling." tech. rep., ar arXiv:2105.00 100.393333333/ [es.LG)) 2021.
[21] Candes and T. Tao, "Decoding by linear programming." IEEE Transactions on Information Theory, 2005.
vol. 51, no. 12, pp.
[22] F. J. Candes, J. K. Romberg, and T. Tao, "Stable signal recovery from incomplete and inaccurate measurements," Communications Pure and Applied Mathematics, vol. 59, no. 8, pp. 1207-1223, 2006.
[23] D. L. Donoho, "For most large underdetermined systems of linear equations the minimal & fi-norm solution is also the sparsest solution," Communications on Pure and Applied Mathematics, vol. 59, no. 6,
pp. 797- 829, 2006.
[24] A. Cohen, W. Dahmen, and R. DeVore, "Compressed sensing and best k-term approximation." Journal of the American Mathematical Society, vol. 22, no. 1, pp. 211-231, 2009.
[25] T. Blumensath and M. E. Davies, "Iterative hard thresholding for compressed sensing," Applied and Computational Harmonic Analysis, vol. 27, no. 3, pp. 265-274, 2009.
[26] 1. Daubechies, R. DeVore, M. Forasier, and C. S. Güntürk, "Iteratively reweighted least squares mini- mization for sparse recovery," Communications on Pure and Applied Mathematics, vol. 63, no. 1. pp. 1-38, 2010.
11

---

[27] C. Kümmerle, C. Mayrink Verdun, and D. Stöger, "Iterativlyy reweighted least squares for basis pursuit with global linear convergence rate," Advances in Neural Information Processing Systems, 2021.
[28] S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive Sensing. Springer New York, 2013.
[29] R. Chartrand, "Nonconvex compressed sensing and error correction." " in IEEE International Conference on Acoustics, Speech and Signal Processing, vol. 3. pp. III-889, 2007.
[30] R. Chartrand, "Exact reconstruction of sparse signals via nonconvex minimization," IEEE Signal Process- ing Letters, vol. 14, no. 10, pp. 707-710, 2007.
[31] R. Chartrand and Staneva, "Restricted isometry properties and nonconvex compressive sensing." " Inverse Problems, vol. 24, no. 3, p. 035020, 2008.
[32] S. Foucart and M.-J. Lai, "Sparsest solutions of underdeterminnd linear systems via E--minimization for
0 < q I 1," Applied and Computational Harmonic Analysis, vol. 26, no. 3, pp. 395-407, 2009.
[33] M. Wang. W. Xu, and A. Tang, "On the performance sparse recovery via e--minimization (0 p 1)," IEEF Transactions on Information Theory, vol. 57, no. 11, pp. 7255-7278, 2011.
[34] Q. Sun, "Recovery of sparsest signals via ey-minimization." Applied and Computational Harmonic Analysis, vol. 32, no. 3. pp. 329-341 2012.
[35] D. Ba, B. Bahadi, P. L.. Purdon, and E. N. Brown, "Convergence and stability of iteratively re-weighted least squares algorithms," IEFE Transactions on Signal Processing, vol. 62, no. 1, pp. 183-195, 2013.
[36] Zhou, Q. Zhang, and A. M.-C. So, "fr. regularization: Error bounds and convergence rate analysis of first-order methods," " in International Conference on Machine Learning, vol. 37, pp. 1501- 1510, 2015.
[37] L. Zheng, A. Maleki, H. Weng. X. Wang. and T. Long. "Does fy-minimization outperform 41 minimization"." IEEE Transactions on Information Theory, vol. 63, no. 11, pp.
2017.
[38] G. Marjanovic and V. Solo, "On & optimization and matrix complction," IEEE Transactions on Signal Processing, vol. 60, no. 11, pp. 5714-5724 2012.
[39] K. Mohan and M. Fazel, "Iterative reweighted algorithms for matrix rank minimization." " The Journal of Machine Learning Research, vol. 13, no. 1, pp.
2012.
[40] F. Nie, H. Huang, and C. Ding. "Low-rank matrix recovery via efficient Schatten p-norm minimization," in AAAI Conference on Artificial Intelligence, 2012.
[41] C. Kümmerle and J. Sigl, "Harmonic mean iteratively reweighted least squares for low-rank matrix recovery," The Journal of Machine Learning Research, vol. 19, no. 1. pp. 1815-1865 2018
[42] P. Giampouras, R. Vidal, A. Rontogiannis, and B. Haeffele, "A novel variational form of the Schatten-p quasi-norm." Advances in Neural Information Processing Systems, vol. 33, pp. 21453-21463, 2020.
[43] C. Kümmerle and C. Mayrink Verdun, "A scalable second order method for ill-conditioned matrix completion from few samples," in International Conference on Machine Learning, pp. 5872-5883, 2021.
[44] Lerman, M. B. McCoy, J. A. Tropp. and Zhang, "Robust computation of linear models by convex relaxation," " Foundations of Computational Mathematics, vol. 15, no. 2. pp. 363-410, 2015.
[45] M. C. Tsakiris and R. Vidal, "Dual principal component pursuit," " Journal of Machine Learning Research, vol. 19, 18, pp. 1-50, 2018.
[46] Z. Zhu, Y. Wang. D. Robinson, D. Naiman, R. Vidal, and M. C. Tsakiris, "Dual principal component pursuit: Improved analysis and efficient algorithms," in Advances in Neural Information Processing Systems, 2018.
[47] Ding, Z. Zhu, T. Ding, Y. Yang, R. Vidal, M. C. Tsakiris, and D. Robinson, "Noisy dual principal component pursuit," in International Conference on Machine Learning, pp. 1617-1625, 2019.
[48] K. Aftah and R. Hartley, "Convergence of iteratively re-weighted least squares to robust M-estimators." in TEEE Winter Conference on Applications of Computer Vision, pp. 480-487, 2015.
[49] W. Dong, X.jj. Wu, and J. Kittler, "Sparse subspace clustering via smoothed E, minimization," " Pattern Recognition Letters, vol. 125, pp. 206-211, 1, 2019.
[50] D. Iwata, M. Waechter, W.-Y. Lin, and Y. Matsushita, "An analysis of sketched IRL. for accelerated sparse residual regression." in European Conference on Computer Vision, pp. 609-66966006 626. 2020.
[51] T. Ding, Y. Yang, Z. Zhu, D. P. Robinson, R. Vidal, L. Kneip, and M. C. Tsaaiiiss "Robust homography estimation via dual principal component pursuit," in IEEE Conference on Computer Vision and Pattern Recognition, pp. 6080-6089, 2020.
[52] H. Yang, P. Antonante, V. Tzoumas, and L. Carlone, "Graduated non-convexity for robust spatial perception: From non- minimalim solvers to global outlier rejection," IEEE Robotics and Automation Letters. vol. 5, no. 2. pp. 1127-1134, 2020.
12

---

[53] Y.-B. Zhao and D. Li, "Reweighted fi-minimization for sparse solutions to underdetermined linear systems," SIAM Journal Optimization, vol. 22, 3. pp. 1065-1088. 2012.
[54] A. Beck, First- First-Order Methods in Optimization. Society for Industrial and Applicd Mathematics, 2017.
[55] Nesterov, Lectures on Convex Optimization. Springer, 2018.
[56] Y. Chen, Y. Chi, J. Fan, and C. Ma. "Gradient descent with random initialization: fast global convergence for nonconvex phase retrieval," Mathematical Progranmingg vol. 176, no. 1, pp. 5-37, 2019.
[57] Y. S. Tan and R. Vershynin, "Phase retrieval via randomized Kaczmarz: Theoretical guarantees." Informa- tion and Inference: A Journal of the IMA, vol. no. 1. pp. 97-123, 2019.
[58] J. Unnikrishnan, S. Haghighatshoar, and M. Vetterli, "Unlabeled sensing with random lincar measure- ments," IEEE Transactions on Information Theory, vol. 64, no. 5, pp. 3237-3253, 2018.
[59] M. Slawski and E. Ben- Ben-David,"" "Linear regression with sparsely permuted data." Electronic Journal of Statistics, vol. 13, no. I, pp. 1-36, 2019.
[60] M. C. Tsakiris, L. Peng, A. Conca, L. Kneip, Y. Shi, and H. Choi, "An algebraic-geometric approach for linear regression without correspondences," IEEE Transactions on Information Theory, vol. 66, no. 8, pp. 5130-5144, 2020.
[61] L. Peng and M. C. Tsakiris, "Linear regression without correspondences via concave minimization," JEEE Signal Processing Letters, vol. 27, pp. 1580-1584 2020.
[62] J. Wright, A. Y. Yang. A. Ganesh, S. Sustry, and Y. Ma, "Robust face recognition via sparse represen- tation." IEEE Transactions on Puttern Analysis and Machine Intelligence, vol. 31, no. 2, pp. 210-227. 2008.
[63] C. Solomon and T. Breckon, Fundamentals of Digital Image Processing: A practical approach with examples in Matlab. John Wiley & Sons, 2011.
[64] Yao, L. Peng, and M. Tsakiris, "Unlabeled principal component analysis." Advances in Neural Information Processing Systems, 2021.
[65] R. Balan, P. Casazza, and D. Edidin, "On signal reconstruction without phase," Applied and Computational Harmonic Analysis, vol. 20, no. 3, pp. 345 356, 2006.
[66] T. Strohmer and R. Vershynin, "A randomized Kaczmarz algorithm with exponential convergence." Journal of Fourier Analysis and Applications, vol. 15, no. 2, pp. 262-278, 2009.
[67] K. Wei, "Solving systems of phaseless equations via Kaczmarz methods: A proof of concept study." Inverse Problems, vol. 31, no. 12, p. 125008, 2015.
[68] O. Dhifallah, C. Thrampoulidis, and Y. M. Lu, "Phase retrieval via linear programming: Fundamental lim- its and algorithmic improvements." in Allerton Conference on Communication, Control, and Computing. 1071-1077, 2017.
pp.
[69] Y. Chen and E. Candes, "Solving random quadratic systems of equations is nearly as easy as solving linear systems," " Advances in Neural Information Processing Systems, vol. 28, 2015.
[70] W.-J. Zeng and H.-C. So, "Coordinate descent algorithms for phase retricval," tech. rep., arXiv: [cs.IT], 2017.
[71] D. Applegate, M. Diaz, O. Hinder, H. Lu, M. Lubin, B. O'Donoghue, and W. Schudy, "Practical large- scale linear programming using primal-dua hybrid gradient," Advances in Neural Information Processing Systems, vol. 34, 2021.
[72] "Gurobi 9.5.1." https://www.gurbi..com
Date Accessed: May 10, 2022.
[73] A. Beck and N. Guttmann- Bartmann-Beck, "FOM a MATLAB toolbox first-order methods for solving convex optimization problems," Optimization Methods and Software, vol. 34, no. 1, pp. 172-193, 2019.
[74] A. S. Georghiades, P. N. Belhumeur, and D. J. Kriegman, "From few many: Illumination cone models for face recognition under variable lighting and pose," JEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 23, no. 6, pp.
2001.
[75] S. S. Chen, D. I., Donoho, and M. A. Saunders, "Atomic decomposition by hasis pursuit," SIAM Review, vol. 43, no. 1, pp. 129-159, 2001.
[76] G. Lerman and T. Maunu, "Fast, robust and non-convex subspace recovery." Information and Inference:
A Journal of the IMA, vol. 7, no. 2. pp. 277-336, 2018.
[77] Q. Qu, Z. Zhu, X. Li, M. C. Tsakiris, J. Wright, and R. Vidal, "Finding the sparsest vectors in subspace: Theory, algorithms, and applications," " tech. rep., arXiv:2001.066 006970 [cs.LG], 2020.
[78] T. Chan and P. Mulet, "On the convergence of the lagged diffusivvtty fixed point method in total variation image restoration." SIAM Journal Numerical Analysis, vol. 36, no. 2, pp. 354-367, 1999.
13

---

[79] A. Beck, "On the convergence of alternating minimization for convex programming with applications to iteratively reweighted least squares and decomposition schemes," SIAM Journal on Optimization, vol. 25, no. 1. pp. 185-209. 2015.
[80] A. Y. Aravkin, J. V. Burke, and D. He. "IRLS for sparse recovery revisited: Examples of failure and a remedy." tech. rep.. 190.0:191 [math.ST]. 2019.
[81] ]] X. Yang. J. Wang. and H. Wang. "Towards an efficient approach for the nonconvex by-ball projection: Al Algorithm and analysis," tech. rep., arXiv:2101 .0135vv5 [math.OC], 2021.
[82] A. Tillmann and M. E. Pfetsch, "The computational complexity of the restricted isometry property, the nullspace property, and related concepts in compressed sensing," IEEE Transactions on Information Theory, vol. 60, no. 2. pp. 1248- 1259, 2013.
[83] S. Arora and B. Barak, Computational Complexity: A Modern Approach.. Cambridge University Press, 2009.
[84] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex optimization. Cambridge University Press, 2004.
[85] Z. Allen-Zhu, "Katyusha: The first direct acceleration of stochastic gradient methods," The Journal of Machine Learning Research, vol. 18, no. 1, pp. 8194-8244, 2017.
[86] G. Lecué and S. Mendelson, "Sparse recovery under weak moment assumptions," Journal of the European Mathematical Society, vol. 19, no. 3, pp. 881-904, 2017.
[87] S. Dirksen, G. Lecué, and H. Rauhut, "On the gap between restricted isometry properties and sparse recovery conditions," " IEEE Transactions on Information Theory, vol. 64, no. 8. pp. 5478-5487, 2016.
[88] M. C. Tsakiris and L.. Peng, "Homomorphic sensing," in International Conference on Machine Learning, 2019.
[89] Peng, B. Wang, and M. Tsakiris, "Hommmorphic sensing: Sparsity and noise," in International Conference on Machine Learning, pp. 8464-8475, 2021.
[90] L. Peng and M. C. Tsakiris, "Homomorphic sensing of subspace arrangements." Applied and Computa- tional Harmonic Analysis, vol. 55, pp. 466-485, 2021.
[91] P. Grohs, S. Koppensteiner, and M. Rathmair, "Phase retrieval: Uniqueness and stability," SIAM Review, vol. 62, no. 2. pp. 301-350 2020.
[92] M. Huang and Y. Wang, "Linear convergence of randomized Kaczmarz method for solving complex- valued phaseless equations." tech. rep., arXiv:2109. 11811 [math.NA], 2021.
[93] Chandra, T. Goldstein, and C. Studer, "PhasePack: A phase retrieval library," International conference on Sampling Theory and Applications. pp. 1 1-5,2019.
[94] A. Pananjady, M. J. Wainwright, and T. A. Courtade, "Linear regression with shuffled data: Statistical and computational limits permutation recovery," IEEE Transactions on Information Theory, vol. 64, no. 5. pp. 3286-3300, 2018.
[95] M. Slawski, G. Diao, and E. Ben-David," "A pscudo-likclihood approach to linear regression with partially shuffled data," Journal of Computational and Graphical Statistics, vol. O, no. 0. pp. 1-31, 2021.
[96] Y. Xie, Y. Mao, S. Zuo, H. Xu, X. Ye, T. Zhao, and H. Zha, "A hypergradient approach to robust regression without correspondence," in International Conference on Learning Representations, 2021.
[97] Li, K. Fujiwara, F. Okura, and Matsushita, "Generalized shuffled lincar regression," in JEEE/CVF
International Conference on Computer Vision, pp. 6474-6483, 2021.
[98] D. Hsu, K. Shi, and X. Sun, "Linear regression without correspondence," in Advances in Neural Informa-
tion Processing Systems, 2017.
[99] C. You, D. Robinson, and R. Vidal, "Scalahle sparse subspace clustering by orthogonal matching pursuit,"
in TEEE Conference on Computer Vision and Pattern Recognition, pp- 3918-3927, 2016.
[100] R. Basri and D. W. Jacobs, "Lambertian reflectance and linear subspaces," IEEE Transactions on Pultern Analysis and Machine Intelligence, vol. 25, no. 2. pp. 218-233, 2003.
[101] M. Fornasier, S. Peter, H. Rauhut, and Worm, "Conjugate gradient acceleration iteratively re- weighted least squares methods," Computational Optimization and Applications, vol. 65, no. 1, pp. 205- 259, 2016.
[102] Gazzola, C. Meng, and J. G. Nagy. "Krylov methods for low-rank regularization," SIAM Journal on Matrix Analysis and Applications, vol. 41, no. 4. pp. 1477-1504, 2020.
[103] A. Chatterjee and V. M. Govindu, "Robust relative rotation averaging," JEEE Transactions Patter Analysis and Machine Intelligence, vol. 40, no. pp. 958 972, 2017.
[104] C. Sidhartha and V. M. Govindu, "It is all in the weights: Robust rotation averaging revisited," in International Conference 3D Vision, pp. 1134 1 1143, 2021.
14

---

[105] L. Peng. M. C. Tsakiris, and R. Vidal, "ARCS; Accurate rotation and correspondences search," in
IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 11153-11163, 2022.
[106] H. Lee and J. Civera, "HARA: A hierarchical approach for robust rotation averaging." in /EEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 15777-15786, 2022.
[107] L. Peng, M. Fazlyab, and R. Vidal, "'emidefinite relaxations of truncated least-squares in robust rotation search: Tight or not," in European Conference on Computer Vision. pp. 0-0, Springer. 2022.
[108] L. Carlone, "Estimation contracts for outlier-robust geometric pereeption." tech. rep., arXiv:2208. 10522 [stat.ML], 2022.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] Furthermore, we made a lot of efforts in comparing the theoretical results of prior works and position our contribution in the literature.
(b) Did you describe the limitations of your work? [Yes] See for example at Lines 268-274 on real phase retrieval. We also discussed several directions for improvements or future research at the end of the paper.
(c) Did you discuss any potential negative societal impacts of your work? [No]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full of assumptions of all theoretical results? [Yes] In particular, we used an entire subsection (Section 3.1) to justify the major assumption, the stable range space property.
(b) Did you include complete proofs of all theoretical results? [Yes] We omitted the proofs of Proposition 1 and 2 as they follow directly from existing results; see Line 83 and Line 155. All other theoretical results are proved in the appendix.
3. If you experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as URL)? [Yes] We attached a copy of our code in the supplemental material.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] In the appendix, we use an entire section to describe our experimental setup. We did not do that in the main paper in interest of space.
(c) Did you report error bars (e.g., with respect the random seed after running experi- ments multiple times)? ? [No] We did not report error bars for two reasons. The main reason is that some of our figures are quite dense, and adding error bars on them makes it harder to read. A slightly subjective reason is that we believe our current figures are informative enough, correctly deliver the message, and are fair for other methods.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [N/A]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [No] In the main paper we used the Extended
Yale B face dataset, which is free for research purposes.
(c) Did you include any new assets either in the supplemental material or as URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're asing/curating? [Yes]
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [No]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
15

---

(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
16


Reviews:
[{"review_id": "L7FWCATRJZX", "reviewer": "Reviewer_yi3J", "summary": "This work analyzes the famous IRLS algorithm with any choice of $p\\in [0,1]$ and the adaptively updated smoothing parameter $\\epsilon_t$. It is able to establish global linear convergence theory for $p=1$ with arbitrary initialization. Unlike previous works, it also establishes local superlinear convergence theory with sufficiently good initialization, for any $p\\in [0,1]$. Experiments confirm the advantages of IRLS with $p<1$.", "questions": "1. Would it be possible to analyze the radius of the convergence zone (RHS of eq. 13) under some probabilistic assumptions. For example, one may consider the entries of $A$, $x^*$ are i.i.d standard normal, $y = Ax^* +r^*$, where $r^*$ is k-sparse. The coordinates of nonzero elements of $r^*$ are randomly sampled, and its nonzero elements follow any smooth distribution with zero mean and finite second order moments. Would the RHS of (13) converge to 0 as $k=\\lambda m$, $n,m\\to \\infty$ where $\\lambda$ is a small constant? I suspect that it is true (since you take a min of a set of random variables that has a growing cardinality), but I would like the authors to confirm or disprove my conjecture.\n\n2. Following my point 1, if my conjecture is true, would it contradict your claim in line 66 that your radius of local convergence does not diminish as $m\\to \\infty$? It probably won't diminish if $k$ is fixed even when $m\\to \\infty$ (not $k=\\lambda m$), but the author should add some discussions about which assumption is more reasonable, and which one is more common in real practice. In some sense, $k$ should grow with $m$, but I welcome the author(s) to provide me with some counterexamples in real applications.\n\n3. As for the experiments, the performance of IRLS-0.1 is very impressive. Is $p=0.1$ typically (approximately) the best choice in all experiments? What if you choose $p=0.05$? Also, how the IRLS-1 and IRLS-0.1 are sensitive to the choice of update rules of smoothing parameters (as you introduced in section 2.2)? Is the current update rule of $\\epsilon$ critical to the success of the presented IRLS-0.1? It might be good to add some simple experiments in the supplementary material.\n\n\n", "limitations": "Yes.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. This work improved the analysis and the result of the previous global convergence theory for IRLS with $p=1$.\n\n2. The authors are able to show local superlinear convergence for any $p\\in [0,1]$, which is not done in previous works.\n\n3. The paper is very well-written. The previous works are reviewed in detail. The contributions are highlighted. The content is well-structured.\n\nWeaknesses:\n\n1. Although the stable RSP condition is satisfied for Wigner matrices, it does not seem to be easy to verify in general (one has to check all subsets of coordinates of size k). I understand that it was also adopted in some previous works, but it would be nice to have a condition that can be more easily verified in a general setting.\n\n2. The radius of the convergence zone (RHS of eq. 13) for the nonconvex IRLS depends on the min over the nonzero residuals. Thus, the radius could be small when k (the number of nonzero elements), n, m are all large (consider an extreme case where $k,n,m\\to \\infty$), unless one can show that the smallest nonzero element is sufficiently far away from zero.\n\n\nThis work has made a nice contribution to the theory of the iteratively reweighted least squares, although there is still room for improvement on the theory side.\n\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "8Pd7jCzr3dy", "reviewer": "Reviewer_qErq", "summary": "This paper studies the classical algorithm of IRLS for solving $\\ell_p$ regression for $p \\in (0,1]$. The algorithm updates between solving a weighted least squares problem and updating the weights, where the authors also use a dynamically updated smoothing (regularization). Under a stable range space property, which is satisfied by Gaussian matrices, the authors prove that IRLS with $p=1$ converges linearly. Furthermore, for $p<1$, the authors show that the algorithm converges at a superlinear rate provided that the algorithm has good local initialization.", "questions": "- What is the point of the $(p,k,\\eta)$ RSP when even for the case of $p<1$ you use the $(1,k,\\eta)$ RSP?\n- Is there a connection between the IRLS algorithm with $p=0$ and a Newton method, which is typical of exhibiting such quadratic rates?\n- Does Proposition 3 hold for sub-Gaussian rather than just Gaussian?\n- Since the regularization update as well as the reweighting scheme have been considered before, what are the “novel variants of IRLS” considered here? What other notions of IRLS with $p<1$ exist in the literature? Some examples are  \n    - Mohan, K., & Fazel, M. (2012). Iterative reweighted algorithms for matrix rank minimization. The Journal of Machine Learning Research, 13(1), 3441-3473.\n    - \"​Lerman, G., & Maunu, T. (2018). Fast, robust and non-convex subspace recovery. Information and Inference: A Journal of the IMA, 7(2), 277-336\n- What is the connection, if any, between RSP conditions and the notion of restricted isometry properties?\n\n", "limitations": "Limitations of the algorithm for $p<1$ are discussed in the example of phase retrieval. Also, the suboptimality with respect to $m$ is discussed in 177-179.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper is very clear and well-written.\n- The authors give an improved result for the convergence of IRLS with $p=1$ over [15]. In particular, the result of the authors works under a generic condition and uses a better smoothing update rule.\n- The result for superlinear convergence with $p<1$ is interesting, and the result extends to the case of $p=0$, which is the first time I have seen this in the literature.\n- The choice of $p<1$ for phase retrieval achieves the theoretical minimum number of samples.\n\nWeaknesses:\n- There is not a novel algorithmic twist developed here, and the work is gauged to purely developing theory. How realistic is the assumption in the linear regression without correspondences that only $k$ of the rows are permuted?\n- It is unclear how stringent the initialization condition in Theorem 2 is.\n- There is not a driving application for this method, and it rather serves as a general-purpose methodology.\n- There is no discussion of the choice of $p$ in practice. It seems like $p<1$ is beneficial. However, in the case of phase retrieval, it seems to do poorly in the balanced regime of $|I^+| = |I^-|$.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "2NgMJWq8fl9", "reviewer": "Reviewer_GXHV", "summary": "In this paper, the authors propose a new variant of iteratively reweighted least square (IRLS) algorithm and studied its theoretical property as well as the practical performance. The authors stated the stable range space property, which plays the fundamental role in the theoretical analysis of the new IRLS method. The authors showed that for the case of $p = 1$, the IRLS can achieve a global linear convergence rate. The authors also presented that for the case of $p \\in (0, 1)$, the IRLS reach a superlinear convergence rate of the order of $2 - p$. Both local and global convergence rates can match the state-of-art results. The authors also conduct numerical experiments on different applications. And the empirical results showed that the performance of the IRLS algorithm is competitive against other methods.", "questions": "One question for the structure of the paper.\n\nCould the authors put the discussion and future work on page 9 of the paper into one independent conclusions section? This can make the paper more clean and clear.", "limitations": "The authors have presented the limitations of this paper in the final page of the submission.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This paper has the following strengths in terms of the originality, significance, quality and clarity:\n\nOriginality and significance: as far as we know, this paper provides the one of the first novel IRLS algorithm that could reach linear rate for the global convergence and superlinear rate for the local convergence. This is the main and independent theoretical contribution of this paper.\n\nQuality and clarity: the authors provides detailed proof in the theoretical section and the consistent empirical results in the numerical experiments section. These empirical results show that the performance of the proposed IRLS method is competitive against other state-of-art algorithms in the application areas of real phase retrieval, regression without correspondences and robust face restoration. The new IRLS algorithm is specially good at handling the outliers, which is the key property in the problem of robust $l_p$ quasinorm regression. The paper is also organized with clear and clean structure and the language is well-organized so that the readers can easily understand the content of this submission.\n\nHowever, despite the above advantages, the paper has the following weaknesses.\n\n The theoretical results only hold under the assumptions that $m$ must be much larger than the $n$ and $k$. However, in the numerical experiments, the IRLS algorithm also performs well under the condition that $m$ is comparable larger than $k$ and $n$. Hence, there is some gap between the theoretical analysis and the numerical experiments. Also the authors should mention in the paper any theoretical results for the case of $m$ is not large enough, i.e.,  assumption (10) doesn't hold. Is there any significant results for this special case?\n\nAnother minor issue is that the authors should add an independent conclusion section after the numerical experiments/applications section for the summary, discussion and future research directions of the paper.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept