Improved Convergence for $\ell_\infty$ and $\ell_1$ Regression via Iteratively Reweighted Least Squares

TL;DR

This paper introduces a simple IRLS algorithm for $\, ext{l}_ ext{infinity}$ and $ ext{l}_1$ regression that converges efficiently with a running time independent of input conditioning, outperforming previous methods.

Contribution

A new natural IRLS algorithm with provable convergence for $\, ext{l}_ ext{infinity}$ and $ ext{l}_1$ regression, featuring a runtime independent of input conditioning and sublinear in $1/\, ext{ extepsilon}$.

Findings

Converges to a $(1+ ext{ extepsilon})$-approximate solution in specified iterations.

Runtime is independent of input matrix conditioning.

Outperforms previous algorithms by at least a factor of $1/ ext{ extepsilon}^2$.

Abstract

The iteratively reweighted least squares method (IRLS) is a popular technique used in practice for solving regression problems. Various versions of this method have been proposed, but their theoretical analyses failed to capture the good practical performance. In this paper we propose a simple and natural version of IRLS for solving $\ell_\infty$ and $\ell_1$ regression, which provably converges to a $(1+\epsilon)$-approximate solution in $O(m^{1/3}\log(1/\epsilon)/\epsilon^{2/3} + \log m/\epsilon^2)$ iterations, where $m$ is the number of rows of the input matrix. Interestingly, this running time is independent of the conditioning of the input, and the dominant term of the running time depends sublinearly in $\epsilon^{-1}$, which is atypical for the optimization of non-smooth functions. This improves upon the more complex algorithms of Chin et al. (ITCS '12), and Christiano et al. (STOC '11) by a factor of at least $1/\epsilon^2$, and yields a truly efficient natural algorithm for the slime mold dynamics (Straszak-Vishnoi, SODA '16, ITCS '16, ITCS '17).

Equations

$$O\left(\frac{m^{1/3}\log (1/\epsilon )}{{\epsilon }^{2/3}}+\frac{\log m}{{\epsilon }^2}\right)$$

$$x^{(t+1)}=\arg \underset{Ax=b}{\min }⟨\frac{1}{|x^{(t)}|},x^2⟩,$$

$$x^{(t+1)}=\arg \underset{Ax=b}{\min }⟨\frac{1}{\sqrt{(x^{(t)}{)}^2+{\eta }^2}},x^2⟩.$$

$$c_i^{(t+1)}$$

$$x^{(t+1)}$$

$$\underset{x:Ax=b}{\min }\Vert x^2{\Vert }_{\infty }=\underset{x:Ax=b}{\min }\underset{r\in {\Delta }_m}{\max }⟨r,x^2⟩=\underset{r\in {\Delta }_m}{\max }\left(\underset{x:Ax=b}{\min }⟨r,x^2⟩\right):=\underset{r\in {\Delta }_m}{\max }{\mathcal{E}}_r(b),$$

$$\frac{{\mathcal{E}}_{r{}^{\prime }}(b)-{\mathcal{E}}_r(b)}{\Vert r^{\prime }-r{\Vert }_1}\ge M^2,$$

$$\underset{x:Ax=b}{\min }\Vert x{\Vert }_1^2=\underset{x:Ax=b}{\min }\left(\underset{c\in {\Delta }_m}{\min }⟨\frac{1}{c},x^2⟩\right)=\underset{c\in {\Delta }_m}{\min }\left(\underset{x:Ax=b}{\min }⟨\frac{1}{c},x^2⟩\right)=\underset{c\in {\Delta }_m}{\min }{\mathcal{E}}_{1/c}(b).$$

$$\frac{\frac{1}{{\mathcal{E}}_{c{}^{\prime }}(b)}-\frac{1}{{\mathcal{E}}_c(b)}}{\Vert c{}^{\prime }-c{\Vert }_1}\ge \frac{1}{M^2}.$$

$$\underset{Ax=b}{\min }\Vert x{\Vert }_{\infty },$$

$$b^{\top }(AD(r{)}^{-1}{A}^{\top }{)}^{+}b=\underset{x:Ax=b}{\min }⟨r,x^2⟩\le ⟨r,(x^{\ast }{)}^2⟩\le \Vert r{\Vert }_1\Vert x^{\ast }{\Vert }_{\infty }^2=\Vert x^{\ast }{\Vert }_{\infty }^2,$$

$$\underset{Ax=b}{\min }\Vert x{\Vert }_1,$$

$$b^{\top }(A\mathbf{D}(c)A^{\top }{)}^{+}b$$

$$=\sum _{i=1}^m\frac{1}{c_i}\cdot {\left(\mathbf{D}(c)A^{\top }(A\mathbf{D}(c)A^{\top }{)}^{+}b\right)}^2$$

$$=\sum _{i=1}^m\frac{1}{c_i}\cdot x_i^2.$$

$$\Vert x{\Vert }_1=\sum _{i=1}^m\frac{|x_i|}{\sqrt{c_i}}\cdot \sqrt{c_i}\le \sqrt{\left(\sum _{i=1}^m\frac{x_i^2}{c_i}\right)\left(\sum _{i=1}^m{c}_i\right)}=\sqrt{b^{\top }(A\mathbf{D}(c)A^{\top }{)}^{+}b}.$$

$${\Vert A^{\top }{\Phi }^{(t^{\prime })}\Vert }_{\infty }={\Vert A^{\top }\cdot \sum _{t\in S}\frac{{\varphi }^{(t)}}{b^{\top }{\varphi }^{(t)}}\Vert }_{\infty }$$

$$={\Vert \sum _{t\in S}\frac{A^{\top }{\varphi }^{(t)}}{b^{\top }{\varphi }^{(t)}}\Vert }_{\infty }\le {\Vert \sum _{t\in S}\left|\frac{A^{\top }{\varphi }^{(t)}}{b^{\top }{\varphi }^{(t)}}\right|\Vert }_{\infty }\le \frac{t^{\prime }}{(1-\epsilon )M}.$$

$$A_{ve}=\{\begin{array}{ll}+1& \text{if }e\in E^{+}(v),\\ -1& \text{if }e\in E^{-}(v),\\ 0& \text{otherwise.}\end{array}$$

$${\mathcal{E}}_r(x)=⟨r,x^2⟩.$$

$${\mathcal{E}}_r(b)=\underset{x:Ax=b}{\min }{\mathcal{E}}_r(x),$$

$${\mathcal{E}}_r(b)$$

$$=\underset{\varphi }{\max }2\cdot b^{\top }\varphi -\sum _{i=1}^m\frac{{\left(A^{\top }\varphi \right)}_i^2}{r_i}$$

$$={\left(\underset{\varphi :b^{\top }\varphi =1}{\min }\sum _{i=1}^m\frac{{\left(A^{\top }\varphi \right)}_i^2}{r_i}\right)}^{-1}.$$

$$x_i=(A^{\top }\varphi {)}_i/r_i.$$

$${\mathcal{E}}_{r{}^{\prime }}(b)\ge {\mathcal{E}}_r(b)+\sum _{i=1}^m{r}_i{x}_i^2\left(1-\frac{r_i}{r_i^{\prime }}\right).$$

$${\mathcal{E}}_{r{}^{\prime }}(b)$$

$$=2\cdot b^{\top }\varphi -\sum _{i=1}^m\frac{{\left(A^{\top }\varphi \right)}_i^2}{r_i}+\sum _{i=1}^m\frac{{\left(A^{\top }\varphi \right)}_i^2}{r_i}\cdot \left(1-\frac{r_i}{r_i^{\prime }}\right)$$

$$={\mathcal{E}}_r(b)+\sum _{i=1}^m\frac{{\left(A^{\top }\varphi \right)}_i^2}{r_i}\cdot \left(1-\frac{r_i}{r_i^{\prime }}\right).$$

$$\frac{1}{{\mathcal{E}}_{1/c{}^{\prime }}(b)}\ge \frac{1}{{\mathcal{E}}_{1/c}(b)}+\frac{1}{{\mathcal{E}}_{1/c}(b{)}^2}\cdot \sum _{i=1}^m{c}_i(A^{\top }\varphi {)}_i^2\left(1-\frac{c_i}{c_i^{\prime }}\right).$$

$${\mathcal{E}}_{1/c{}^{\prime }}(b)\le \sum _{i=1}^m\frac{1}{c_i^{\prime }}{x}_i^2=\sum _{i=1}^m\frac{1}{c_i}{x}_i^2+\sum _{i=1}^m\frac{1}{c_i}{x}_i^2\cdot \left(\frac{c_i}{c_i^{\prime }}-1\right)={\mathcal{E}}_{1/c}(b)+\sum _{i=1}^m\frac{1}{c_i}{x}_i^2\left(\frac{c_i}{c_i^{\prime }}-1\right).$$

$$\frac{1}{{\mathcal{E}}_{1/c{}^{\prime }}(b)}$$

$$\ge \frac{1}{{\mathcal{E}}_{1/c}(b)}+\frac{1}{{\mathcal{E}}_{1/c}(b{)}^2}\cdot \sum _{i=1}^m{c}_i(A^{\top }\varphi {)}_i^2\left(1-\frac{c_i}{c_i^{\prime }}\right),$$

$${\mathcal{E}}_{r/\Vert r{\Vert }_1}(d)\ge (1-\epsilon {)}^2{M}^2.$$

$$\frac{{\mathcal{E}}_{r^{(t+1)}}(d)-{\mathcal{E}}_{r^{(t)}}(d)}{\Vert r^{(t+1)}-r^{(t)}{\Vert }_1}\ge M^2.$$

Full text

In this paper we propose a simple and natural version of IRLS for solving $\ell_{\infty}$ and $\ell_{1}$ regression, which provably converges to a $(1+\varepsilon)$-approximate solution in $O(m^{1/3}\log(1/\varepsilon)/\varepsilon^{2/3}+\log m/\varepsilon^{2})$ iterations, where $m$ is the number of rows of the input matrix. Interestingly, this running time is independent of the conditioning of the input, and the dominant term of the running time depends sublinearly in $\varepsilon^{-1}$, which is atypical for the optimization of non-smooth functions.

This improves upon the more complex algorithms of Chin et al. (ITCS ’12), and Christiano et al. (STOC ’11) by a factor of at least $1/\varepsilon^{2}$, and yields a truly efficient natural algorithm for the slime mold dynamics (Straszak-Vishnoi, SODA ’16, ITCS ’16, ITCS ’17).

1 Introduction

Regression problems are fundamental primitives in scientific computing. Among these, $\ell_{\infty}$- and $\ell_{1}$-regression are their hardest instantiations, since through standard reductions they can be shown to be equivalent to linear programming.

While the series of works on these topics is truly extensive and diverse, the simpler methods have pervaded into the realm of practical applications. Among these, an extremely popular scheme known for its simplicity is the iteratively re-weighted least squares (IRLS) method. The idea behind it is to reduce optimization problems to iteratively solving a series of weighted $\ell_{2}$-minimization problems, where the weights are adaptively chosen in such a way that the resulting solutions from the sequence of least-squares problems converge to the sought optimal point. In particular, due to its relevance in signal processing, $\ell_{1}$ regression is a very important application of IRLS [C*+*06, CY08].

Despite the fact that various versions of this method have been studied ever since the 60’s [Law61, Osb85] theoretical understanding of their convergence has lacked. Recent works have attempted to fill this gap, and offer provable guarantees [DDFG10, SV16c, SV16b, SV16a], some of them inspired from the interpretation of this method as a dynamical system. In particular, we note the Physarum dynamics, which have been studied in a completely different context [IJNT11, JZ12, TKN07, BMV12, BBD*+*13] in order to justify an experiment which revealed that a unicellular organism, the slime mold, could solve the shortest path problem in a maze [NYT00]. The fact that these dynamics are essentially just a version of the IRLS method was observed in [SV16c].

Returning to the more classical world of algorithm design and analysis, it is worth observing that existing analyses of IRLS methods fall into one of the following two categories: (i) they show convergence only when the problem is properly initialized, or (ii) the guaranteed running time is prohibitive in the sense that it is highly dependent on how the input is conditioned, or it has a high polynomial dependency on the desired solution accuracy.

In this paper, we focus on analyzing simple versions of IRLS which overcome both aforementioned obstacles. In particular, our methods always converge to $1+\varepsilon$ multiplicative approximation for the objectives $\min_{\boldsymbol{\mathit{x}}:\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\|\boldsymbol{\mathit{x}}\|_{p}$, $p\in\{1,\infty\}$, in $\widetilde{O}(m^{1/3}/\varepsilon^{2/3}+1/\varepsilon^{2})$ iterations111We use $\widetilde{O}$ notation to suppress polylogarithmic factors in $m/\varepsilon$. of solving a weighted least squares problem, where $m$ is the dimension of the sought vector $\boldsymbol{\mathit{x}}$.

Inspiration for our methods is heavily drawn from the work of [CKM*+*11], which offered a ground-breaking result by showing that in undirected graphs, a $(1+\varepsilon)$-approximate maximum flow can be found in $\widetilde{O}(m^{1/3}/\varepsilon^{11/3})$ iterations (subsequently the $\varepsilon$ dependence was improved to $1/\varepsilon^{8/3}$, see [CMMP13]) of solving a weighted least squares problem – which in conjunction with efficient Laplacian system solvers, broke a longstanding barrier for fast graph algorithms. While these algorithms generalize to arbitrary $\ell_{1}$ and $\ell_{\infty}$ regression problems, they are somewhat involved, in particular due to the fact that they are the product of combining the multiplicative weights update method with a regularization technique, and a second potential function.222To be more specific, Christiano et al. solve the approximate maximum flow problem, which is a specific instance of $\ell_{\infty}$ regression. Chin et al. build on this work to solve $\ell_{1}$ regression with block structure; the block structure is relevant for their specific applications, but is a direct extension of the method, so solving vanilla $\ell_{1}$ regression is still the main problem tackled there.

Instead, our method attempts to directly solve the non-smooth objective while tracking a single potential function. The number of iterations looks surprising, since the dominant term is $\widetilde{O}(m^{1/3}/\varepsilon^{2/3})$, whenever $\varepsilon\geq m^{-1/4}$, while classical techniques for optimizing non-smooth functions require a number of iterations that depends on the product between the function’s parameters (such as Lipschitz constant of the gradient or radius of the domain), and $1/\varepsilon$ in the best case, when accelerated methods are used; see [Nes05] for more details.

333We emphasize that using off-the-shelf methods, without further assumptions on the input, the number of iterations of any standard optimization method would be $\Omega(\sqrt{m})$ even for the very special instances where the affine constraint corresponds to a flow satisfying a given demand in unweighted graphs, and in general will depend on how the input matrix is conditioned, since this conditioning determines the magnitude of the subgradients or the diameter of the domain we are optimizing over. The breakthrough of Christiano et al. was the first work that managed to reduce this dependence for maximum flow, which is a specific instance of the $\ell_{\infty}$ regression problem.

Interestingly, a line of work that yielded results very similar in spirit to ours is that of approximately solving positive linear and semidefinite programs [You01, AZO15, AZLO16], where the goal was to produce a first order optimization method that can be implemented in a number of iterations independent of the conditioning of the input. Improving the $\varepsilon$ dependence to $o(1/\varepsilon^{2})$ is an important open problem in this subfield.

We believe that understanding the connection between these results can pave the way for designing new efficient optimization primitives.

1.1 Main Theorem

We state the main theorem of this paper. It follows from the correctness proofs described in Sections 3.1 and 3.2, and the convergence proofs from Lemmas 4.6 and 4.8.

Theorem 1.1**.**

There exist algorithms $\ell_{\infty}$-Minimization and $\ell_{1}$-Minimization such that, on input $(\boldsymbol{\mathit{A}},\boldsymbol{\mathit{b}},\varepsilon,M)$, where $\boldsymbol{\mathit{A}}\in\mathbb{{R}}^{n\times m}$ is a matrix, $\boldsymbol{\mathit{b}}\in\mathbb{{R}}^{n}$ is a vector which lies in the span of $\boldsymbol{\mathit{A}}$’s columns, $\varepsilon$ is an accuracy parameter, and $M$ is a target value:

$\ell_{\infty}$-Minimization* returns a solution $\boldsymbol{\mathit{x}}$ such that $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$, and $\|\boldsymbol{\mathit{x}}\|_{\infty}\leq(1+\varepsilon)M$, or certifies that $\min_{\boldsymbol{\mathit{x}}:\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\|\boldsymbol{\mathit{x}}\|_{\infty}\geq(1-\varepsilon)M$.* 2. 2.

$\ell_{1}$-Minimization* returns a solution $\boldsymbol{\mathit{x}}$ such that $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$, and $\|\boldsymbol{\mathit{x}}\|_{1}\leq(1+\varepsilon)M$, or certifies that $\min_{\boldsymbol{\mathit{x}}:\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\|\boldsymbol{\mathit{x}}\|_{1}\geq(1-\varepsilon)M$.*

Furthermore both algorithms finish in

[TABLE]

iterations, each of which can be implemented in the time required to solve a linear system of the form $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{D}}\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}}=\boldsymbol{\mathit{b}}$, where $\boldsymbol{\mathit{D}}\in\mathbb{{R}}^{m\times m}$ is an arbitrary nonnegative diagonal matrix.

While our theorem statements are concerned with approximately solving a decision problem which requires a guess $M$ on the value of the objective, it follows from standard techniques that this can be used to find a good approximation to the optimal solution without paying more than a $\widetilde{O}(1)$ overhead in the number of iterations. For completeness, we provide the details in Section 6.

1.2 Relation to Previous IRLS Methods and Slime-Mold Dynamics

A popular method for solving $\ell_{1}$ minimization is the iteratively re-weighted least squares method (IRLS). This is essentially based on the observation that whenever $\boldsymbol{\mathit{x}}^{*}=\arg\min_{\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\|\boldsymbol{\mathit{x}}\|_{1}$, one also has that this is the minimizer of the least squares problem $\arg\min_{\boldsymbol{\mathit{x}}:\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\langle{1}/{\boldsymbol{\mathit{x}}^{*}},\boldsymbol{\mathit{x}}^{2}\rangle$.444Throughout the paper we use the convention that $0/0=0$. Hence one approach that has been employed ever since the 60’s [Law61, Osb85, DDFG10] is to iteratively adjust the weighting of the coordinates and re-solve the least squares problem, until $\boldsymbol{\mathit{x}}$ converges to a stationary point. This is rigorously described by the iteration

[TABLE]

We abuse notation by applying scalar operations to vectors, with the meaning that they are applied element-wise.

Subsequent works attempted to rigorously analyze this iteration and prove convergence bounds. Oftentimes this relied on specific structure, such as $\boldsymbol{\mathit{x}}$ being sparse [DDFG10]. A recent series of works drew inspiration from convergence proofs for the slime-mold dynamics – a method which essentially solves $\ell_{1}$ minimization, based on a model used to describe the evolution of a slime mold (Physarum polycephalum) as it spreads through its environment in order to optimize its access to food sources [NYT00, TKN07]. Based on the intuition that these dynamics yield a method for solving the transportation problem, Straszak and Vishnoi proved in a series of works [SV16c, SV16b, SV16a] that this is as a matter of fact equivalent to the IRLS method, and provided a rigorous convergence analysis for a damped version of it:

[TABLE]

Unfortunately their convergence proof shows that this method is highly inefficient, and the time to convergence has a high polynomial dependence in the desired accuracy, and the structure of the linear constraint.

By comparison, what we describe in this work is an IRLS method where the weights are updated according to a thresholding rule. Given a guess $M$ for the optimal value, we perform an iteration equivalent to:

[TABLE]

where $\psi$ is a thresholding operator i.e. $\psi_{b}(u)=u$, if $u\geq b$, and $\psi_{b}(u)=1$ otherwise. Intuitively, this increases the weights $c_{i}$ only for the elements where the corresponding component $x_{i}^{2}/c_{i}$ of the quadratic objective contributes significantly, therefore we want to favor increasing it even more in the future by decreasing the weight $1/c_{i}$ we place on this coordinate.555Another way to think of this is that, ignoring the thresholding operator, the update would simply be $c_{i}^{(t+1)}=(x_{i}^{(t)})^{2}/c_{i}^{(t)}\cdot\gamma$, where $\gamma$ is some normalization factor. What thresholding achieves here is to decide whether the contribution of a particular coordinate to the energy of the system is sufficiently large compared to the contributions of the entire vector $\boldsymbol{\mathit{x}}$.

2 Preliminaries

2.1 Basic Notation

Sets.

We let $\mathbb{{R}}$ be the set of real numbers. For any natural number $n$, we write $[n]:=\{1,\dots,n\}$. We denote by $\Delta_{m}$ the $m$-dimensional simplex i.e. \Delta_{m}=\{\boldsymbol{\mathit{p}}\in\mathbb{{R}}^{m}:\sum_{i=1}^{m}\boldsymbol{\mathit{p}}_{i}=1,\boldsymbol{\mathit{p}}_{i}\geq 0\textnormal{ for all i}\}.

Vectors.

We let $\mathbf{0},\mathbf{1}\in{\mathbb{{R}}}^{n}$ denote the all zeros and all ones vectors, respectively. When it is clear from the context, we apply scalar operations to vectors with the interpretation that they are applied coordinate-wise.

Matrices.

We write matrices in bold. We use $\boldsymbol{\mathit{I}}$ to denote the identity matrix. Given a vector $\boldsymbol{\mathit{x}}$ we let $\mathbf{D}(\boldsymbol{\mathit{x}})$ be the diagonal matrix whose entries are given by $\boldsymbol{\mathit{x}}$. For a symmetric matrix $\boldsymbol{\mathit{A}}$, we let $\boldsymbol{\mathit{A}}^{+}$ be its Moore-Penrose pseudoinverse, i.e. $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{A}}^{+}=\boldsymbol{\mathit{A}}^{+}\boldsymbol{\mathit{A}}=\boldsymbol{\mathit{I}}_{\textnormal{Im}(\boldsymbol{\mathit{A}})}$. The pseudoinverse can be thought of as replacing all the nonzero eigenvalues of $\boldsymbol{\mathit{A}}$ with their reciprocals.

Inner products.

When it is convenient, we use $\langle\cdot,\cdot\rangle$ notation to denote inner products. Given two vectors $\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}$ of equal dimensions, we let $\langle\boldsymbol{\mathit{x}},\boldsymbol{\mathit{y}}\rangle=\boldsymbol{\mathit{x}}^{\top}\boldsymbol{\mathit{y}}$.

Norms.

Given a vector $\boldsymbol{\mathit{x}}$, we denote the $\ell_{p}$ norm of $\boldsymbol{\mathit{x}}$ by $\|\boldsymbol{\mathit{x}}\|_{p}=(\sum x_{i}^{p})^{1/p}$. When the subscript is dropped, we refer to the $\ell_{2}$ norm. From this definition, we can also see that $\|\boldsymbol{\mathit{x}}\|_{\infty}=\max_{i}|x_{i}|$.

2.2 Proof Technique

Let us first understand the idea behind our $\ell_{\infty}$ minimization algorithm. The problem we aim to solve is $\min_{\boldsymbol{\mathit{x}}:\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\|\boldsymbol{\mathit{x}}\|_{\infty}$. Letting $\Delta_{m}$ be the $m$-dimensional unit simplex, we can write our objective equivalently as

[TABLE]

where the second identity follows from Sion’s theorem [Sio58], which allows us to interchange min and max. The quantity between the parentheses has a very natural interpretation, in the case of electrical networks: it is precisely the electrical energy required to route a demand $\boldsymbol{\mathit{b}}$ through an electrical network encoded in $\boldsymbol{\mathit{A}}$. Furthermore, we have an easy way to lower bound how this energy increases whenever resistances are increased, which is a finer quantitative version of Rayleigh’s monotonicity principle. More precisely, we can easily certify a lower bound on the increase in energy determined by increasing a single coordinate of $\boldsymbol{\mathit{r}}$. Using this observation, which we make more precise in Section 4.2, we can identify a set of coordinates of $\boldsymbol{\mathit{r}}$ to increase, guaranteeing that if $\boldsymbol{\mathit{r}}\mkern 2.0mu\vphantom{r}^{\prime}$ is the new vector with perturbed resistances, we have

[TABLE]

for a fixed parameter $M$. In the case when no coordinates of $\boldsymbol{\mathit{r}}$ can be increased, while preserving this property, this yields a certificate that $\boldsymbol{\mathit{r}}$ is as a matter of fact (close to) optimal, and thus we are done (Lemma 4.5). Hence our goal becomes that of guaranteeing that $\|\boldsymbol{\mathit{r}}\|_{1}$ increases very fast. Indeed, since the ”electrical energy” increases at the right rate relative to $\|\boldsymbol{\mathit{r}}\|_{1}$, after the latter has increased sufficiently, we can safely guarantee that $\mathcal{E}_{\boldsymbol{\mathit{r}}}(\boldsymbol{\mathit{b}})/\|\boldsymbol{\mathit{r}}\|_{1}\geq(1-\varepsilon)M$, since the increase in $\|\boldsymbol{\mathit{r}}\|_{1}$ cancels out most of the initial error introduced by starting with a potentially poor solution.

The $\ell_{1}$ minimization algorithm relies on squaring the objective, and then writing it equivalently as

[TABLE]

For the first identity we used the fact that $\|\boldsymbol{\mathit{x}}\|_{1}^{2}=\min_{\boldsymbol{\mathit{c}}\in{\Delta}_{m}}\langle 1/\boldsymbol{\mathit{c}},\boldsymbol{\mathit{x}}^{2}\rangle$, achieved at $\boldsymbol{\mathit{c}}=\boldsymbol{\mathit{x}}/\|\boldsymbol{\mathit{x}}\|_{1}$; see [Owe07, SZ12] for further use of this trick.666Interestingly, this can also be thought of as achieving tightness for reverse Hölder’s inequality whenever we are considering the dual ‘norms’ $\ell_{-1}$ and $\ell_{1/2}$. The second identity follows from joint convexity w.r.t. $\boldsymbol{\mathit{c}}$ and $\boldsymbol{\mathit{x}}$, which can be verified by computing the Hessian of the function in $(\boldsymbol{\mathit{x}},\boldsymbol{\mathit{c}})$. So completely oppositely from the previous case, the objective of our problem becomes minimizing electrical energy with respect to a set of inverse resistances, which we will call conductances. Note that in this case the quantity that is invariant under scaling $\boldsymbol{\mathit{c}}$ by a constant is $\mathcal{E}_{1/\boldsymbol{\mathit{c}}}\cdot\|\boldsymbol{\mathit{c}}\|_{1}$. Therefore, equivalently, our goal will be to find the set of conductances $\boldsymbol{\mathit{c}}\geq 0$ for which $\left(\mathcal{E}_{1/\boldsymbol{\mathit{c}}}\right)^{-1}/\|\boldsymbol{\mathit{c}}\|_{1}\geq\frac{1}{(1+\varepsilon)M}$. Similarly to the $\ell_{\infty}$ case, in this case we make progress by iteratively increasing conductances from $\boldsymbol{\mathit{c}}$ to $\boldsymbol{\mathit{c}}^{\prime}$ in such a way that

[TABLE]

Just as before, we can prove that unless the value of the objective can not be made smaller than $M$, then $\boldsymbol{\mathit{c}}$ can be increased while enforcing this invariant (Lemma 4.7). Hence we can prove fast convergence by arguing that $\|\boldsymbol{\mathit{c}}\|_{1}$ increases very fast.

2.3 Approximate Solutions and Infeasibility Certificates

$\ell_{\infty}$ minimization

We consider the formulation

[TABLE]

for which we seek an approximate solution in the following sense. Given a target value $M$, we aim to find one of the following:

an approximate solution $\boldsymbol{\mathit{x}}$ in the sense that $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$ and $\|\boldsymbol{\mathit{x}}\|_{\infty}\leq(1+\varepsilon)M$, 2. 2.

an approximate infeasibility certificate $\boldsymbol{\mathit{r}}$ in the sense that $\boldsymbol{\mathit{r}}\in\Delta_{m}$ and $\boldsymbol{\mathit{b}}^{\top}(\boldsymbol{\mathit{A}}\mathbf{D}(\boldsymbol{\mathit{r}})^{-1}\boldsymbol{\mathit{A}}^{\top})^{+}\boldsymbol{\mathit{b}}\geq(1-\varepsilon)^{2}M^{2}$.

We prove in Lemma 2.1 that the latter is indeed an infeasibility certificate.

Lemma 2.1**.**

Let $\boldsymbol{\mathit{x}}^{*}$ be the solution to the problem defined in Equation 2.3, and let $\boldsymbol{\mathit{r}}\in\Delta_{m}$. Then $\|\boldsymbol{\mathit{x}}^{*}\|_{\infty}^{2}\geq\boldsymbol{\mathit{b}}^{\top}(\boldsymbol{\mathit{A}}\mathbf{D}(\boldsymbol{\mathit{r}})^{-1}\boldsymbol{\mathit{A}}^{\top})^{+}\boldsymbol{\mathit{b}}$.

Proof.

Using Lemma 4.2 we can write

[TABLE]

which gives us what we needed. ∎

$\ell_{1}$ minimization

We consider the formulation

[TABLE]

for which seek an approximate solution in the following sense. Given a target value $M$, we seek one of the following:

an approximate infeasibility certificate $\boldsymbol{\mathit{\phi}}\in\mathbb{{R}}^{n}$ in the sense that $\frac{\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\phi}}}{\|\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}}\|_{\infty}}\geq(1-\varepsilon)M$, 2. 2.

an approximate feasibility certificate $\boldsymbol{\mathit{c}}$ in the sense that $\boldsymbol{\mathit{c}}\in\Delta_{m}$ and ${\boldsymbol{\mathit{b}}^{\top}(\boldsymbol{\mathit{A}}\mathbf{D}(\boldsymbol{\mathit{c}})\boldsymbol{\mathit{A}}^{\top})^{+}\boldsymbol{\mathit{b}}}\leq(1+\varepsilon)^{2}M^{2}$, which yields an approximately feasible solution $\boldsymbol{\mathit{x}}=\mathbf{D}(\boldsymbol{\mathit{c}})\boldsymbol{\mathit{A}}^{\top}(\boldsymbol{\mathit{A}}\mathbf{D}(\boldsymbol{\mathit{c}})\boldsymbol{\mathit{A}}^{\top})^{+}\boldsymbol{\mathit{b}}$ in the sense that $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$ and $\|\boldsymbol{\mathit{x}}\|_{1}\leq(1+\varepsilon)M$.

The fact that the former is an approximate infeasibility certificate follows from convex duality. Indeed, one can see that the dual of the minimization problem is $\max_{\boldsymbol{\mathit{\phi}}:\|\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}}\|_{\infty}\leq 1}\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\phi}}$, so exhibiting a solution as above implies that the value of this objective is at least $(1-\varepsilon)M$. A proof for the fact that the latter is indeed an approximate feasibility certificate, and that it yields an approximately feasible solution can be found in Lemma 2.2.

Lemma 2.2**.**

Given $\boldsymbol{\mathit{c}}\in\Delta_{m}$, the vector $\boldsymbol{\mathit{x}}=\mathbf{D}(\boldsymbol{\mathit{c}})\boldsymbol{\mathit{A}}^{\top}(\boldsymbol{\mathit{A}}\mathbf{D}(\boldsymbol{\mathit{c}})\boldsymbol{\mathit{A}}^{\top})^{+}\boldsymbol{\mathit{b}}$ satisfies $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$, and $\|\boldsymbol{\mathit{x}}\|_{1}^{2}\leq\boldsymbol{\mathit{b}}^{\top}(\boldsymbol{\mathit{A}}\mathbf{D}(\boldsymbol{\mathit{c}})\boldsymbol{\mathit{A}}^{\top})^{+}\boldsymbol{\mathit{b}}$.

Proof.

The fact that $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$ follows directly by substitution, and using the fact that $\boldsymbol{\mathit{b}}\in\textnormal{Im}(\boldsymbol{\mathit{A}})$. Using Lemma 4.2 and the definition in (4.1) we write

[TABLE]

We can use this identity inside the following upper bound, which we obtain by applying Cauchy-Schwarz:

[TABLE]

This yields our claim. ∎

3 The Algorithms

Having introduced the necessary notation, we can describe our simple IRLS routine. We prove convergence in Section 4.

3.1 The $\ell_{\infty}$ Minimization Algorithm

We first present the algorithm for the $\ell_{\infty}$ version of the problem, since it is the most intuitive. The method attempts to find a weighting of the columns of $\boldsymbol{\mathit{A}}$ i.e. a vector $\boldsymbol{\mathit{r}}\in\mathbb{{R}}^{m}$ for which the corresponding least squares solution has a small $\ell_{\infty}$ norm; more precisely $\|\boldsymbol{\mathit{x}}\|_{\infty}/\|\boldsymbol{\mathit{r}}\|_{1}\leq(1+\varepsilon)M$ for some chosen target value $M$.

Then the weighting is updated via the following simple thresholding rule. Elements for which the corresponding coordinate of the least squares solution $x_{i}$ is below the desired target value are left unchanged. The others are scaled exactly by the amount by which the square of the corresponding coordinate $x_{i}$ violates the desired threshold i.e. $x_{i}^{2}/M^{2}$.

Note that the iteration defined here simply attempts to construct an infeasibility certificate for the problem defined in Equation 2.3. Building the feasible solution involves maintains a solution obtained by uniformly averaging a subset of the iterates $\boldsymbol{\mathit{x}}$ witnessed so far. These are used to return the approximately feasible solution in case the algorithm fails to quickly produce an (approximate) infeasibility certificate. The details referring to how and why we perform this specific set of updates are explained in the convergence proof. The steps involved in building this feasible solution are written in blue. They can be ignored if the goal is simply that of returning a yes/no answer.

Correctness.

We notice that Algorithm 1 has two possible outcomes. Either it returns a primal approximately feasible vector (lines 11and 15), or returns a dual certificate (line 20). In the former case, it is clear from the description of the algorithm that the returned vector is indeed approximately feasible: line 11  returns a uniform average of vectors satisfying the linear constraint with small $\ell_{\infty}$ norm; line 20  returns the $\boldsymbol{\mathit{x}}^{(t)}$ computed within the corresponding iteration, whenever $\boldsymbol{\mathit{\alpha}}^{(t)}=\vec{1}$, i.e. $\|\boldsymbol{\mathit{x}}^{(t)}\|_{\infty}<(1+\varepsilon)M$.

Also, note that in case none of these stopping conditions is triggered, the algorithm returns a dual certificate on line 20  after a finite number of iterations. Indeed, note that every iteration where ${\alpha}_{i}^{(t)}\neq\vec{1}$, at least one element of $\boldsymbol{\mathit{r}}^{(t)}$ gets increased by a factor of at least $(1+\varepsilon)^{2}$, due to way $\boldsymbol{\mathit{\alpha}}^{(t)}$ is defined. Since the algorithm stops when $\|\boldsymbol{\mathit{r}}^{(t)}\|_{1}=1/\varepsilon$, no element of $\boldsymbol{\mathit{r}}$ can be scaled more than $O(\log_{(1+\varepsilon)}(m/\varepsilon))$ times, hence the total number of iterations is very roughly upper bounded by $O(m\log(m/\varepsilon)/\varepsilon)$. We will see in Section 4 that we can prove a much finer upper bound.

Finally, we need to argue that whenever the algorithm returns on line 20, it returns an infeasibility certificate as per Lemma 2.1. We defer the proof to Lemma 4.5 in Section 4.

3.2 The $\ell_{1}$ Minimization Algorithm

The $\ell_{1}$ version is very similar. As a matter of fact, it can be re-derived simply by attempting to solve the convex dual of the problem from (2.3), which is an $\ell_{\infty}$ minimization problem, by using the routine from Figure 1. However, since the reduction requires several, and previous works attempted to solve this directly using various versions of IRLS, we provide a natural iteration which does not involve any reductions.

Correctness.

We notice that Algorithm 2 has two possible outcomes. Either it returns an approximate infeasibility certificate (lines 11  and 15), or returns an approximately feasible solution (line 20).

Let us verify that in the former case the returned vector is indeed an approximate infeasibility certificate. Line 11  returns ${\boldsymbol{\mathit{\Phi}}}^{(t^{\prime})}=\sum_{t\in S}\frac{\boldsymbol{\mathit{\phi}}^{(t)}}{\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)}}$, where we know that $S$ is a set for which

[TABLE]

Since $\boldsymbol{\mathit{b}}^{\top}{\boldsymbol{\mathit{\Phi}}}^{(t^{\prime})}=t^{\prime}$, we see that returned vector ${\boldsymbol{\mathit{\Phi}}}^{(t^{\prime})}/t^{\prime}$ is an approximate infeasibility certificate, as defined in Section 2.3. If the algorithm returns on line 15, we get that $\left\|\frac{\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)}}{\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)}}\right\|_{\infty}\leq\frac{1}{(1-\varepsilon)M}$, hence $\boldsymbol{\mathit{\phi}}^{(t)}$ is an approximate infeasibility certificate.

Also, note that in case none of these stopping conditions is triggered, the algorithm returns a solution on line 20  after a finite number of iterations. Indeed, just as in the $\ell_{\infty}$ case, in every iteration some conductance gets increased by a factor of at least $\Omega(1+\varepsilon)$, hence the algorithm must stop in finite time. We provide a rigorous analysis of the time required for convergence in Section 4.

Finally, we need to argue that whenever the algorithm returns a solution on line 20, it is indeed an approximately feasible solution. We defer the proof to Lemma 4.7 in Section 4.

4 The Algorithm Analyses

4.1 The Flow/Potential Interpretation

While we study a very general problem, it is very useful to develop intuition based on the case where $\boldsymbol{\mathit{A}}$ is the vertex-edge incidence matrix of a graph. In this case we will always think of the sought solution $\boldsymbol{\mathit{x}}$ as a flow on the graph’s edges. The corresponding dual object is a set of potentials $\boldsymbol{\mathit{\phi}}$ defined on the graph’s vertices.

To be more precise, we consider the following setting. Let $G=(V,E)$ be an undirected graph. For each edge we choose an arbitrary orientation, and define $E^{+}(v)$ be the set of arcs leaving vertex $v$, and $E^{-}(v)$ the set of arcs entering vertex $v$, for all $v$.

Letting $m=|E|$, $n=|V|$, we consider the matrix $\boldsymbol{\mathit{A}}\in\mathbb{{R}}^{n\times m}$ where

[TABLE]

One can easily verify that given a vector $\boldsymbol{\mathit{x}}\in\mathbb{{R}}^{m}$ defined on the arcs of the graph (which we will think of as a flow), after applying the operator $\boldsymbol{\mathit{A}}$ we obtain the demand routed by this flow $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}\in\mathbb{{R}}^{n}$, which lives in the space of potentials defined on the graph’s vertices.

Therefore the $\ell_{\infty}$ minimization problem from 2.3 can be interpreted as finding the flow $\boldsymbol{\mathit{x}}$ with minimum congestion which routes the demand $\boldsymbol{\mathit{b}}$, while the $\ell_{1}$ minimization problem from 2.4 corresponds to finding the minimum cost flow routing the demand $\boldsymbol{\mathit{b}}$.

With this interpretation in mind, we proceed to define some objects that in the case of electrical networks correspond to energy and electrical flows.

We use weightings of $\boldsymbol{\mathit{A}}$’s columns $\boldsymbol{\mathit{c}}\in\mathbb{{R}}^{m}$ which we refer to as conductances. We equivalently refer to the reciprocals $\boldsymbol{\mathit{r}}\in\mathbb{{R}}^{m}$, with $\boldsymbol{\mathit{r}}_{i}=1/\boldsymbol{\mathit{c}}_{i}$, which we call resistances. Our analysis is exclusively based on tracking a potential function which corresponds to the electrical energy of a flow.

Definition 4.1** (Energy of a flow).**

Given a flow $\boldsymbol{\mathit{x}}\in\mathbb{{R}}^{m}$, along with a vector of resistances $\boldsymbol{\mathit{r}}\in\mathbb{{R}}^{m}$, we let the energy of $\boldsymbol{\mathit{x}}$ be

[TABLE]

Overloading this notation, given a vector $\boldsymbol{\mathit{b}}\in\mathbb{{R}}^{n}$, let the electrical energy be

[TABLE]

in other words this is the minimum energy over all flows satisfying $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$. We drop the argument whenever $\boldsymbol{\mathit{b}}$ is clear from the context.

4.2 Preliminaries on Electrical Energy

Throughout the paper, our analyses will rely on a potential function, which in the case of resistor networks corresponds to the electrical energy. In this section we provide a few useful facts.

Lemma 4.2** (Characterization of Electrical Energy).**

Given a vector of resistances $\boldsymbol{\mathit{r}}\in\mathbb{{R}}^{m}$, we have the following equivalent characterizations for the electrical energy.

[TABLE]

Furthermore, if $\boldsymbol{\mathit{x}}$ is the minimizing flow for the expression in (4.1), and $\boldsymbol{\mathit{\phi}}$ is the maximizing set of potentials for the expression in (4.3), then for all $i$:

[TABLE]

Since the proof is standard, we defer it to Section A.1.

As a corollary, we can derive a lower bound on the increase in energy after increasing resistances.

Lemma 4.3**.**

Let $\boldsymbol{\mathit{r}},\boldsymbol{\mathit{r}}\mkern 2.0mu\vphantom{r}^{\prime}$, and let $\boldsymbol{\mathit{x}}=\arg\min_{\boldsymbol{\mathit{x}}:\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\langle\boldsymbol{\mathit{r}},\boldsymbol{\mathit{x}}^{2}\rangle$. Then, one has that

[TABLE]

Proof.

We use the characterization from Equation 4.3 for characterizing electrical energy. Let $\boldsymbol{\mathit{\phi}}$ be the argument that maximizes $(\ref{eq:energy_max_char})$ for resistances $\boldsymbol{\mathit{r}}$. We certify a lower bound on $\mathcal{E}_{\boldsymbol{\mathit{r}}\mkern 2.0mu\vphantom{r}^{\prime}}(\boldsymbol{\mathit{b}})$ using $\boldsymbol{\mathit{\phi}}$ as follows:

[TABLE]

Finally substituting the relation between flows and potentials from Lemma 4.2, Equation (4.5), we obtain the desired claim. ∎

We can derive a similar lower bound on the inverse energy, after increasing conductances.

Lemma 4.4**.**

Let $\boldsymbol{\mathit{\phi}}=\arg\min_{\phi:\langle\boldsymbol{\mathit{b}},\boldsymbol{\mathit{\phi}}\rangle=1}\langle\boldsymbol{\mathit{c}},(\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}})^{2}\rangle$. Then one has that

[TABLE]

Proof.

We use the following basic inequality: for $x,x^{\prime}>0$ one has $\frac{1}{x^{\prime}}\geq\frac{1}{x}+\frac{x-x^{\prime}}{x^{2}}$, which follows from $(x-x^{\prime})^{2}\geq 0$. Also, from the definition of energy in (4.1), we obtain an upper bound on the new energy, after perturbing conductances. Let $\boldsymbol{\mathit{x}}=\arg\min_{\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}}\langle 1/\boldsymbol{\mathit{c}},\boldsymbol{\mathit{x}}^{2}\rangle$, i.e. the electrical flow corresponding to conductances $\boldsymbol{\mathit{c}}$. We therefore have:

[TABLE]

Using the fact that by optimality, $x_{i}=c_{i}(\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}})_{i}$ (per Lemma 4.2), and combining with the previous inequality we obtain

[TABLE]

which is what we wanted. ∎

4.3 Convergence Proof for $\ell_{\infty}$ Minimization

Having put together all these tools, we are ready to analyze the algorithms presented in Section 3. We first prove that $\ell_{\infty}$-Minimization returns a correct infeasibility certificate, whenever it returns on line 20. This lemma is key to understanding the intuition behind the algorithm.

Lemma 4.5**.**

Whenever $\ell_{\infty}$-Minimization returns on line 20, $\boldsymbol{\mathit{r}}/\|\boldsymbol{\mathit{r}}\|_{1}$ is a correct approximate infeasibility certificate in the sense that

[TABLE]

Proof.

First notice that by Lemma 2.1, the lower bound on energy is indeed an approximate infeasibility certificate. Now we proceed to prove that throughout the iterations of the algorithm, energy increases at the right rate.

We show that every iteration satisfies the invariant

[TABLE]

This is easy to verify using Lemma 4.3, which lower bounds the increase in energy after perturbing resistances. We see that using the perturbation rule defined on line 13  of the algorithm, energy increases as follows

[TABLE]

For every coordinate of $\boldsymbol{\mathit{r}}^{(t)}$ that has changed we see that the ratio between the contribution to above lower bound of that specific coordinate, and the increase in resistance is

[TABLE]

Therefore, summing up over all coordinates we obtain the desired inequality. Finally, we notice that initially $\mathcal{E}_{\boldsymbol{\mathit{r}}^{(0)}}(\boldsymbol{\mathit{d}})\geq 0$, and $\|\boldsymbol{\mathit{r}}^{(0)}\|_{1}=1$. So once $\|\boldsymbol{\mathit{r}}^{(t)}\|_{1}\geq\frac{1}{\varepsilon}$, one has that, using (4.6),

[TABLE]

and thus

[TABLE]

and equivalently:

[TABLE]

which implies what we needed. ∎

Knowing that the algorithm is correct, we can now proceed and prove that it converges fast (convergence rate can be slightly improved by using a more careful schedule for $M$ and $\varepsilon$; we defer this improvement to Section 5).

Lemma 4.6**.**

The algorithm $\ell_{\infty}$-Minimization returns a solution after $O(m^{1/3}\log(1/\varepsilon)/\varepsilon+\log(m/\varepsilon)/\varepsilon^{2})$ iterations.

Proof.

We show that unless the algorithm returns an approximately feasible solution on lines 11  or 15, then there exists a coordinate $i\in[m]$ for which $r_{i}$ increases very fast.

Suppose the algorithm has run for $T$ iterations without returning an approximately feasible solution. Consider the partial sum of iterates obtained so far $\boldsymbol{\mathit{s}}^{(t^{\prime})}$ for some $t^{\prime}\leq T$. Since the algorithm did not return on line 11, we know that $\|\boldsymbol{\mathit{s}}^{(t^{\prime})}\|_{\infty}/t^{\prime}\geq(1+\varepsilon)M$. Therefore there exists a coordinate $i\in[m]$ for which $\boldsymbol{\mathit{s}}^{(t^{\prime})}_{i}\geq(1+\varepsilon)Mt^{\prime}$. In other words, letting $I$ be the set of iterates that have contributed to $\boldsymbol{\mathit{s}}^{(t^{\prime})}$, one definitely has that

[TABLE]

and thus

[TABLE]

where we used the fact that for each iteration $t\in I$ one has that $\sqrt{\alpha_{i}^{(t)}}=|x_{i}^{(t)}|M$ due to the perturbation rule defined on line 13. This implies that restricting ourselves only to iterations where $\alpha_{i}$ increased the corresponding resistance $r_{i}$, we have that

[TABLE]

By the condition on line 7  we see that for all iterations $t\in I$, one has

[TABLE]

Also since we only consider the iterations $t\in I$ with $\alpha_{i}^{(t)}>1$, the rule from line 13  also enforces that for all these iterations

[TABLE]

Equations (4.7), (4.8) and (4.9) suggest that the product $\prod_{t\in I,\alpha_{i}^{(t)}>1}\sqrt{\alpha_{i}^{(t)}}$ increases very fast: intuitively the worst case should occur either when all the factors contribute equally, either all of them are as small as possible (i.e. $1+\varepsilon$, or as large as possible, i.e. $m^{1/3}$). We formalize this intuition in Lemma A.1, which implies that

[TABLE]

Hence setting

[TABLE]

suffices to lower bound this product by $\sqrt{m/\varepsilon}$. Since each iteration a resistance ${r}_{i}^{(t)}$ gets multiplied by the corresponding $\alpha_{i}^{(t)}$, and all resistances are initially $1/m$, this lower bound implies that $r^{(t)}_{i}\geq 1/\varepsilon$. But this means that the algorithm will finish execution after the current iteration, according to the condition on line 5.

Finally, we need to upper bound the number of iterations not in $I$; these correspond to those iterations where $\|\boldsymbol{\mathit{x}}^{(t)}\|_{\infty}\geq m^{1/3}\cdot M$, so there exists some index $i$ for which $\alpha^{(t)}_{i}\geq m^{2/3}$. Therefore some resistance gets multiplied by $m^{2/3}$. Since all resistances are initially $1/m$, in the worst case, each such iteration increases one resistance from $1/m$ to $m^{-1/3}$. Therefore this can happen at most $m^{1/3}\log(1/\varepsilon)/\varepsilon$ times, before the sum of resistances becomes at least $1/\varepsilon$, and the algorithm finishes.

Combining these two cases, we obtain our bound. ∎

We can prove the convergence bound for $\ell_{1}$-Minimization similarly. The main difference is that this time we maintain conductances, and the potential function that enables us to prove convergence is $1/\mathcal{E}_{1/\boldsymbol{\mathit{c}}}$.

4.4 Convergence Proof for $\ell_{1}$-Minimization

Lemma 4.7**.**

Whenever $\ell_{1}$-Minimization returns on line 20, $\boldsymbol{\mathit{c}}/\|\boldsymbol{\mathit{c}}\|_{1}$ is a correct approximate feasibility certificate in the sense that

[TABLE]

Proof.

By Lemma 2.2, this also yields a solution $\boldsymbol{\mathit{x}}$ such that $\boldsymbol{\mathit{A}}\boldsymbol{\mathit{x}}=\boldsymbol{\mathit{b}}$ and $\|\boldsymbol{\mathit{x}}\|_{1}\leq\sqrt{\mathcal{E}_{\|\boldsymbol{\mathit{x}}\|_{1}/\boldsymbol{\mathit{c}}}}\leq M(1+\varepsilon)$.

In order to prove that at the end of the execution the $\ell_{1}$ norm of this solution is small enough, this time we track as potential function the inverse energy $1/\mathcal{E}_{1/\boldsymbol{\mathit{c}}}$. More precisely, we show that every iteration satisfies the invariant

[TABLE]

This is easy to verify using Lemma 4.4, which lower bounds the increase in inverse energy after perturbing conductances. We see that using the perturbation rule defined on line 13  of the algorithm, inverse energy increases as follows

[TABLE]

For every coordinate of $\boldsymbol{\mathit{c}}^{(t)}$ that has changed we see that the ratio between the contribution to above lower bound of that specific coordinate, and the increase in conductance is

[TABLE]

where we used the fact that $\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)}=\mathcal{E}_{1/\boldsymbol{\mathit{c}}^{(t)}}$ (Lemma 4.2).

Therefore, summing up over all coordinates we obtain the desired inequality. Since $\mathcal{E}_{1/\boldsymbol{\mathit{c}}^{(0)}}\geq 0$ and $\|\boldsymbol{\mathit{c}}^{(0)}\|_{1}=1$, we know that once $\|\boldsymbol{\mathit{c}}^{(t)}\|_{1}\geq 1+\frac{1}{(1+\varepsilon)^{2}-1}=O(\frac{1}{\varepsilon})$, one has that, using (4.10),

[TABLE]

and thus,

[TABLE]

and equivalently:

[TABLE]

which is what we needed. ∎

Next we prove that the algorithm converges fast. Convergence rate can be slightly improved by using a more careful schedule for $M$ and $\varepsilon$, which we defer to Section 5.

Lemma 4.8**.**

The algorithm $\ell_{1}$-Minimization returns a solution after $O(m^{1/3}\log(1/\varepsilon)/\varepsilon+\log(m/\varepsilon)/\varepsilon^{2})$ iterations.

Proof.

The proof follows the lines of the proof we used for Lemma 4.6: unless the algorithm returns an approximate infeasibility certificate on lines 11  or 15, then there exists a coordinate $i\in[m]$ for which $c_{i}$ increases very fast.

Suppose the algorithm has run for $T$ iterations without returning an approximate infeasibility certificate. Consider the partial sum of iterates obtained so far $\boldsymbol{\mathit{s}}^{(t^{\prime})}$ for some $t^{\prime}\leq T$. Since the algorithm did not return on line 11, we know that $\|\boldsymbol{\mathit{s}}^{(t^{\prime})}\|_{\infty}/t^{\prime}\geq\frac{1}{(1-\varepsilon)M}$, therefore there exists a coordinate $i\in[m]$ for which $s^{(t^{\prime})}_{i}\geq t^{\prime}\cdot\frac{1}{(1-\varepsilon)M}$. In other words, letting $I$ be the set of iterates that contributed to $\boldsymbol{\mathit{s}}^{(t^{\prime})}$, one has that

[TABLE]

and thus, since by definition $\alpha_{i}^{(t)}=\left(\frac{(\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)})_{i}}{\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)}}\right)^{2}\cdot M^{2}$,

[TABLE]

Therefore, restricting ourselves only to iterations where $\alpha_{i}$ increased the corresponding conductance $c_{i}$, we have that

[TABLE]

By the condition on line 7  we see that for all iterations $t\in I$, one has

[TABLE]

So considering only the iterations $t\in I$ with $\alpha_{i}^{(t)}>1$, the rule from line 13  also enforces that for all these iterations

[TABLE]

Combining Equations (4.11), (4.12), and (4.13), and applying Lemma A.1, exactly the same way we did in the proof of Lemma 4.6 implies that

[TABLE]

So if

[TABLE]

once again we have that this product is lower bounded by $\sqrt{m\cdot\left(1+\frac{1}{(1+\varepsilon)^{2}-1}\right)}$, therefore the corresponding conductance $c_{i}^{(T)}\geq 1+\frac{1}{(1+\varepsilon)^{2}-1}$, since its initial value was $1/m$. Since we can only control the total number of iterations $T$, we can lower bound $t^{\prime}$ by showing that the number of iterations not in $I$ can not be too large. Just as before, we lower bound the number of iterations where $\left\|\frac{\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)}}{\boldsymbol{\mathit{b}}^{\top}\boldsymbol{\mathit{\phi}}^{(t)}}\right\|_{\infty}\geq m^{1/3}/M$. Note that whenever this happens, there exists an index $i$ for which $\alpha_{i}^{(t)}\geq m^{2/3}$. Therefore some conductance gets multiplied by $m^{2/3}$. Again, using an identical argument to the one from the proof of Lemma 4.6, we see that this can not happen more than $O(m^{1/3}\log(1/\varepsilon)/\varepsilon)$ times. Combining this with the sufficient number of iterations required by the other case, we obtain our bound. ∎

5 Using Phases to Improve the Iteration Count

In this section, we show that via minor modifications to our algorithms, we can improve the number of iterations to $O\left(\frac{m^{1/3}\log(1/\varepsilon)}{\varepsilon^{2/3}}+\frac{\log m}{\varepsilon^{2}}\right)$ thus obtaining the bound promised by Theorem 1.1. This relies on the observation that the entire difficulty of the problem is concentrated on improving the quality of a solution from $(1-2\varepsilon)M$ to $(1-\varepsilon)M$. For conciseness, let us focus on the $\ell_{\infty}$ case, and consider the convergence argument described in Sections 4.3. Our goal there is to increase the sum of resistances to $1/\varepsilon$, since our argument assumes that the initial energy could be arbitrarily small.

However, if we assume that we warm start the method with a set of resistances $\boldsymbol{\mathit{r}}_{0}$, $\|\boldsymbol{\mathit{r}}_{0}\|_{1}=1$, for which the corresponding energy is already large enough, $\mathcal{E}_{\boldsymbol{\mathit{r}}_{0}}\geq(1-2\varepsilon)^{2}M^{2}$, we only need to iterate until we obtain a set of resistances $\boldsymbol{\mathit{r}}$ such that $\|\boldsymbol{\mathit{r}}\|_{1}=3$ (rather than $1/\varepsilon$) in order to certify that the current energy/resistance ratio is as large as desired, i.e. $\mathcal{E}_{\boldsymbol{\mathit{r}}}/\|\boldsymbol{\mathit{r}}\|_{1}\geq(1-\varepsilon)^{2}M^{2}$. This in turn improves the number of iterations the algorithm needs before it returns. We expand these ideas in what follows.

Now, suppose we have a set of resistances $\boldsymbol{\mathit{r}}_{0}$, such that $\|\boldsymbol{\mathit{r}}_{0}\|_{1}=1$ and $\mathcal{E}_{\boldsymbol{\mathit{r}}_{0}}\geq(1-2\varepsilon)^{2}M^{2}$. Let us analyze the number of iterations of the method described in Section 4.3 that we require before we can return $\boldsymbol{\mathit{r}}$ such that $\mathcal{E}_{\boldsymbol{\mathit{r}}}/\|\boldsymbol{\mathit{r}}\|_{1}\geq(1-\varepsilon)^{2}M^{2}$ or a solution $\boldsymbol{\mathit{x}}$ such that $\|\boldsymbol{\mathit{x}}\|_{\infty}\leq(1+\varepsilon)M$.

First, we claim that if each update satisfies the invariant from Equation (4.6), then we can stop iteration once $\|\boldsymbol{\mathit{r}}\|_{1}=3$. Indeed, in this case, one has that

[TABLE]

whenever $\|\boldsymbol{\mathit{r}}-\boldsymbol{\mathit{r}}_{0}\|_{1}\geq 3$.

The remaining analysis is carried over almost identically, except that the threshold set on line 7  is changed to $\rho=\varepsilon m$, and our goal is to get $\prod_{t\in I,\alpha_{i}^{(t)}>1}\alpha_{i}\geq\sqrt{3m}$.

For the iterations that satisfy this threshold, by applying Lemma A.1 we see that it is sufficient to witness a small number $t^{\prime}$ of such iterations such that

[TABLE]

so $t^{\prime}=\Theta\left(\frac{\rho}{\varepsilon}\cdot\frac{\log m}{\log\rho}+\frac{\log m}{\varepsilon^{2}}\right)$ suffices.

For the iterations that do not satisfy the threshold, in the worst case, each of them increases one resistance from $1/m$ to $\rho^{2}/m$ so this can happen at most $O(m/\rho^{2})$ times.

Setting $\rho=\left(\varepsilon m\right)^{1/3}$ we get that the total number of iterations is at most $O\left(\frac{m^{1/3}}{\varepsilon^{2/3}}\log(1/\varepsilon)+\frac{\log m}{\varepsilon^{2}}\right)$.

All of this holds assuming we have a good warm start for resistances. We obtain it by recursively invoking the same method for target value $(1-1.75\varepsilon)^{2}M^{2}$, and $.25\varepsilon$ accuracy. In case of failure, this returns a vector $\boldsymbol{\mathit{x}}$ which certainly satisfies $\|\boldsymbol{\mathit{x}}\|_{\infty}\leq M$, so this concludes the entire run on the algorithm; otherwise, it returns a certificate consisting of resistances for which the ratio between the corresponding energy and $\ell_{1}$ norm is at least $(1-2\varepsilon)^{2}M^{2}$, so they can be used as a warm start.

Recursion ends once $\varepsilon\geq 1/2$. We note that since the desired accuracy gets increased by a constant factor after each level of recursion, the total number of iterations is dominated by those performed at the top level (i.e. for the lowest $\varepsilon$). Hence our result.

Note that this method can also be implemented slightly more naturally by running Algorithm 1 with a varying schedule for $M$ and $\varepsilon$.

Improving the number of iterations for $\ell_{1}$ minimization is done analogously.

6 From Approximate Decision to Approximate Optimization

Our algorithms are designed to solve an approximate decision problem, given a guess for the value of the optimal solution. While this follows from a standard reduction, for the sake of completeness we prove here that this is sufficient to optimize the problem approximately without paying more than an additional constant overhead in running time.

To be more specific, let us first focus on $\ell_{\infty}$ minimization. Theorem 1.1 states that given a guess $M$ and accuracy $\varepsilon$, the algorithm either returns an approximately feasible solution with value $\|\boldsymbol{\mathit{x}}\|_{\infty}\leq(1+\varepsilon)M$, or an infeasibility certificate certifying that $\|\boldsymbol{\mathit{x}}^{*}\|_{\infty}\geq(1-\varepsilon)M$. Hence this restricts the search interval for the true value either within the interval $[0,(1+\varepsilon)M]$ or $[(1-\varepsilon)M,\infty)$.

We initialize our search interval to $[\|\boldsymbol{\mathit{x}}_{0}\|_{2}/\sqrt{m},\|\boldsymbol{\mathit{x}}_{0}\|_{\infty}]$ where $\boldsymbol{\mathit{x}}_{0}$ is the initial iterate obtained with uniform resistances. Using Lemma 2.1 we easily verify that $\|\boldsymbol{\mathit{x}}_{0}\|_{2}/\sqrt{m}$ is indeed a lower bound on $\|\boldsymbol{\mathit{x}}^{*}\|_{\infty}$, since energy lower bounds the squared optimal value.

Given a search interval $[L,U]$, we let $M=\sqrt{LU}$, $\tilde{\varepsilon}=\min\left\{\frac{1}{2},\left(\frac{U}{L}\right)^{1/6}-1\right\}$. We invoke Theorem 1.1 for target value $M$ and accuracy $\tilde{\varepsilon}$. Depending on the outcome we update the search interval to $[L,(1+\tilde{\varepsilon})M]$ or $[(1-\tilde{\varepsilon})M,U]$.

When $U/L\leq 1+\varepsilon/4$ we stop the search, call the algorithm for target value $U(1+\varepsilon/4)$ and accuracy $\frac{\varepsilon/4}{1+\varepsilon/4}$, then output the approximately feasible iterate returned by the algorithm. The fact that this call indeed returns an approximately feasible iterate follows from the fact that $U$ is certainly feasible, since this is an invariant maintained by our search, and that if the algorithm were to return an infeasibility certificate it must have needed that $U(1+\varepsilon/4)(1-\frac{\varepsilon/4}{1+\varepsilon/4})<U$, which is false. Thus we know that the returned solution has value at most $U(1+\varepsilon/4)(1+\frac{\varepsilon/4}{1+\varepsilon/4})\leq L(1+\varepsilon)$, so it satisfies the desired approximation guarantee.

Finally, we analyze the cost of the search. Note that each iteration of the search reduces $\log U-\log L$ be a constant factor, and it stops whenever it becomes at most $\log(1+\varepsilon/4)=\Theta(\varepsilon)$. For as long as $U/L>(3/2)^{6}$, the algorithm is invoked with accuracy $1/2$, and $\log U-\log L$ gets reduced by a constant factor, so this happens at most $O(\log\log m)$ times. Once $U/L$ becomes small enough, i.e. $\log U-\log L<6\log(3/2)$, we use accuracy $\exp((\log U-\log L)/6)-1=\Theta(\log(U/L))$. Note that from Theorem 1.1 we know that the number of iterations of the algorithm for a single invocation depends on $1/\tilde{\varepsilon}^{c}$, where $c$ is a fixed constant; due to our schedule for choosing $\tilde{\varepsilon}$, the total cost of this sequence of invocations is dominated by the cost of the final one, where $\tilde{\varepsilon}=\Theta(\varepsilon)$.

So letting $\mathcal{T}(\varepsilon)$ be the number of iterations required by the algorithm from Theorem 1.1 to solve the approximate decision problem to accuracy $\varepsilon$, we have that solving the approximate optimization problem requires $O\left(\mathcal{T}(1/2)\log\log m+\mathcal{T}(\varepsilon)\right)$ iterations.

The $\ell_{1}$ minimization problem is treated similarly, so we omit its description.

7 Experiments

We test both our resistance/conductance update schemes in order to verify that the resulting algorithms converge fast in practice. We slightly modify the schemes such that they always update their target value $M$ depending on the value of the objective they have achieved so far. We stop when given the history of witnessed iterates, we can certify a sufficiently small duality gap. For solving linear systems, we used the conjugate gradient implementation from the $\ell_{1}$-MAGIC optimization suite [CR].

We test both algorithms while varying $\varepsilon$, and varying $m$. We consider both the update scheme given by our algorithms from Section 3, and one where we attempt to double the length of the step for as long as the invariants from (2.1) and (2.2), respectively, are maintained. We notice that in general, using this long step strategy, we improve both the number of iterations and the running time.

The plots corresponding to the standard update scheme (short-steps) are drawn in red, those corresponding to the long-step version are drawn in blue.

The experiments where we varied $\varepsilon$ are reported in figures 1(a), 1(b), 1(e), and  1(f). For all these experiments, the input consists a random $150\times 200$ matrix $\boldsymbol{\mathit{A}}$ with orthogonal rows, and a vector $\boldsymbol{\mathit{b}}$ obtained from applying $\boldsymbol{\mathit{A}}$ to a $\pm 1$-vector of sparsity $15$. We plot the number of iterations/running time of the algorithm for $\varepsilon=1/2^{k}$, where $k\in\{1,\dots,12\}$.

We notice that for these experiments, the number of iterations for the short-step version does indeed scale linearly with $\varepsilon^{-1}$; the long-step version makes significant gains in the $\ell_{\infty}$ case.

The experiments where we varied $m$ are reported in figures 1(c), 1(d), 1(g), and 1(h). For all these experiments, the input consists of a random $150\times(200\cdot k)$ matrix $\boldsymbol{\mathit{A}}$ with orthogonal vectors, and a vector $\boldsymbol{\mathit{b}}$ obtained from applying $\boldsymbol{\mathit{A}}$ to a $\pm 1$-vector of sparsity $15$, and a fixed accuracy $\varepsilon=.01$. We plot the number of iterations required by the algorithm for $k\in\{1,\dots,30\}$.

We notice that for these experiments, both the number of iterations and the running time scale significantly better than by $m^{1/3}$, which suggests that this polynomial dependence in $m$ depends on the input structure, and can be avoided in practice.

Acknowledgements

AE was partially supported by NSF CAREER grant CCF-1750333 and NSF grant CCF-1718342. AV was partially supported by NSF grant CCF-1718342.

Appendix A Deferred Proofs

A.1 Proof of Lemma 4.2

Proof.

We can write the formulation from (4.1) as an unconstrained optimization problem using Lagrange multipliers:

[TABLE]

By making the gradient with respect to $\boldsymbol{\mathit{x}}$ equal to [math], we see that the inner minimization problem is optimized at $2r_{i}\cdot x_{i}=2(\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}})_{i}$ for all $i$, and equivalently $x_{i}=(\boldsymbol{\mathit{A}}^{\top}\boldsymbol{\mathit{\phi}})_{i}/r_{i}$. Plugging this back into the maximization objective w.r.t. $\boldsymbol{\mathit{\phi}}$ we obtain

[TABLE]

where for the last equality we used that by optimality conditions one must have $(\boldsymbol{\mathit{A}}\mathbf{D}(\boldsymbol{\mathit{r}})^{-1}\boldsymbol{\mathit{A}}^{\top})\boldsymbol{\mathit{\phi}}=\boldsymbol{\mathit{b}}$.

Finally, we prove (4.4) by using the fact that for any symmetric matrix $\boldsymbol{\mathit{L}}$ and vector $\boldsymbol{\mathit{b}}$ one has that

[TABLE]

which can be seen by observing that both expressions are optimized at $\boldsymbol{\mathit{\phi}}=\boldsymbol{\mathit{L}}^{+}\boldsymbol{\mathit{b}}$, then applying (4.3).

∎

A.2 Lower Bound Lemma

Lemma A.1**.**

Let a set of nonnegative reals $\beta_{1},\dots,\beta_{k}$ such that $1+\varepsilon\leq\beta_{i}\leq\rho$ for all $i$, and $\sum_{i=1}^{k}\beta_{i}\geq S$. Then, for any $k$, one has that

[TABLE]

Proof.

Consider a fixed $k$, and let us attempt to minimize the product of $\beta_{i}$’s subject to the constraints. Equivalently we want to minimize $\sum_{i=1}^{k}\log(\beta_{i})$, which is a concave function. Therefore its minimizer is attained on the boundary of the feasible domain. This means that for some $0\leq k^{\prime}\leq k-1$, there are $k^{\prime}$ elements equal to $1+\varepsilon$, $k-1-k^{\prime}$ equal to $\rho$, and one which is exactly equal to the remaining budget, i.e. $S-k^{\prime}(1+\varepsilon)-(k-1-k^{\prime})\rho$, which yields the product $(1+\varepsilon)^{k^{\prime}}\rho^{k-k^{\prime}-1}(S-k^{\prime}(1+\varepsilon)-(k-1-k^{\prime})\rho)$. This can be relaxed by allowing $k$ and $k^{\prime}$ to be non-integral. Hence we aim to minimize the product $(1+\varepsilon)^{k^{\prime}}\rho^{k-k^{\prime}}$, subject to $(1+\varepsilon)k^{\prime}+\rho(k-k^{\prime})=S$. Finally, we observe that we can always obtain a better solution by placing all the available mass on a single one of the factors, i.e. we lower bound either by $(1+\varepsilon)^{S/(1+\varepsilon)}$, or $\rho^{S/\rho}$, whichever is lowest. ∎