Improved $\ell_{p}$ Regression via Iteratively Reweighted Least Squares
Authors: Alina Ene
Ta Duy Nguyen
Adrian Vladu
Venue: ICLR 2026 Poster (2026)
Keywords: lp regression
optimization
Primary area: optimization
TL;DR: We introduce fast practical algorithms for solving $\ell_{p}$ regression, outperforming the best known practical algorithm and matching the state of the art theoretical guarantee.

Abstract: We introduce fast algorithms for solving $\ell_{p}$ regression problems using the iteratively reweighted least squares (IRLS) method. Our approach achieves state-of-the-art iteration complexity, outperforming the IRLS algorithm by Adil-Peng-Sachdeva (NeurIPS 2019) and matching the theoretical bounds established by the complex algorithm of Adil-Kyng-Peng-Sachdeva (SODA 2019, J. ACM 2024) via a simpler lightweight iterative scheme. This bridges the existing gap between theoretical and practical algorithms for $\ell_{p}$ regression. Our algorithms depart from prior approaches, using a primal-dual framework, in which the update rule can be naturally derived from an invariant maintained for the dual objective. Empirically, we show that our algorithms significantly outperform both the  IRLS algorithm by Adil-Peng-Sachdeva and MATLAB/CVX implementations.

Published as a conference paper at ICLR 2026
IMPROVED ep REGRESSION REWEIGHTED LEAST SQUARES
VIA ITERATIVELY
Alina Ene Department of Computer Science Boston University aene@bu. edu
Ta Duy Nguyen
Adrian Vladu CNRS & IRIF
Department of Computer Science Boston University taduy@bu. edu
Université Paris Cité v vladue Arif.fr
ABSTRACT
We introduce fast algorithms for solving ép regression problems using the iteratively reweighted least squares (IRLS) method. Our approach achieves state-of-the-ws iteration complexity, outperforming the IRLS algorithm by Adil-Peng-Sachdeww (NeurIPS 2019) and matching the theoretical bounds established by the complex algorithm of Adil-Kyng Kyng- Sacheeav (SODA 2019, ACM 2024) via a simpler lightweight iterative scheme. This bridges the existing gap between theoretical and practical algorithms for l regression. Our algorithms depart from prior approaches, using primal-duc framework, in which the update rule can be naturally derived from an invariant maintained for the dual objective. Empirically, we show that our algorithms significantly outperform both the IRLS algorithm by Adil-Peng- Sachdeva and MATLAB/CVX implementations.
INTRODUCTION
In this paper, we study the lp regression problem defined as follows. The input the problem is matrix A _dx", a vector E R2 that lies in the column span of A, and an accuracy parameter e. The goal is to approximately solve the problem min Az=b ||>]|p- i.e., find a solution x E IR" such that Ar and ||~||p- S (1 +>]l=" pl.. where x is an optimal solution to the problem, and Il-llp denotes the Ep norm. Solving lp regression for all values of p is fundamental problem in machine learning with numerous applications and has been studied in long line of research beyond the classical least squares regression with p 2. Ly- Lp-norm regression problems with general p arise in several areas, including supervised learning, graph clustering, and wireless networks. Examples of applications include fn-norm based algorithms in semi-supervised learning (Alaoui 2016 Liu and Gleich 2020), k-clustering with E, e,-norm (Huang and Vishnoil 2020), robust regression and robust clustering (eng and Mahoney 2013 Huang et al.])222))
For this general class of convex optimization problems, designing provably fast iterative algorithms to obtain high accuracy solutions with empirical efficiency is an important question. General convex programming methods such as interior point methods are usually slow in practice. In theory, Bubeck ct al. (2018) show that algorithms based on interior point methods cannot improve beyond O ((T) iterations for any p & {1,2, so). Breaking this barrier and finding iterative algorithms that faster than interior point methods both in theory and practice is the goal of this line work.
Recent developments have led new algorithmic approaches such as homotopy method (Bubeck et al., 2018), and an iterative refinement approach (Adil et al. 2019a b. 2024) for lp regression with p & (T, oo). We highlight the notable works by Adil al. (2019a 2024). On the one hand, the algorithm with the best known theoretical runtime is given by Adil et al. 2019a 2024 with
Forr simplicity the introduction, we assume that d (zx). In regime when n d, the IPM iteration complexity improves ((v)).

---

Published as a conference paper at ICLR 2026
O(p2n log (2)) calls to a linear system solver. This algorithm. however, relies on complex subroutines and includes theoretical choices for several hyperparameters. In practice, to obtain an efficient implementation, hyperparameters require tuning. Due to these reasons, this theoretical algorithm by Adil et al. (2019a 2024) does not provide a practical implementation. On the other hand, an algorithm known as p- o-IRLS by Adil et al. (2019b) has been shown to have significant speed up over standard solvers such CVX. This algorithm is implemented based
Iteratively Reweighted Least Squares (IRLS) method, which is general iterative framework for solving regression problems. The key element of an IRLS method is solving weighted least squares regression problem in each iteration. This is equivalent to solving linear system of the form minggR=::An diagonal matrix, which can be computed very efficiently in practice with the advance of numerical solvers. IRLS algorithms are favored in practice (Burrus) 2012). but designing IRLS algorithms with strong convergence guarantees is challenging. In particular, to obtain the efficiency, the algorithm by p-2
Ra. where R is
by log (2)) linear system Adil et al. (201 9b) sacrifices the theoretical guarantee, requiring Opp (y"n"" solves. This brings forth the question:
Can we design an algorithm that retains the empirical efficiency of an IRLS approach while achieving the state-of-the-ar theoretical runtime?
In this work, we give a positive answer to this question. We provide a new algorithmic framework for Pp regression based IRLS approach for all values of p (1, oc). We propose an algorithm that
uses Opp'n 2 log n linear system solves, matching the state-of-thh-arr theoretical algorithm
by Adil et al. (219a). and improving upon the guarantee of Op²n log (²)) for the p-IRLS
algorithm by Adil et al. (2019b). We experimentally compare our algorithm with the p-IRLS algorithm (Adil et al. 2019b) and CVX solvers, and we observe significant improvements in all instances.
1.1
OUR CONTRIBUTIONS
For the simplicity of the exposition, we study the lp regression problem in both low and high precision regimes for p 2 2.
B we show a simple reduction for the more general problem minz: Ar=b || Na: rllp to the form min x; År= ||>||p with the dependence of the runtime the number of rows of N instead of the dimension of x. We also show in Appendix Cla reduction for the case p to the case p >
Remark 1.1. In Appendix
In the low precision regime when the runtime dependence on e is poly (t). we have the following theorem.
Theorem 1.1. For any p 2 2, there is an iterative algoritium for the Ep regression problem minzeR": Az=b ||~:|pp that solves O (log log n + log (1/c)) subproblems, each of which makes
se2 4p by sy14 log )) calls solve a linear system of the form ADA o = b, where D is an arbitrary non non-negative diagonal matrix.
(((!) " + nor (2)
Remark 1.2. When p oc, each subproblem makes 0 + n. log(2)) calls to linear system solver.
Prior approaches for solving f regression problem in the low precision regime commonly use the Taylor expansion of ||=]| which then allows for deriving and bounding the updates. In contrast to this, our algorithm relies on primal-dual approach using the dual formulation of the squared objective min,: Ar=b ||=]| = min.: Ar=b ||=2 max, E(r) where bg is the dual norm of f// [|rll,
and {r) = A (r, 2²) The term E(r) is often referred to the energy. The high level idea of our approach is as follows. Starting with an initial solution r for the dual problem, we will increase the coordinates of much possible so that the increase in the energy E(r) relative the increase
|| | I1=) - for finding 2 such that ||| S min ||=11 ||=
Thhe original result is O (pn &o log
c. This translates to Opn # log
O Opp" f log (?)) for finding ã such that | S
(1 Ar=b ||=||, for (0) initialized to min |1//|2-
2

---

Published as a conference paper at ICLR 2026
of ||lll, is also sufficiently large, until we can obtain a (1 1 c) optimal dual solution and whereby recover approximately optimal primal solution. This template is close to the approach for Loo regression by [ne and Vladu 2019). However, lp regression does not have the readily decomposable structure along the coordinates lo regression and novel technique is required in the design of the algorithm. Our approach is also reminiscence of the width- -ndeeeieneeennnnt multiplicative weights update method for solving mixed packing covering linear program, where in each step the algorithm updates the coordinates the maximize the bang-for-buce ratio (Quanrud] 2020). In contrast to MWU, we do not use a mirror map or regularize lp norms make them smooth as in standard approaches. Our scheme allows our method to take much longer steps, where in each step. the coordinates of the dual solution are allowed to change by large polynomial factors and thereby achieve faster running time.
To obtain faster algorithms in the high accuracy regime with logarithmic dependence on the accuracy, we adapt the iterative refinmeen approach of Adil et al. (2019a) and obtain improved running times.
Theorem 1.2. For any p 2 2, there is an iterative algorithm for the Ep regression problem p-2
minzeR": Az=b [|<|pp that solves O (p² log n log (H) subproblems, each of which makes O(1)p=)) calls to solve a linear system of the form ADA & = where D is arbitrary non-negative diagonal matrix, A is a matrix obtained from A by appending a single row, and z a vector obtained from the all-zero vector by appending single non-zero coordinate.
Using the iterative refinement template by (Adil et al. 2019a bJ 2024), we instead use an IRLS solver for the residual problems with improved runtime. The residual solver solves mixed ep lg problem in the form min,: Az=b |a=1ß (0, 2²) only to a constant approximation. Here the challenge lies in the fact that the E2 term makes the dual problem no longer scale-free and thus our low precision solver is not immediately usable. However, by an appropriate initialization of the dual solution and
careful adjustments to the step size, our algorithm achieves the desired (n 5y-3 bound. Since regularized 6p + lg regression problems arise in many applications in machine learning and beyond, our algorithm for the mixed E, + l2 objective is of independent interest.
Finally, we experimentally evaluate our high prcchspreiiioon algorithm. Our algorithm significantly outperforms the p-IRLS algorithm (Adil et al. 2019a) both in the number of linear system solves as well as the overall running time. Our algorithm is significantly faster than CVX solvers and is able to run on large instances, which is not possible for CVX solvers within time constraint.
1.2 RELATED WORK
lp regression problems have received significant attention. Here we summarize the results that are closest to our work. The surveyed algorithms are iterative algorithms where the running time of each iteration is dominated by single linear system solve.
Algorithms based interior point methods use Ô (vr) iterations for any p E [o]
Nesterov
and Nemirovskii 1994p, which was improved to Ö (Vd) iterations for E {1, Lee and
Sidford 2014). Subeck-Cohen-LLe-l (Bubeck et al., 2018) show that this iteration bound is
generally necessary for interior point methods and propose homotopy- based algorithm that uses
Óppoly( (-=+) 111/2-1/pl) iterations for any p & {1, Adil et al. 2019a 2024 introduced
an iterative refinement framework that uses 0(p" n 5 log(2) (2) iterations for any p 2. Using Lewis weight sampling, Janbulapati--iu-SSune Sidford (Jambulapati et al. 2022) improve the method
by Adil et al. 2019a 2024 to (pp - polylo(()) for overconstrained regression problems mingzeR | bll, where A E R and n is much larger than d (the iteration complexity of the prior algorithms will still depend on the larger dimension in this case). Bullins (2088) gives a faster algorithm for minimizing structured convex quartics, which implies an algorithm for é regression with Ō(z)) iterations. Building on the work of CCristiano et al. 2001 1); Chin et al. (2013) for maximum flows and regression, Enne and Vladu (2019) give an algorithm for & and foo regression using (( 2/3 log(1/s) + loga) iterations. This work also uses primal-dual framework but the algorithm and analysis are specific to the special structure of the l and E norm and work only in the low precision regime with p) convergence.
2/3
3

---

Published as a conference paper at ICLR 2026
Algorithm 1 b2 minimizztionnA.....
Input:
Matrix A E Rdx=, vector b IRd accuracy
Output: Vector æ such that Ax = b and ||0|| S (1 + A |
Initialize 20() ming:Az=b |
max : e)²
min {i: (1:(1++ e)¹ |,x0)12
while L < U:
P
3 (1 (1+ + E)
if SubSolver [A, b, M) is infeasible then
-=P+1
else
Let ((+1) be the output SubSolver(A, b, e, M)
== P; 1+1
end if
end while return (1)
Algorithm 2 SubSolver(A,b,,,) M)
Input: Matrix A E _dxa, vector Rd. accuracy e, target value M Output: Vector such that Ar = b and ||0|| S (1 + c)M, approximate infeasibility certificate ||r||, 1.
n-7 L' 0, s(()) =0 while (r(())
1 - 0, (0)
x(1)
arg min Aa o (r(() x²)
a? 1rll
if 2 1irll
7,
(0)
Ma
(1++)M²
for all
1
otherwise
if 70) then return (1)) end if
D Case 1
(0 = (,(1) ((+))) = (t) (1)
if (1) < n पंं (t) ThpT then stt'+1) = s(t') +0(0; 1' end if
if and s(t)/// S (1 + e)M then return s(11/(' end
> Case
t=t+1
end while return (4)
D. Case 3
OUR ALGORITHM WITH poly (t) CONVERGENCE
In this section, we present our algorithm with guarantee provided in Theorem 1.1
Before describing the algorithm, we first introduce some basic notations. For constant a R, abuse the notation and use E R to denote the vector with all entries equal to (the dimension will be clear from context). When it clear from the context, we apply scalar operations vectors with the interpretation that they applied coordinate-wise. For p 2 we let g be such that q and by is the dual norm of the tp norm.
2. OUR ALGORITIM
For ease of notation, it is convenient consider the following equivalent formulation of the problem: For p 2 1, we solve ming:Ar=b || | ming:Ar=5 ||0211, to + c) multiplicative error. We provide our algorithm in Algorithms and We give an overview of our approach and explain the intuition in the following section.
4

---

Published as a conference paper at ICLR 2026
2.2 OVERVIEW OF OUR APPROACH
Our algorithm is based on a primal-dual approach, starting with the following dual formulation of the problem. Using q as the dual norm of p and by duality, write
((r)
min |0ll m:Ar=b
min amAr=b
||
min max (r, (r,x²) brrlllllll,,
max r20: Iill, SI m:rr=b
min (r, x²) = max
20 |r|]
where we defined E(r) = An x²)) The main part of our algorithm is the subroutine shown in Algorithm which takes as input guess M for the optimum value ||0 112 To find an
(1 + <) approximation of the optimum value, the main Algorithm | performs a binary search as we can show that ||==||p is contained in the
follows. Since x(()) is initialized to min ||
range
The algorithm performs binary search over the indices such that (1 + e)²
in that range. Note that the main algorithm only needs perform at most log logg iterations,
each of which makes one call to the subproblem solver.
We now focus on the subproblem when we are given guess M and a target precision c. The goal is to find primal solution that satisfies ||
2p < M (1 + e) dual solution (infeasibility certificate)
which can certify that mitl;:A==b 121 Tr ! )2 This lower bound on the optimal value of the problem tells us that we can increase the guess M.
E(r)
The objective function (r) has very useful monotonicity property: it increases when inereases. The overall strategy of our algorithm is to start with an initial dual solution y((0) (which initialize uniformly to 7. and increase it while maintaining the following invariant
8(rlt+1))-88 - E(r(0) 2 M M - |/000112
or equivalently,
E(r(++1)) - E(r(0) ||+t11)11 ||00 |1g
2 M²
The telescoping property of both sides of ① will guarantee that, if the algorithm outputs a dual ((r)
solution with sufficiently large ||rlly: this solution will satisfy E(r) 2 M
(+te)
||rll i.e
|r >
M I++
To maintain the invariant we have two useful bounds for the change in the objective and
dual solution:
(
E(r(t+2)) E(r(()) 2 Σ 100
Ti
1
(2)
((+1)
911rc0111011 ||((t)))| -|10001 E, (1((tt)) (v!0)"
1
Both inequalities allow to decompose the invariant along the coordinates. That is, maintain the invariant by ensuring for each coordinate that we increase that
(0)
(0)
ri
x
(c)
M
(0)9
(++1)
In order to do this, we update each r multiplicatively, via the term 7,
M
To guarantee fast convergence, we want to increase r!!) as much as possible, by setting a target threshold on 7/ (1), if 7 exceeds the threshold, we update y (t+1) = r((() (7/" (t) 1/9
otherwise, r
5

---

Published as a conference paper at ICLR 2026
remains unchanged. When sure that we have found the corresponding primal solution x with small norm. During the course of the algorithm, we also keep track of iterations with small increases in and use the uniform average over the corresponding primal solutions to obtain an approximately feasible primal solution, in case the algorithm fails to return an infeasibility certificate quickly enough.
longer increase while preserving the invariant, can be
We note that our update approach derived in a completely different way from standard iterative frameworks such as multiplicatives weights updates and, generally, mirror descent. In contrast to these standard approaches, we do not use mirror map or regularize Ep norms to make them smooth. Our update scheme allows our algorithm take much longer steps, and the coordinates of the dual solution are allowed to change by large polynomial factors in each step. This allows us to obtain fast convergence rate.
We outline the necessary lemmas needed to prove Theorem 1.1 before providing complete analysis and proof in Appendix D
Correctness of Algorithm? 2 There two possible outcomes of Algorithm Either it returns
primal solution (Case and Case 2) or a dual certificate (Case 3). the former two cases, Case immediately gives approximate solution. We show in Lemma 2.2 that the returned vector in Case achieves the target approximation guarantee. In Case 3, we use the invariant shown in Lemma 2.1 show that the returned dual solution an infeasibility certificate.
We formalize these statements in the lemmas below.
Lemma 2.1 (Invariant). For all t, we have that if (()) + then 2(r(11)) 10 |
2 M
|
Lemma 2.2 (Case 1). Let r be a dual solution and x arg ming: 22). If
||-1 p S (1+ M then ||>l| S M(1 c).
Lemma 2.3 (Case 3). the algorithm returns (T),, then ((r(T) p
M2
(1+e)2
Convergence of Algorithm We run the algorithm for T iterations. The algorithm terminates if at any point it finds a solution that satisfies the desired bound (otherwise it is unable further increase the dual solution). Otherwise, we show that it must finish very fast. Suppose we run for T = The Tio iterations. Let the iterations in The correspond to those where at least single
coordinate of was scaled by S श (1) क Let Tio be the remaining iterations. The following lemmas give an upperbound on Th and Tio
Lemma 2.4. We have Thi 392
Lemma 2.5. We have Ti, < ((
81/2 Is
क log (=
Since = n 2 (1) 35F we obtain the following convergence guarantee:
at terminates in O ((tt)******
2+29
Lemma 2.6. Algorithm
1o++00 n
(t)
log ( iterations.
Equipped with these lemmas, we give the proof for Theorem
Proof of Theorem 1.1
Returning the problem min Armb ||
p. have the main algorithm
a 10]
executes binary search over the power of (1 + c) in the range
2x (1)
|| so the total
number of calls the subroutine solver O (log log + log !) By Lemma 2.6 the subroutine
solver requires O (((2)
+ デ
log
linear system solves, where q
P
p-2
the dual norm of p/2. Substituting the value of g, obtain the conclusion.

---

Published as a conference paper at ICLR 2026
Algorithm 3 Iteratively Reweighted Least Squares
Input:
Matrix A E _dxx, vector b E IR
Output: Vector such that Ar b that minimizes 1|=||p Initialize z000 arg min
| 12
- 0;
if p S log n I- otherwise
O
16p
P-2
while M 2 [6p(1+e)
|0.0011p
(()) (z)(t) (1);; (p(())
2
Ä < ResidualSolver
P. (o())T
A
[0, Ar(e)
(M - }
if Å is infeasibility certificate or R((1)
} 2M(()) then
M(++) /2, (1++)) (1))
else
M(+1) c x(())) = x0()
64pM
end if
t ++11
end while return (1()
OUR ALGORITHM WITH log (t) CONVERGENCE
3.1 ALGORITIM
In this section, we present our algorithm with guarantee provided in Theorem 1.2 For the ease of the exposition, we consider a slight variation of the problem: for p 2 2, we solve millad ||>||| I| to
(1 + E) multiplicatiie error. We show our algorithm in Algorithms 3 and
3.2 OVERVIEW OF OUR APPROACH
At the highest level, the main algorithm relies on simple yet powerful observation by Adil al.) (2019a), which is that the Ep minimization problem we are attempting to solve supports iterative refinement. Adil et al. (2019a) show that having access to a weak solver which gives a constant factor multiplicative approximation to mixed objective of Ep and & norms suffices to boost the multiplicative error to + e while only making Öp 1/e) calls to the solver. This reduces the entire difficulty of the problem to implementing the weak solver.
More precisely, starting with an initial solution (set to arg ming:AA=s || 12). we maintain M as an upper bound for the function value gap, ie. ||=="0]]/" | S 16pM(()) We show this invariant in Lemma B.2 In each iteration, the algorithm makes call to solver for the residual problem which approximates the function value progress ||==I|²² - residual solution tells us either the progress is too small, in which case we can improve the upperbound on the suboptimality gap by reducing M (t), or the progress is at least Q (M)) in which case we perform the update and obtain new solution. This new solution improves the function value gap
Alp if we update the solution & < I A. The
by at least a factor 1 - S (2(2) and thus the algorithm requires only O p plog |100000 11
calls
to the residual solver. We show this guarantee Lemma
E.2
We give the pseudocode for the residual solver in Algorithm 2024) give algorithms for this solver either via width-reduces multiplicative weights update method
Prior works by Adil et al. 20119 b
'Note that while the residual solver takes as input the original matrix A augmented with an extra row, the least squares problems required by the residual solver reduce to least squares problems involving only A, using the Sherman-Morrison formula. This guarantees that we only require linear system solver for structured matrices of the form A" DA, for non-negative diagonal D).
7

---

Published as a conference paper at ICLR 2026
Algorithm 4 ResidhalSSolver(p,A M)
Input: Matrix A Rd*" vector b R target value M. weight 0 Output: Vector æ such that Ar = b, || I 2M and (0, x²) I min;:Aa=b || + +0,,
approximate infeasibility certificate ||rll 1.
if p
S
log
then
log
arg Ai + 0, x²)
if
r 2p
S 2M then return E else return r end
else
1 - 0, (1)
24-1 1'' = 0, s()) = 0
2
while
( 2
xx arg min Ar= (r(1) + x2)
if =111-11- 2 2M²
(0
M
for all
otherwise
(0)
1/9
a.
if a then return (1) end if
Case
r (1+1) a गंं then st'' +1) s(t') + (0;; (1;; = + end if
(t)
if (!) (1)
if It > 0 and s((')/('
< 2M then return s(t)/ end if
D Case 2
2p
(=1+1
end while
end if
return (*)
D Case
which achieves the state-of-the-an theoretical runtime but does not support a practical implementation via practical IRLS method with suboptimal theoretical guarantee. In contrast, we build ideas from the low precision IRLS solver we have shown in the previous section and design new IRLS algorithm that attains the best of both worlds.
Our residual solver outputs an approximate solution to a constant factor to the objective of the form
min a:Am=b
||ا2 (0.z²)
for 2 and a positive weight vector 0 E IR. We also start with the dual formulation of the problem
min max (1,2²) + (0,=3)= x:A==br: rrllllll
max min (r 0, = max & =00
T ||rl|,
0
r20
where q is the dual norm to and E(r + 0) = ming:Az= + +0,x2) 2). Given a target M, our goal to find a primal solution x that satisfies ||=r|l,, I 2M and (0, x²) < mill::Amm | + (0. x²)
or a dual solution r E I" (infeasibility certificate) which can certify that ming:AA=b || 2p 2 E ( + 0) 2 2r² where & is value set as shown in Algorithm
2x
We distinguish between two regimes: when is sufficiently small, 1 S S log
for which we
log n-I
Ing will show that we can obtain a solution by 0(1) calls to the linear solver. and when log Bogn-T to which we need to pay more attention. In the latter case, similarly to Algorithm[2] we want to maintain invariant
E(r(t+1) + 0) E(r()) + 0) | 11, |r(0
2 M².
Notice the differences between this objective and the problem min || which we solve in the previous section. The &g term (0,x²) makes this objective no longer scale-free. However, this &
8

---

Published as a conference paper at ICLR 2026
(a) x (1-50) 8
(b) size=500 400
(c) size=n x (Tc 100)
(d) size=n x (n2 100)
(e) size=n x (n 50), (f) size=500 x 400
(g) size=n x (n 100) (h) size x (n 100)
Figure 1: Performance random matrices: min ||Ar bll with 10010 10, We compare
algorithm with CVX using SDPT3 and SeDuMi solvers and --IRLS by Adil et al. (2019)). Figures ()).(b)(ee)()) plot the average and standard deviation of number of iterations and time taken by the solvers to find solution over 10 runs. Figures (g). measure over 5 runs.
term does affect the lower bound Σer, (1)
(I)(1)
2
1 Terr
in the change in the objective
Exrl' (2f))²
503 14
2 M in order to guarantee the invariant (0)))+)++))----- (1) 2 M At the same time, if we maintain ||]||, S 1, we can show | I1, -1.11
(eq. (2); thus it suffices to maintain
|
[[(0)
that if the algorithm outputs a primal solution x, the & term (0,=²) S ming: || l|, + (0, x²). This requires us initialize with sufficiently small ||rll,- Algorithm 4 then follows similarly to Algorithm [ with the note that it suffices to obtain only constant approximation. We give the correctness and convergence of Algorithm 4 in Lemma E.1 whose proof is based on the same idea as the analysis for Algorithm 2
The complete analysis of our algorithm is provided in Appendix
4 EXPERIMENTAL EVALUATION
On synthetic data. We follow the experimental setup in Adil et al. (2019b). and build on the provided code
[4] We evaluate the performance of our high-precision Algorithm 3 on the problem min || || blf on two types of instances: (1) Random matrices. The entries of and b are generated uniformly at randomly between 0 and 1; (2) Random graphs. We use the procedure in Adil al.) (2019b) to generate random graphs and the corresponding A and b. The generated graph is weighted graph, where the vertices are generated by choosing a point in [0, 1]10 uniformly at random, each vertex is connected to the 10 nearest neighbors. The edge weights are generated by gaussian type function (by vector. Let B be the edge-vertes, adjacency matrix, W be the diagonal matrix with edge weights. We generate A W1/p B, b B(: + klgg
Calder-Lerman) & (around 10) nodes are labeled in [0, 1] and let g be the label
We vary p and the size of the matrices and graphs, while keeping the error e = 10-10, All imple- mentations were done on MATLAB 2024a on a MacBook Pro M2 with 16GB RAM. We measure the number of iterations and running time for each algorithm and report them in Figures In the appendix, we provide additional experimental results when 1 < p < 2 and when varies.
*The code is available

---

Published as a conference paper at ICLR 2026
Seranty
(a)
8
(b) Number of nodes=500
(c) nodes
(d) nodes
(c) (e)p= (f) Number nodes=500
(g) nodes
(h) nodes
Figure 2: Performance random graph instances: min |r bll with 10 10, We compare algorithm with CVX using SDPT3 and SeDuMi solvers and p-IRLS by Adiil et al. (201bb). Figures (e) measure over runs. Figures (g).( measure over runs.
Table 1: Performance of our algorithm against p-IRLS on six real-world datasets for = 8, e 10 10-10 10
Power Consump-
KEGG
Buzz in
Song Year Pre-
Metabolic
Social Media
Protein Property
CT slices
tion
Naeem
diction
Graf et al.
Hebrail
and
Kawala
Rana
Bertin
(201
and
Asghar (2011)
al.
(2013
Mahieux (2011)
Berard
(2013)
(2006)
48150
57248
1844352
524925
46381 1
Size
41157x9
x385
x27
x11
x77
x90
p-IRLS
48
50
45
50
44
45
iters
Ours
36
42
36
42
36
36
time
p- IRRS
14.3
2.5
32.
28.
1.6
22.5
(s)
Ours
9.2
1.7
15.7
18.1
1.1
13.3
On real-world datasets. We test our algorithm against p-IRLS on six regression datasets from the UCI repository. CVX has excessive runtime and hence is excluded from the comparison. Results are provided in Table
Remark 4.1. Regarding the correctness of the algorithm, we use the output by CVX as the baseline. In all experiments, our algorithm has error within the margin compared with the objective value of the CVX solution (see appendix).
On smaller instances, we compare our algorithm with CVX using SDPT3 and Sedumi solvers and the p-IRLS algorithm by Adil al. (2019)). While CVX solvers generally need fewer iterations find solution, they are significantly slower on all instances than our algorithm and p-IRLS. Our algorithm also significantly outperforms p-IRLS in both the number of iterations (calls a linear system solver) and running time. When the size of the problem and the value of increases, the gap between our algorithm and p-IRLS also increases. On average, our algorithm is 1-2.6 times faster than p-IRLS.
ACKNOWLEDGEMENT
AE was supported in part by Alfred P. Sloan Research Fellowship. AV partially supported by the French Agence Nationale de la Recherche (ANR) under grant ANR-21- ANN-21-CE48-0011 (project COMCOPT).
10

---

Published as a conference paper at ICLR 2026
REPRODUCIBILITY STATEMENT
For the reproducibility purpose, we submitted the source code in the supplementary material. We included the MATLAB implementation by Adil et al. (2019b).
REFERENCES
Deeksha Adil, Rasmus Kyng, Richard Peng, and Sushant Sachdeva. Iterative refinement for & -y-nrrm regression. In Proceedings of the Thirlieth Annual ACM-SIAM Symposium Discrete Algorithms, pages 1405-1424. SIAM, 2019a.
Deeksha Adil, Richard Peng, and Sushant Sachdeva. Fast, provably convergent irls algorithm for p-norm linear regression. Advances Neural Information Processing Systems, 32, 2019b.
Deeksha Adil, Rasmus Kyng, Richard Peng, and Sushant Sachdeva. Fast algorithms for l regressssoon. J. ACM, ((5):34:1-34:44 2024. URL UtL/ttpp:///bl..org//00-
Ahmed EI Alaoui. Asymptotic behavior of \Ceel pl))-based laplacian regularization in semi- supervised learning. In COLT, volume 49 of JMLR Workshop and Conference Proceedings, pages 879-906. JMLR.org. 2016.
T. Bertin-Mahieux. Year Prediction MSD. UCI Machine Learning Repository, 2011. DOI: https://bb..rr/11.244/////3
Sébastien Bubeck. Michael Cohen, Yin Tat Lec, and Yuanzhi Li. An homotopy method for Ep regression provably beyond self-concordance and in input-sparsity time. In Proceedings of the 50th Anmual ACM SIGACT Symposium on Theory of Computing, pages 1130-1137. ACM, 2018.
Brian Bullins. Fast minimization of structured convex quartics. arXiv preprint 1812.11:1 2018.
C Sidney Burrus. Iterative reweighted least squares. OpenStax CNX. Available online: http://enx. 12(2012):6, 2012.
Hui Han Chin, Aleksander Madry, Gary L. Miller, and Richard Peng. Runtime guarantees for regression problems. In Robert D. Kleinberg, editor, Innovations in Theoretical Computer Science, ITCS '13, Benkeley, CA. USA, January 9-12, 2013, pages 269 ACM, 2013. ISBN 978-1-4503- /ttpr///....1115222222
1859-4. doi: 10.1 D0.1445524223362424444 URL 2422469
Paul Christiano, Jonathan A. Kelner, Aleksander Madry, Daniel A. Spielman, and Shang- Teng. Electrical flows, laplacian systems, and faster approximation of maximum flow in undirected graphs. In Lance Fortnow and Salil P. Vadhan, editors, Proceedings of the 43rd ACM Symposium on Theory of Computing, STOC 2011, I, San Jose, CA, USA, 6-8 June 2011, pages 273-282 ACM, 2011. ISBN 978- 1778-145503-069-11 doi: 10. 10.1145/1993636.1990 19936744 URL http:ttps:///mm /Rthttps:/ddi..rrg I0.1145/1993636..
Alina Ene and Adrian Vladu. Improved convergence for E and too regression via iteratively reweighted least squares. In International Conference on Machine Learning, pages 1794-1801, 2019.
Graf, H.-P. Kriegel, M. Schubert, S. Poelsterl, and A. Cavallaro. Relative location of CT slices on axial axis. UCI Machine Learning Repository, 2011. DOI:
Georges Hebrail and Alice Berard. Individual Houschold Electric Power Consumption. UCI Machine Learning Repository, 2006. DOI: http://......r//00223333/C53.
Lingxiao Huang and Nisheeth K. Vishnoi. Coresets for clustering in euclidean spaces: importance sampling nearly optimal. In STOC, pages 1416-1429. ACM, 2020.
Lingxiao Huang, Shaofeng H.-. Jiang, Jianing Lou, and Xuan Wu. Near-optimal coresets for robust clustering. In ICL.R. Reeeeeee,, 2023.
11

---

Published as a conference paper at ICLR 2026
Arun Jambulapati. Yang P Liu, and Aaron Sidford. Improved iteration complexities for overcon- strained p-norm regression. In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 529-542, 2022.
Franois Kawala, Ahlame Douzal, Eric Gaussier, and Eustache Diemert. Buzz in social media UCI Machine Learning Repository, 2013. DOI: htpps//oi.rr///2225//666666
Yin Tat Lee and Aaron Sidford. Path finding methods for linear programming: Solving linear programs in (vrank) iterations and faster algorithms for maximum flow. In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 424 433. IEEE, 2014.
Meng Liu and David Gleich. Strongly local p-norm-cut algorithms for semi-supervised learning and local graph clustering. In NeurlPS, 2020.
Xiangrui Meng and Michael Mahoney. Robust regression on mapreduce. In International Conference on Machine Learning. pages 888-896. PMLR, 2013.
Muhammad Naeem and Sohail Asghar. KEGG Metabolic Relation Network (Directed). UCI Machine Learning Repository, 2011. DOI: https///ab..rg//222335/55555
Yurii Nesterov and Arkadi Nemirovskii. Interior-point polynomial algorithums in convex program- ming. SIAM, 1994.
Kent Quanrud. Nearly linear time approximations for mixed packing and covering problems without data structures randomization. In Symposium on Simplicity in Algorithms, pages 69-80. SIAM, 2020.
Prashant Rana. Physicochemical Properties of Protein Tertiary Structure. UCI Machine Learning Repository, 2013. DOI: http://oo..rgg102233//C5V333
A PROPERTY OF THE ENERGY FUNCTION
We recall the definition energy function and its properties used in the algorithms.
Definition A.1. (Energy function). Given a vector r E R", we let the electrical energy be E(r) = min xr,
Lemma A.1. (Computing the energy minimizer) Given E IRd and r E RR. the least squares problem mine:A==b (r, x²) can be solved by evaluating -1A" (AI(r)-)AAT) where D(r) is the diagonal matrix whose entries given by
The following lemma gives lower bound on the increase in electrical energy when increase Lemma A.2. Given r' 2 and letting I arg min Armm ( xr,, one has that
E(r') (r) N Σ ri
Proof. This inequality follows from the standard lower bound for E(r') - -((r,, which the reader can find Ene and Vladu (0019).
B REDUCING GENERAL REGRESSION PROBLEMS TO THE AFFINE-CONSTRAINED VERSION
In this section we show that the affine constrained version of the problem we consider is in full generality. Formally, we show that any l, regression problem of the form min Az ||N +llp can be reduced the form consider.
12

---

Published as a conference paper at ICLR 2026
Lemma B.1. Let A E R
R. N Rmx", v E Rn and consider the optimization objec-
tive min NN - w1lp: Les
I
be a (1 + +E) approximate solution to the affine-constrained
2
regression problem
min
1
N mmxm
[[
A
xmm
Then is a (1 + approximate solution to the original objective. Furthermore, each least squares subproblem can be solved using two calls to linear system solver for N RN,, and one call to linear system solver for A (NT RN) + AT.
Proof. We augment the dimension of the iterate by introducing m additional variables encoded in vector 2 1². Hence one can equivalently enforce the constraints
Nx - 2 =
Ax
and simply seek to minimize ||2|p instead of ||Ar bl which is the suitable formulation required by our solver. Note that while we do not have any weights on the x iterate, the analysis goes through normally, since in fact tolerates solving a more general weighted l, regression problem.
To solve the corresponding least squares problem, we need to compute
1 (r. (N. u)²) min Ax=b
1 T RNz (N Ru, + " Ru
min Az 2
2
2
max max z N RNI - (N RV, x))
_o Rv + (b - 1-Ax,1) v)
1 T T RNx (N Ru + A v.
max (0.v) + min 2"
Z"
Ro.
where is the diagonal matrix whose entries are given by r. The inner problem is minimized at
x x= (NT RN)* (NT Rv + 1)
which simplifies the problem to max (b, v) (N Rv + AT)) (N (N)* (NT Ro + ATx)
1 = Ro
2
max (D. A (NT (TTRRN" + N Rv.y) 2" A RN)* AT
y
1 JUTRA RN (NTRN)* NT .T.
2" Rv.
2
which is maximized
y= (A(NTRN)* ( (6-A(NTRNY) NT Rw)
so
I = (N) T Rv + (NTRR) RN)* + AT (A RN)* AT)*
(AN)+AA)* (b (N-N(NNT (N RR) + NT Ru)
= RNTR RN AT (A(NTRN)* AT)" Ab-A(NT RN)* N NTRO)) R
We observer that to execute this step we require two calls to a solver for N RN, and one call to solver for A (NT RN) AT.
13

---

Published as a conference paper at ICLR 2026
C SOLVING Ep REGRESSION FOR 1 < p <
In this section show that while our solvers are defined for l regression when p 2 they also provide solutions Eg regression for S < 2. This follows directly from exploiting duality. See Adil al. 2019a), section 7.2 for proof detailed. Here briefly explain why this is the case. Let p, q such that 1 + = 1, S q < 2, and consider the &g regression problem, along with its dual
min ==A==b
||ا|,
max I| vll,
(b,
We can use our solver to provide a high precision solution to the dual maximization problem, which we then show can be used to read off primal nearly optimal solution. Indeed, we can equivalently solve
min (b.2))
|TTTyll wlpp
A
to high precision c no based on which we construct the nearly-feasible primal solution
(b.y) yAT
(ATy)=11
x
To see why this is a good solution, let us assume that we achieve exact gradient optimality for which means that for some scalar X,
A y b-
(5)
First let verify that is feasible. Using
5 we see that:
Ax = A ( (b. y)
(AT (ATy))-1
(b, v)
A A (ATy)
2) b. ||AT ylly"
(b, v)
(ATT vll
= A y/
Additionally can also 5 again obtain that
A L
(v. (.. A (A v) ) = (v, X,
which allows to conelude that
Ax= - b,
so a is feasible. Finally, we can measure the duality gap by calculating
||=1|, (A v)
1 1
[ (A
A
=
(y, 6)
||A ||ATTT]
w11
(y, b)
vllp"
which certifies optimality for b. While in general we do not solve the dual problem exactly, which yields a slight violation the demand for the primal iterate x, this can be fixed by adding to æ a flow H = AT (AAT)* (b Ax) that routes the residual demand. This affects the ly norm only slightly since the residual demand is guaranteed be very small due to the near-optimality of the dual problem. Then can proceed to bounding the duality gap by following the argument sketched above, while also carrying the polynomially small error through the calculation. We refer the reader to Adil al. (2019a) for the detailed error analysis. We have the following theorem.
Theorem C.I. For any p S 2, there an iterative algorithm for the l, regression problem mingeR"; Az=b ||æ||p that solves O (q² log log (2)) subproblems, each of which makes O 3-2 calls to solve a linear system of the form ÄDĀ & = 2, where q = -". D is an arbitrary non- negative diagonal matrix, À is matrix obtained from A by appending a single row. and z is a vector obtained from the all-zero vector by appending single non-zero coordinate.
14

---

Published as a conference paper at ICLR 2026
D PROOF OF THEOREM 1.1
Proof of Lemma 2.7
First we show
3
1
|r((+1)
|| - 11lo11ll-
This is equivalent to show
|
(q 1)
|1-0001 2 1 (0-1)
which can easily be obtained from AM-GM inequality.
Using 3 and Lemma A.2 we have
9 111-00111-1
E (1) (2. (0)))
(1 )
(1))
E(r((-1)) - E(r(1) |((+t) 6 ||/00/1
(((+1)
Σ, 11 (1(41)) -(1(0) (T( (
9 -
Em (1)
[2)-) 2
(0
ri
En (re (1((t1))² ()((1)² 9
((+1) (t) (9)
For such that of > 1. we have ri a
thus
/ |||
1- p(()) (zx (1 1
1
||(00]11111 4
(
1
1
(c- )²=1)²= (1!0)
1
2 ( M
1
(4)
= -22
where the first inequality is due to a(a"-1) conclusion from here.
a(a-1)
2 aa.
1 for a > 1. We can then obtain the desired
Proof of Lemma 2.2 If
x²
18 1-11 9
<
1
for all we have
a7 S (1 + e)² M
| I
which gives
200 < (1 + c) 29 22 r
|r|
We obtain
||0|| 2p < (1 + 2)20 M2p
needed.
15

---

Published as a conference paper at ICLR 2026
Proof of Lemma 2.3
We have that
E(r(T))
E(r(0)) 20 ( -8(((0))
||(2) 1g
||(21111
8(1(0)) + Σi= (1|-(ct2111 - -1|(((111 -M2
(due to the invariant)
||r(7)11,
(| 11, 1) M2
1 (1-
1
= .2.
||r!T)|1
| 11
(since | 2 2
M². cc)
M
(1 c1+c)2"
Proof of Lemma 2.4 Suppose the contrary. Then we claim that the perturbations that scale the dual
solution by 2 S will have increased it lot to the point where ||rl|l > zे Indeed, since is initialized 1 to n//9 in the worst case each perturbation in Thi touches a different coordinate i. Therefore this establishes lower bound of T 8a
||rl1 As this must be at most ca. since otherwise
obtained a good solution per Lemma 2.3
we obtain the conclusion.
Before showing the proof of Lemma 2.5 we claim that we can either look at the history produced in
Tie and obtain an approximately feasible solution, a single coordinate of must have increased lot. Lemma D.1. Consider the set of iterates (f(tt) (11)) used for the iterates in Tio If
1
Σ 2(1)
> M(1 +e)
Tio LET
2p
then there exists a coordinate i for which
Tiot +t1
Σ
a
2
("
>1
Proof. Suppose that
1
Σ x(1)
> >M1+
Tio
rET
2p
Note that by the update rule,
(0) 9-1
(0)g
(1) 9-1
(t)
a
S (1 t)
10,1
M
||r(0]|
(,(00) 9-1
a (t)a (1) 9-1
S (1+ + )
+ lcti1
\
||100
|r(0o119
Hence can write
at0 (0) 9-1
(n(t))
Σ 2
(1+ 9 () Σ
M
S
V
|| 1-
||r(0 ||
TE
(
16

---

Published as a conference paper at ICLR 2026
(1(o))
ai
ai
(0u
(8)
(-)
1-b
< (1 2
Σ LETi.
Σ
||:00||1-1
1
+
r(()) 11
9
r
>1
(by triangle inequality)
(1)
(t) 9-1
9-1
(1+ + ;) Tio +
Σ (ETioraf")>?
a
2
||(0]
9
2p
We obtain
(t)
(1) 9-1
9-1
ai
c
Tio
||((0||
2
TETion 11 21
12p
On the other hand, we have
2p
(t)9
(
(1)
9-1
(1)
9-1
a
r
Σ
9--1
||((t)
11g
rETiona
rETiona!
2p
2p
(T)
Vag
(t)
(T)
max
Σ
Val
(8)
(ETiona
(ETim.a
2p
max
Σ
Vag
(Elio
Therefore there exists such that
2p
Σ
Va
(1)
TT
2p
e
2
which gives us
Tiot 2+1
Σ
Va
>
2
IETio.a (9)
Now we show the proof of Lemma 2.5
Proof of Lemma 2.5 From Lemma D.I we know that there exists coordinate i for which
Tipe zي
Va
(6)
2
tETr. (9)
Furthermore by definition for all iterates in Tio we have that pointwise (1 ++) S (a (0)² < 5
This enables us to lower bound the final value of (r (T)) 9 which is a lower bound on || More
17

---

Published as a conference paper at ICLR 2026
precisely, we have
8
(a. (o.)
(ot"))²
(7)
rETio:
rETug (1 2
Now we can proceed lower bound this coodinate i.e. we lower bound the product in 7 using the lower bound we have in (6).
Intuitively, the worst case behavior i.e. slowest possible increase in (r. (T(T)) 9 two extreme cases:
is achieved in one of the
(i) the og (t) are all minimized i.e.
(a. ) (1 e) in which case + ( log (+) such terms
1 sufficient to make their product 2 ev which means that we are done, since then we have ||| ( ||g
(r:")) e setting The
+ (1 + c) 2a i log (#) i.e Tio θ log (tt is
sufficient to make this happen;
(ii) all the entries maximized, i.e. di
(1) in which case have that their product power q is
at least S 3172
if we set gie _+و S172 In S 2 log [ ie., Tio θ
S1/2 log
7 _
we guarantee that the corresponding r; increases to a value larger than cases capture the slowest possible increase is shown in Lemma F.
es The fact that these two
Therefore we can set
Tio = ((=
1//2 1 g S
log
PROOF OF THEOREM 1.2
First, we give guarantee for the subproblem solver (Algorithm proof follows subsequently)
ifp
log log T
Lemma E.1. For p 2 1, K
Algorithm either returns x such that Ar b.
otherwise
| < 2M and (0, x²)
< ming:Au=t
||~0|, + (0,x²) or certifies that ming:A==b ||
(0,x2) M in O (n calls solve a linear system of the form ADAT ø b, where D is an arbitrary non-negatie diagonal mairix.
The next lemma provides guarantees on the iterate progress in the main algorithm (Algorithm 3
UpS
2logn
Lemma E.2. For p 2 2
logn-I Algorithm maintains that |z 1| ||>01|
--N otherwise
16pM(( and that if (++1) # (9) then
|
([ ||
(1 213pr\ )
(1=c0)1²-- 11=011
Finally, we show the proof of Theorem
Proof. Algorithm terminates when (1) < 16p[1+e |0.0 This gives ||x 1
(1)(0) PP
n. which implies ||= < (1 +c) || 1| and thus [|a-(0) (1 + | 1p. Hence,
Top[1+e]
(1)) is (1 +)) approximate solution. Since Top || 16p(1+c | lp. the number of times
| = O (plog 2). By Lemma E.2 the number of times the
(!) can be reduced O (log |. 1
18

---

Published as a conference paper at ICLR 2026
aa (0)
11
iterate makes progress is
2 pk log
O (p² log r log +) where K Ollog n).
Therefore the total number of calls to the subroutine solver is O (p² log log ") By lemma E. the
subroutine solver makes O (12011) O
(
p-2
calls to a linear system solver. This concludes
the proof.
E. PROOF LEMMA
E.T
We let OPT = minl: Ar
x2
(0, x²) and x² = arg mille::emm ||~011, (0,,=)) We consider
two cases: when p S lemmas:
log logn-1
and when p
log log--I We will prove for each using the following
log
Lemma E.3. For S S log m I- Algorithm either returns x such that Ax = b, ||011 S 2M
and (0, x²) < OPT or certifies that OPT
in O(L) call to solve a linear system.
Lemma E.4. For p Tog n-
log
Algorithm
4 either returns such that Ax = b, | S 2M and
(0,x²) < OPT or certifies that OPT
2
M in O n .द calls to solve a linear system.
To start, we have the following lemma that controls the &g term the objective
Lemma E.5. For such that ||rll, s 1, suppose arg Ar=or + 0, 2²) Then we have (0,x²) S OPT.
Proof. For r with ||rll, S 1, we have
(0, x, & (r + +0, ( (by definition of x)
< <1 "1l + 10,00, (x < OPT.
Now, let turn to the first case when p S
log n log oggn-1
We give the proof for Lemma E.3
Proof of Lemma E.3
When 1 S p 5
log
we have g = -
log n. Algorithm
computes
log n X-I
f
min (r + 0, m:Az=b
where ri for all i.
Since |rlll = 1, if | & 2M, by Lemma E.5 we immediately have 1711 S 2M and (0,x²) s OPT.
Assume that | 12 > 2M. We have
|x)) + (0,()))) (,(()))))+ (0,
OPT
(0 + (x²)²) > (0 (+))())))
1 - | 11 - 2
(since | 2 13
n
2 |
(since g 2 log m)
M
21/² >
For the case when Algorithm 2
logn log n- the proof for Lemma E.4 follows similarly to the analysis of We proceed showing the following invariant.
19

---

Published as a conference paper at ICLR 2026
Lemma E.6 (Invariant), For all t, we have that if + 1 then E( 40)-
+0)
2 M
|
11,
|
Proof. Using Lemma A.2 we have
9. 1|100111111
Σ3 Pi + 8.
(1)
2
+0000
E(r((t1) + 0) - -E(r(() + +0)
+0
[ |r(ct2/|| 1 -11-0-1-1---- |r(t)
Σ,
(1( ((((t))
(0)
9. 1111111111
Σ (110)) +0,
((41) (11))
Σ,
((tt1)
-
Σ
(e+1) (1))
(0
((+1)
111111111
(1)
2 (0) ((+1) - r!(1)
1 (
Σ ΣingD>
(1!1))
where in the second inequality we use
1 +0.
r لل+ा
for (t+1)
ri
2 ri 0 0. For i such that
(t)
ri
(0) > 1, we have r: ((-1)
a; (,,),)
thus
g. ||:0001111 l 2
(1)
(())
ri - -r!
[[(4t))
1
((()))
(0th (1) M
(1)
9
7:
1
2 71
(1)
M
(1)
M
where the first inequality is due gla-1) conclusion from here.
2 1 for a We can then obtain the desired
ala--1) as
Lemma E.7 (Case 1). Let r be a dual solution and arg mig + 0, If
2
< 2M then [1llla s 2M and (0,x2) OPT.
x²)
--
Proof. If
x²
|m||
S
2MM
---1
for all i we have
r--
x7 < 2
|||| 9
which gives
2 < 22p 22 r?
||r11
We obtaín
|| | 2r < 22r MP
as needed. The second claim comes directly from Lemma
E..
20

---

Published as a conference paper at ICLR 2026
Lemma E.8 (Case 3), If the algorithm returns r(P) then &
[ + )0 2 2q
p(t)
Proof. We have that
E(r(T) +0)
E(r(0)
0 + Ei=o (E(r((+1)
+ +0) - -E(r() + 0))
|r(2))
||r(()
Σi=
(|r(ttD)|,
||~(60 [1,)
-22
(due to the invariant)
|||(T]]
1, - ||
-12
|
11
24-1
M
24
(since
|-000 1
2q 1
|
29
M
11,
(since
2 1).
29
Finally since ||r1011,
21
(())
E(r(T) + +0) |,|3)|1,
M²
E
|r(T]]
24
Convergence Analysis We run the algorithm for T iterations. The algorithm terminates if at any point finds a solution that satisfies the desired bound (otherwise it is unable to further perturb the dual solution). Otherwise, show that must finish very fast.
Suppose we run it for - Thi + Tio iterations. Let the iterations in The correspond to those where at least single r; was scaled by 2 S = उपका Let Tio be the remaining iterations.
Lemma E.9. We have Th I
2n
34
Proof. Suppose the contrary. Then we claim that these perturbations alone will have increased r lot to the point where ||1l 2 1. Indeed, let r; be the current value of coordinate and ry be
its value after being increased, and assume that r r; 2 S. Since is initialized 2q-1 2q n in the
worst case each perturbation in The touches a different Therefore this establishes lower bound of Ths [2g-1 2 Thi
2n on ||- As this must be at most 1, since otherwise we obtained E.8 we obtain the conclusion.
2q
good solution per Lemma
Now we claim that we can either look at the history produced in Ti and obtain an approximately feasible solution, single coordinate r; must have increased lot.
Lemma F.10. Consider the set of iterates (x(0),x((1))
x((1) used for the iterates in Tio- If
1
Σ
(0
> 2M
Tio
LET
2p
then there exists coordinate i for which
Σ
Tio
a
4
(0>1
21

---

Published as a conference paper at ICLR 2026
Proof. Suppose
1
Σ
x(t)
> 2M. Note that by the update rule,
The
2p
(1(1)) 9-1
(2/0 (T(())
9-1
(t)
a
Ij
S V2
|||00 9-1 + 1nts
M
||r00111/-
Hence we can write
(1!0)) 8-8
(1(t))
agt
Σ
22 Σ
S
M
||(t] 1
|
BETi
2p
rEToma (1
(rit)) 9-1 8-1
109
(0(0))
S 2/ Σ
+
(
||(((]
1
9--
te
rETioin!
2p
(by triangle inequality)
(0)
(x()) 9-1
a
= /22Tio +
Σ
||1(1) |1g-1
We obtain
(T(00) 9-1
a
Tio
Σ
2 2 v2)
Tio 2
||((0]
110
2
LETiora
(0
>1
1p
On the other hand, we have
2p
2p
(t)g
(1)
9-1
at (0) (1 (1(tt1) 9-1
ai
a
r;
Σ
Σ
||:0011011
|||0 11:
(ETiona!"
>1
(ETio:a >1
(T) 9
2p
2p
|r( 4
(t)
Σ
Vala
|
a
(0)
9
max
("
g
(ETiv.n
(1
2p
<
2q
max
Σ
Va
(since |,00
29
2q
2g
(ETiona >1
2p
< max
Σ
(t) Valt
(since q 2 1)
{" >1
Therefore there exists i such that
Σ
(t)
Tio
2p
Va
> 2
2
(ETio.a(")>1
>1
22

---

Published as a conference paper at ICLR 2026
which gives us
Tio 1
2 Too
Σ
1
since 2 1.
2
2
4
( >1
This lemma enables us to upper bound Tio-
Lemma E.11. We have Tio S @
1112 In
Inn + In
Proof. From Lemma E.10
know that there exists coordinate i for which
Ti
Σ
a
(8)
4
rETroin
(412)
Furthermore by definition for iterates in Tio we have that pointwise of
S and
af - -3,
1/9 2 24. This enables us to lower bound the final value of (r!")) which is lower
bound on |r(r))1 More precisely, we have
(C+1)
a
(t)
thus
(
II
(a.²)²
29 29-11
1
II
a (a)))²
2g
tETro:a (1 >1
tETio:a!">1
Now we can proceed to lower bound this r, ie. we lower bound the product in using the lower bound we have in
r (T)
Similarly to the previous section, the worst case behavior i.c. slowest possible increase achieved in one of the two extreme cases:
is
(i) the a; are minimized i.c. af 2 in which case (Inn) such terms are sufficient to make their product 2 2n 2 2qm 24-1 which means that we are done, since then we have ||
1; so setting Ti 2 θ (In ) is sufficient to make this happen;
(ii) the entries are maximized, i.e. a
which case we have that their product to power g
at least S 4g//2 > 2n 2
Ti.g In S 2 In 2n, ie, Tio 2 2--I. so if we set 151/z
2qm
851/2 In
luin), we guarantee
that the corresponding r; increases to value larger than 2. The fact that these two cases capture the slowest possible increase shown in Lemma
Therefore we can set
Tio = 0 ( 1//
Inn + lnn)
InS
Finally, by the choice S = I we obtain the runtime guarantee.
Lemma E.12. Algorithm terminates in O T
iterations.
Proof of Lemma
E.4
The proof of Lemma
immediately follows from Lemmas E.7
E.8
and
E.2
PROOF LEMMA E.2
Proof of Lemma
E.2
We define the function res,, follows
res (A) = (g. A RR, A2) A1
23

---

Published as a conference paper at ICLR 2026
where g |z| p-2 R 2 |xr We use the following property of this function from Adil et al.
(2019a 2024 For = 16p and for all A
|0
1
2 res (A);
(10)
P
A
|
pp
S Ares (A). 2
(11)
p
We prove the claim by induction.
|| 12
|0x00 1,"011--10=01
For = 0, we have (0)) = 16p
16p
Now assume that we have ||~|00
||
|1x
16 We have two cases.
Case ResidualSolver returns an infeasibility certificate or ResidualSolven returns a primal solution Å such that (889,7) (M In both scenarios, using Lemma E.1 we have
min AA=0
|
+ (M 27"
2 ())
(9 ,A) M(5
Hence for all A such that AA 0, (9 A) MM
2 either || 2 (M + || 2 M())
or (M 2-p (R(*) A?) 2 (MI( write (g A) M) for some constant E R. We obtain either | 2 a M()
(R(), A 2 ( For all A such that AA 0,
(R(9, A2) 2 a² (()) and thus for all A
1 5a res,()) (A) S (1)
M 1
min {a [a2,aP]]
2 = M(+11)
2011) - for X= = 16p. Using property A/p
We write A
11 of the res2, we have
((+2))
|
(1)
| 11
(0((0) I|"
(1)
x A
p p
S Ares, (A)
(6p
Case 2. We have (R, (R,
< (M and ||||
- 4VK(M(()) and (9.4) M00
|o100/1 p
(111) 1|" = |,000 1|3
(1)
A
64pM
Â
res
64k
2
R,
64x
64k
64k
(1)
(A)
(4)
2
211
24P
M(
M()
(1)
(since p 2 2, 2 1)
27 K
2 K
25
24

---

Published as a conference paper at ICLR 2026
M
2 2m (10t011 - 11 1|0-11²))
2
2
from which we obtain
((+1)
||00 "Ip
|| |1 pp
| |xx00|/ -110= 11
(1 (1=20)
||00113
||
21spk)
as needed.
F LOWER BOUND LEMMA
Lemma F.1. Let a set of nonnegative reals BI,... 3r such that 1+e S 3, S S, and Σ 8 2 K,
where 2. Then for any k one has that
k
II min Sath (1 + e) a)a i=1
Proof. Consider a fixed and let us attempt minimize the product of B subject to the constraints. want to minimize Σ log(3;)) which is concave function. Therefore its minimizer is attained on the boundary of the feasible domain. This means that for some 0 S 1/' s k - 1, there are k' elements equal 1 1 k 1-1 1- k' equal to S, and one which is exactly equal to the remaining budget, i.e. (K- k'(1 + e) (k: - k))'/r), which yields the product (1 + e)*² s -1 (K - k( + e)2// - (k - 1 - k)sttr) This can be relaxed by allowing & and k' to be non-integral. Hence aim minimize the product (1 + e)* gk-k'-1 S//r K.
W.l.o.g. have Σ B K. Equivalently
subject to k1(1 + / (k
Finally, we observe that we can always obtain a better solution by placing the available mass on single one the factors, i.c we lower bound either by S S or (1 + e) e) whichever is lowest.
&
G ITERATIVE REFINEMENT
In this section provide general technique for solving optimization problems to high-precision, by reducing then to an adaptive sequence of casier optimization problems. which only require approximate solutions. This formalizes the minimal requirements for the iterative refinement scheme employed in Adil et al. (2019a b) to go through. We state the main lemma below.
Lemma G.1. Let D C R be a convex set, and let f f D be a convex function. Let nj 2 be scalar, and suppose that for any E 'D there exists function h that approximates the Bregman divergence at the sense that
1 he (gu) Sf f()) ((t)) (Vf (z), h2 (6)
r)
Given access to oracle that for any direction can provide x x-approximate minimizers (v, 6) by (0) in the sense that it returns 6 such that + E 'D and
min (v,5) + hz (8) (6)) U+6ED
(v, + hr (5t) S
along with an initial point xo E D, O x In f(x)))((=))
calls the oracle one can obtain
"
point such that f (x) S f (x*) + e, where x² arg min (x:).
Proof. Let 5 be the x-approximate minimizer of (VS 54) h (6#), which by definition satisfies:
1
(V/ (z) .t) + h1p (5t) <
min (V) ( 5) + ho (0)
(12)
k
V+EED
25

---

Published as a conference paper at ICLR 2026
Updating iterate to x' + 5² we can bound the new function value as
(x + 0
= =x + V (z), )*) + he (6*)
(Bregman divergence upper bound)
(x) + "
(Vi (3) xx - 2) + 1 hr (n (x*
(n(x" (x(x=))) -x))
(using 12
= ( + (VV) (z), 1- - -z)+ (f(x2) (x*)- - (x) - (Bregman divergence lower bound)
= f 2 (())) k
from where we equivalently obtain that
f (2+ f (x*) k
(10))
Therefore to reduce the initial error f (re) f (x") to it suffices to iterate O K In ((a)) (xx") times.
The following lemma provides sandwiching inequality for the Bregman divergence of ||==1p--
Lemma G.2 (Adil al. (019b), Lemma B.1). For any x, 5 and p 2 2, we have for r p-2 and
g pr
(1,8²) 2 |1 1 + 816 ||>1|/ - (g. 8) S 2p² (r, P" 11816
As corollary we see that the function hr (6) 2p² (1p-2 82) + p" |811 satisfies the inequality required by Lemma G.1 for 7 4p approximately minimizes mixed &g lp regression objectives, one can efficiently generate a high precision solution.
1. We can thus conclude that given access to an oracle that
Corollary G.1. Consider the Ep regression problem miny:B f=d |SI Given access to an oracle that can compute x-approximate minimizers the optimization problem
:= min A5=
(pp211 AA) + 2p² ( A + r" |S
in the sense that it returns A satisfying B Af = and
| < Lv=
(pf If)) + 2p² (f 2 A + p2
along with an initial point fo. satisfying T = d. in Ep - 1
calls to the oracle one
can obtain a point f such that [l. < Il + where f E arg min f=nd Il
Proof. Using Lemma G.2 verify that the function hf (Af) 2p² (fr Af2 + p" | satisfies
- hy (nAS) < [I + + 1S (pf (Af
for 7) 4p Therefore by Lemma G.I can need O (spln [LFoll. -6ol 1If" E-additive error to the regression problem.
iterations to obtain an
2
H ADDITIONAL EXPERIMENTAL RESULTS
Correctness of solution. In Figure 3 plot the error of the solutions outputted by our algorithm
and p-IRLS against CVX in the random matrices and random graphs instances for e = 10 10 In all cases, the error is below
26

---

Published as a conference paper at ICLR 2026
(a) matrix size=n X (n (b) malrix size=500 x 400 Graph of n nodes, p 8(d) Graph of n 500 nodes
50). p 8
Figure 3: Error of the solution against CVX/SDPT3 solution in log 10 scale.
13
(a) matrix size=500 400 (b) matrix size=500 x 400 (c) Graph of n 500 nodes(d) Graph of n 500 nodes
(e) matrix size=2500x 24000() matrix size=2500 2400(g) Graph of 10000(h) Graph of 10000 nodes
nodes
Figure 4: Performance when varying random matrices and random graphs instances.
When varying e. In Figure we plot iteration complexity and runtime in seconds of our algorithm, p-IRLS and CVX when varying e. Note that, CVX does allow varying this parameter. In all experiment, we fix p 8. For large instances, we only consider our solution against p-IRLS.
For 1 < < 2. In Figure we plot iteration complexity and runtime in seconds our algorithm, p-IRLS and CVX on random matrices of size n x (n 100). We fix 10-1h We test with 1.1 and = 1.9.
25
16 10c 100
(a) 1.1
(b) 1.1
(c) (c)p= 1.9
(d) 1.9
Figure 5: Performance when = 1.1 and p = 1.9 on random matrices of size X (n - 100).
27


Reviews:
[{"review_id": "kCqQFWmyYQ", "reviewer": "Reviewer_rkw2", "summary": "The paper studies the classic ℓp regression problem, which generalizes least squares (p = 2) and appears in robust regression, clustering, and semi-supervised learning. The authors propose a new Iteratively Reweighted Least Squares (IRLS)–based algorithm for ℓp regression that: Achieves the same theoretical iteration bound as the complex SODA 2019 / JACM 2024 algorithm of Adil–Kyng–Peng–Sachdeva,\nBut retains the simplicity and practical efficiency of the p-IRLS (NeurIPS 2019) method.", "questions": "The authors should include a comparison with existing IRLS algorithms and their variants, along with an explanation of why the proposed approach outperforms these methods.\n\nThe paper presents two versions of the algorithm—one for the low-accuracy regime and another for the high-accuracy regime. A discussion on their practical usage would be valuable. In particular, it would be helpful to clarify under what conditions each version is preferable, and whether a hybrid strategy (e.g., using the low-accuracy regime for initialization followed by the high-accuracy regime for refinement) would offer additional benefits.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "The IRLS method is widely used in practice but typically lacks strong theoretical guarantees. This paper proposes a variant that offers both rigorous theoretical foundations and favorable computational efficiency. The authors reformulate IRLS within a primal–dual invariant framework, deriving the update rule from a monotonic invariant of the dual objective. The algorithm adaptively rescales weights to preserve this invariant, achieving provably fast convergence for all values of p. The proposed approach is novel and potentially impactful beyond this specific setting, given the broad applicability of IRLS across various domains.", "weaknesses": "There should be more discussion on the intuition behind why this reformulation leads to improved convergence, as the underlying mechanism is not immediately clear—especially in contrast to the standard IRLS method, which is more intuitive.\n\nThe case of p<1 is particularly important due to its relevance in applications such as robust regression. However, this scenario is not addressed in the main paper and is only briefly mentioned in Appendix C, where the conclusions remain unclear. Given its significance, this case should be discussed in the main text, with the key findings and implications explicitly stated.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "extra_scores": {}}, {"review_id": "HvJd371ySI", "reviewer": "Reviewer_3oUi", "summary": "This paper presents new algorithms for solving ℓₚ regression problems using the Iteratively Reweighted Least Squares (IRLS) framework. The proposed methods achieve state-of-the-art iteration complexity while maintaining practical simplicity.", "questions": "The algorithm’s performance under modern hardware acceleration is not discussed. I'm more interested in the parallel algorithm.", "limitations": "", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper establishes a solid theoretical foundation by matching the best-known asymptotic complexity of Adil et al. (JACM 2024) with a simpler iterative structure. Its primal-dual, invariant-based design offers a clean and principled acceleration of IRLS, effectively bridging theory and practice through provable guarantees and strong empirical performance. Experiments on synthetic and real datasets confirm notable gains in speed and accuracy, supported by clear and well-organized mathematical exposition.", "weaknesses": "None", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "extra_scores": {}}, {"review_id": "19N8eg6WgW", "reviewer": "Reviewer_eXNj", "summary": "In this paper, the authors present an algorithm for solving $\\ell_p$ regression problems using the IRLS method. The algorithm is based on a clever update rule motivated by the dual formulation of an equivalent problem, and it achieves better complexity than other versions. In particular, the number of calls to the linear system solver matches the state-of-the-art theoretical order.", "questions": "I don't understand why to include the number of iterations of cvx sedumi/SDPT3 in Fig. 2. The nature of the iterations is quite different. I think it would be better to include them only in the time comparison.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is very well written in general. Although my expertise lies far from this area, I was still able to understand the main ideas and appreciate the contributions. The experimental results are very good.", "weaknesses": "This is not a weakness of the paper itself, but rather of this type of work in general. I believe that one of the strongest contributions of the paper lies in Theorems 1.1 and 1.2, but the proofs are relegated to the appendix (and I did not check them). I have always found that such papers are more appropriate for journal publication, where these key aspects can be better appreciated.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "extra_scores": {}}, {"review_id": "mY5LKvXrLX", "reviewer": "Reviewer_BjD5", "summary": "The paper gives an algorithm for the p-norm regression problem for $p> 2$. The known algorithms include the algorithm by AKPS’ JACM 2024 which gives the best theoretical guarantees but a complicated algorithm, and a practical algorithm by APS’19 which outperforms all existing implementations at the cost of a worse iteration complexity. This paper bridges the gap by giving an algorithm that matches the guarantees of AKPS as well as outperforms the algorithm by APS’19. \n\nTheir algorithms is an IRLS style algorithm similar to the approach of APS’19 but their algorithm and analysis are quite different. The analysis follows a more primal dual approach similar to the works of Ene-Vladu ICML’2019. They first use an approach similar to EV’19 to give a low accuracy solver and then use the iterative refinement framework from AKPS to boost the accuracy. In order to do so they require modifying their algorithms to work for a regularized version of the problem.", "questions": "1. The algorithm looks more complicated than APS’19 and also of a similar level of complexity as AKPS JACM’24. It looks like a different approach, and maybe the authors could say something about the fundamental difference with previous works. I am not sold on the simplicity of the algorithm, regardless I think it is interesting as a standalone algorithm.\n2. Also the difference in iteration counts in the experiments seems to be some fixed constant factor of maybe 2 or something. I am wondering if this is because APS’19 is not implemented well enough or this algorithm is genuinely faster due to the improved rates of convergence?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Interesting theoretical result which can be of independent interest.\n2. Bridges the gap between theory and practice for this problem.\n3. They also test their algorithm on some real-world data sets which I dont think has been done before for this problem.", "weaknesses": "See questions. No major weaknesses!", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)