A Laplacian Approach to $\ell_1$-Norm Minimization

TL;DR

This paper introduces a differentiable Laplacian-inspired reformulation of the basis pursuit problem, enabling new gradient-based methods and analyzing their iteration complexities for solving $ ext{l}_1$ minimization tasks.

Contribution

It presents a novel Laplacian-based reformulation of $ ext{l}_1$ minimization that facilitates gradient methods and provides complexity bounds for iterative algorithms.

Findings

Derived bounds on iteration complexity for multiplicative weights update scheme.

Analyzed accelerated gradient scheme complexity for the reformulated problem.

Connected the approach to IRLS methods, offering new theoretical insights.

Abstract

We propose a novel differentiable reformulation of the linearly-constrained $\ell_1$ minimization problem, also known as the basis pursuit problem. The reformulation is inspired by the Laplacian paradigm of network theory and leads to a new family of gradient-based methods for the solution of $\ell_1$ minimization problems. We analyze the iteration complexity of a natural solution approach to the reformulation, based on a multiplicative weights update scheme, as well as the iteration complexity of an accelerated gradient scheme. The results can be seen as bounds on the complexity of iteratively reweighted least squares (IRLS) type methods of basis pursuit.

Tables

Table 1. Table 1. Worst-case iteration complexity of recent IRLS methods for ℓ 1 subscript ℓ 1 \ell_{1} -norm minimization

Algorithm	Iteration complexity
Ref. [19]	$\tilde{O}\left(m^{1/3}{ϵ}^{-8/3}\right)$
PGS – Theorem 5.4	$\tilde{O}\left(m^2{ϵ}^{-3}\right)$
AGS – Theorem 5.7	$O\left(m^2{ϵ}^{-2}\right)$
Ref. [22]	$\tilde{O}\left(m^{1/3}{ϵ}^{-2/3}+{ϵ}^{-2}\right)$

Equations

$${\Vert s\Vert }_1$$

$$As=b,s\in {\mathbb{R}}^m.$$

$$s_{2j-1}^{\prime }=\{\begin{array}{ll}{u}_j/2& \text{ if }{u}_j\ne 0,\\ +1& \text{ if }{u}_j=0,\end{array},s_{2j}^{\prime }=\{\begin{array}{ll}{u}_j/2& \text{ if }{u}_j\ne 0,\\ -1& \text{ if }{u}_j=0.\end{array}j=1,\dots ,m.$$

$$f_0(x)$$

$$f(x)$$

$$\mathrm{dom}f\stackrel{\mathrm{def}}{=}\{x\in {\mathbb{R}}^m:f(x)<+\infty \}.$$

$$\frac{\partial L^{-1}}{\partial x_j}=-L^{-1}\frac{\partial L}{\partial x_j}{L}^{-1}.$$

$$\frac{\partial b^{\top }{L}^{-1}b}{\partial x_j}=-b^{\top }{L}^{-1}\frac{\partial L}{\partial x_j}{L}^{-1}b=-b^{\top }{L}^{-1}{a}_j{a}_j^{\top }{L}^{-1}b=-(a_j^{\top }{L}^{-1}b{)}^2.$$

$$[{\nabla }^{2}f(x){]}_{ij}=2(b^{\top }{L}^{-1}{a}_i{a}_i^{\top }{L}^{-1}{a}_j{a}_j^{\top }{L}^{-1}b)=2d_i(x)d_j(x)a_i^{\top }{L}^{-1}{a}_j.$$

$${\lambda }_k(SM{S}^{\top })={\theta }_k{\lambda }_k(M).$$

$$\Pi (x)=(A{X}^{1/2}{)}^{\top }(AX{A}^{\top }{)}^{-1}(A{X}^{1/2}).$$

$${\Vert d(x)\Vert }_{\infty }\le (\underset{i=1,\dots ,m}{\min }{x}_i{)}^{-1}(c_{A,b}{)}^{1/2}.$$

$$c_{A,b}^{1/2}\le {\Vert s^{\ast }\Vert }_1\le (m\cdot c_{A,b}{)}^{1/2}.$$

$${\Vert d(x)\Vert }_{\infty }={\Vert T(x)s\Vert }_{\infty }\le {\Vert T(x)s\Vert }_2.$$

$$c_{A,b}={\Vert u\Vert }_2^2\le {\Vert s^{\ast }\Vert }_2^2\le {\Vert s^{\ast }\Vert }_1^2\le {\Vert u\Vert }_1^2\le m{\Vert u\Vert }_2^2=m\cdot c_{A,b}.\text{qed}$$

$${\Vert \nabla f(x)\Vert }_{\infty }\le 1+(\underset{i=1,\dots ,m}{\min }{x}_i{)}^{-2}{c}_{A,b}.$$

$${\lambda }_{\max }({\nabla }^{2}f(x))\le 2(\underset{i=1,\dots ,m}{\min }{x}_i{)}^{-3}\cdot c_{A,b}.$$

$${\nabla }^{2}f(x)=2D(x)T(x)D(x),$$

$${\lambda }_{\max }({\nabla }^{2}f(x))=2\theta {\lambda }_{\max }(T(x))$$

$$\theta \le {\lambda }_{\max }(D(x{)}^2)={\Vert d(x)\Vert }_{\infty }^2\le (\underset{i}{\min }{x}_i{)}^{-2}{c}_{A,b},$$

$${\Vert s\Vert }_1$$

$$As=b,s\in {\mathbb{R}}^m.$$

$$\frac{1}{2}{1}^{\top }x+\frac{1}{2}{b}^{\top }(AX{A}^{\top }{)}^{-1}b$$

$$x\in {\mathbb{R}}_{>0}^m.$$

$$s^{\top }{X}^{-1}s$$

$$As=b,s\in {\mathbb{R}}^m.$$

$$s^{\top }{X}^{-1}s=(q+r{)}^{\top }{X}^{-1}(q+r)=q^{\top }{X}^{-1}q+2r^{\top }{X}^{-1}q+r^{\top }{X}^{-1}r\ge q^{\top }{X}^{-1}q,$$

$$q^{\top }{X}^{-1}q=(b^{\top }{L}^{-1}{A}^{\top }X)X^{-1}(X{A}^{\top }{L}^{-1}b)=b^{\top }{L}^{-1}L{L}^{-1}b=b^{\top }{L}^{-1}b.\text{qed}$$

$${\Vert q(x)\Vert }_1$$

$$\le (x^{\top }{X}^{-1}x{)}^{1/2}\cdot (q^{\top }{X}^{-1}q{)}^{1/2}$$

$$=(1^{\top }x{)}^{1/2}\cdot (b^{\top }{L}^{-1}(x)b{)}^{1/2}$$

$$\le \frac{1}{2}{1}^{\top }x+\frac{1}{2}{b}^{\top }{L}^{-1}(x)b$$

$$=\frac{1}{2}f(x),$$

$$\frac{1}{2}f(x)\le \frac{ϵ}{2}+\frac{1}{2}\left(1+\frac{1}{1-ϵ}\right){\Vert s^{\ast }\Vert }_1.$$

$$1^{\top }x={\Vert s\Vert }_1\le {\Vert s^{\ast }\Vert }_1+ϵ.$$

Full text

1. Introduction

An important primitive in the areas of signal processing and optimization is that of finding a minimum $\ell_{1}$-norm solution to an underdetermined system of linear equations. Specifically, for some $n\leq m$, let $\hat{s}\in\mathbb{R}^{m}$ represent an unknown signal, $b\in\mathbb{R}^{n}$ a measurement vector, and $A\in\mathbb{R}^{n\times m}$ a full-rank matrix such that $A\hat{s}=b$. In some circumstances, the unknown signal $\hat{s}$ can be recovered by computing a minimum $\ell_{1}$-norm solution to the system $As=b$; in other words, solving the following optimization problem:

[TABLE]

This $\ell_{1}$-minimization problem is known as basis pursuit. It is a central problem in the theory of sparse representation and arises in several applications, such as imaging and face recognition. Through a standard reduction, it also captures the $\ell_{1}$-regression problem used in statistical estimation.

The convex optimization problem (BP) can be cast as a linear program and thus could be solved via an interior-point method. Another popular approach to $\ell_{1}$-minimization is the iteratively reweighted least squares (IRLS) method, which is based on iteratively solving a series of adaptively weighted $\ell_{2}$-minimization problems. IRLS methods are popular in practice, due to their simplicity and the fact that they do not require preprocessing nor special initialization rules. Despite this, theoretical guarantees for IRLS methods in the literature are not common, particularly in terms of global convergence bounds.

This work contributes to developing the understanding and design of IRLS-type methods for basis pursuit. We propose a novel exact reformulation of (BP) as a differentiable convex problem over the positive orthant, which we call the dissipation minimization problem. A distinguishing feature of this approach is that it entails the solution of a single differentiable convex problem. The reformulation leads naturally to a new family of IRLS-type methods solving (BP).

We exemplify this approach by providing global convergence bounds for discrete IRLS-type algorithms for (BP). We explore two possible routes to the solution of the dissipation minimization problem, and thus of (BP), where we use the established framework of first-order optimization methods to derive two provably convergent iterative algorithms. We bound their iteration complexity as $O(m^{2}/\epsilon^{3})$ and $O(m^{2}/\epsilon^{2})$, respectively, where $\epsilon$ is the relative error parameter. These methods are in the IRLS family since each iteration can be reduced to the solution of a weighted least squares problem. Both methods are very simple to implement and the first one exhibits a geometric convergence rate in numerical experiments.

Our dissipation-based reformulation of (BP) may be of independent interest. It is rooted in the Laplacian framework of network theory: it generalizes concepts such as the Laplacian matrix and the transfer matrix, which were originally developed to express the relation between electrical quantities across different terminals of a resistive network. (Many of our formulas have simple interpretations when the constraint matrix $A$ is derived from a network matrix).

This paper is organized as follows. In Section 2, we present the dissipation minimization reformulation of basis pursuit and some of its structural properties. In Section 3 we prove the equivalence between basis pursuit and dissipation minimization. In Section 4 we look at the continuous dynamics obtained by applying mirror descent to the dissipation minimization objective and connect them with existing literature. In Section 5, we analyze a discretization of these dynamics that yields an iterative IRLS-type method for the solution of the dissipation minimization problem and, hence, of basis pursuit; this method can be seen as an application of the well-known multiplicative weights update scheme, and its iteration complexity is $O(m^{2}/\epsilon^{3})$. Then, by leveraging Nesterov’s accelerated gradient scheme, we present and analyze an improved IRLS-type method with iteration complexity $O(m^{2}/\epsilon^{2})$. In Section 6, implementations of the two methods are compared against existing solvers from the l1benchmark suite [39].

Related literature. Given its central role in the areas of sparse representation and statistics, the literature on the basis pursuit problem and $\ell_{1}$-regression is extensive; see for example [11, 16, 18, 24] and references therein. Several algorithms for basis pursuit are reviewed in Chapter 15 of [24]; for an experimental comparison and an application to face recognition, see [39].

Various versions of IRLS schemes have been studied for a long time [26, 32] and, as already mentioned, the methods have been popular in practice due to their simplicity and experimental performance [17]. On the other hand, theoretical guarantees for IRLS-type algorithms are few and far between [7, 21, 35]. A recent IRLS algorithm stands out in the context of this paper, as it applies to the basis pursuit problem and comes with a worst-case guarantee: a $\tilde{O}(m^{1/3}\epsilon^{-8/3})$ iterations algorithm due to Chin et al. [19, Theorem 5.1], derived by further developing the approach of Christiano et al. [20]. In this context, our approach breaks the $\epsilon^{-8/3}$ bound for an IRLS method (at the cost of a worse dependency on $m$). We nevertheless emphasize that the goal of this work is not to establish the superiority of a specific algorithm, but rather to highlight a new approach that, already when coupled with off-the-shelf optimization methods, offers a principled way to derive IRLS-type algorithms with competitive theoretical performance. Subsequently to the first appearance of our results (on arXiv), an improved bound of $\tilde{O}(m^{1/3}\epsilon^{-2/3}+\epsilon^{-2})$ iterations for a more sophisticated IRLS-type algorithm for (BP) has been derived by Ene and Vladu [22] (again building on the ideas of [20] and [19]). While this algorithm has a rather more favorable worst-case dependency on the parameters, in practice it requires roughly $1/\epsilon$ iterations [22, Section 4]; in contrast, as we observe in Section 6, the experimental convergence rate of our approach is geometric, that is, the iterations required are linear in $\log(1/\epsilon)$, suggesting that a much stronger theoretical bound may hold in our setting.

Our reformulation of basis pursuit is new, though it is in part inspired by the Laplacian framework [34]. In particular, the definition of the dissipation function is based on a generalization of the Laplacian potential of a network. This reinforces the idea from Chin et al. [19] that concepts originally developed for network optimization can be fruitful in the context of $\ell_{1}$-regression. The dissipation-minimizing dynamics considered in Section 4 are an application of the mirror descent (or natural gradient) dynamics [30, 2, 3, 29] to our new objective function. In Section 5.1 we show, in particular, how the algorithmic framework of Lu, Freund and Nesterov [29] (see also [6]) can be applied to the dissipation minimization problem. The improved algorithm discussed in Section 5.2 is instead based on Nesterov’s well-known accelerated gradient method [31].

The dynamics studied in Sections 4 and 5 bear some formal similarity to the so-called Physarum dynamics, studied in the context of natural computing, which are the network dynamics of a slime mold [37, 14, 36, 35, 10]. The fact that Physarum dynamics are of IRLS type was first observed in [35]. In this context, our result can be seen as the derivation of a Physarum-like dynamics purely from an optimization principle: dissipation minimization following the natural gradient. A relevant difference is that the specific dynamics we study is a gradient system, while the dynamics studied in [35, 10] is provably not a gradient system. This is precisely what enables us to apply the machinery of first-order convex optimization methods, and acceleration in particular.

We note that a different proof of Theorem 3.1 has been independently provided by Facca, Cardin and Putti [23] in the context of the Physarum dynamics.

Notation. For a vector $x\in\mathbb{R}^{m}$, we use $\mathrm{diag}(x)$ to denote the $m\times m$ diagonal matrix with the coefficients of $x$ along the diagonal. The inner product of two vectors $x,y\in\mathbb{R}^{m}$ is denoted by $\langle x,y\rangle=x^{\top}y$. The maximum (respectively, minimum) eigenvalue of a diagonalizable matrix $M$ is denoted by $\lambda_{\mathrm{max}}(M)$ (respectively, $\lambda_{\mathrm{min}}(M)$). For a vector $x\in\mathbb{R}^{m}$, $\left\|x\right\|_{p}$ denotes the $\ell_{p}$-norm of $x$ ($1\leq p\leq\infty$), and $\left|x\right|$ denotes the vector $y$ such that $y_{i}=\left|x_{i}\right|$, $i=1,\ldots,m$. Similarly, $x^{2}$ denotes the vector $y$ such that $y_{i}=x_{i}^{2}$, $i=1,\ldots,m$. With a slight overlap of notation, which should nevertheless not cause any confusion, we instead reserve $x^{k}$ with a symbolic index $k$ to denote the vector produced by the $k$th step of an iterative algorithm.

2. Basis pursuit and the dissipation minimization problem

2.1. Assumptions on the basis pursuit problem

We make the following assumptions on (BP):

1. (A.1)

the matrix $A$ has full rank and $n\leq m$; 2. (A.2)

the system $As=b$ has at least one solution $s^{\prime}$ such that $s^{\prime}_{j}\neq 0$ for each $j=1,\ldots,m$.

Proposition 2.1**.**

Assumption (A.2) is without loss of generality, given (A.1).

Proof.

If the basis pursuit instance $(A,b)$ satisfies (A.1) but not (A.2), form a new instance $(A^{\prime},b)$ where $A^{\prime}$ is obtained from $A$ by duplicating every column. Observe the following about the two instances:

• •

$A^{\prime}$ has full rank and $n^{\prime}=n\leq m\leq 2m=m^{\prime}$.

• •

For any solution to $(A,b)$, there is a solution to $(A^{\prime},b)$ with the same cost.

• •

Let $u=A^{\top}(AA^{\top})^{-1}b$ be the least-square solution to $As=b$. There is at least one solution to $A^{\prime}s^{\prime}=b$ with $s^{\prime}_{j}\neq 0$ for each $j=1,\ldots,2m$, given by

[TABLE]

• •

No optimal solution to the instance $(A^{\prime},b)$ is such that $s^{\prime}_{2j-1}\cdot s^{\prime}_{2j}<0$ for some $j$: if that was the case, one could form a solution of lesser cost by replacing each of $s^{\prime}_{2j-1}$ and $s^{\prime}_{2j}$ with their average. Thus, any optimal solution $s^{\prime}$ to $(A^{\prime},b)$ can be transformed back into a solution $s$ to $(A,b)$ by taking $s_{j}=s^{\prime}_{2j-1}+s^{\prime}_{2j}$ for each $j=1,\ldots,m$. Such a solution satisfies $\left\|s\right\|_{1}=\left\|s^{\prime}\right\|_{1}$ and thus must be optimal for $(A,b)$.

∎

Remark 2.1*.*

A special case of (BP) is when $A$ is derived from a network matrix. Specifically, consider a connected network with $n+1$ nodes and $m$ edges, and suppose edge $j$ connects node $u$ to node $v$. Define $b_{j}\in\mathbb{R}^{m}$ as $(b_{j})_{u}=1$, $(b_{j})_{v}=-1$, and all other entries [math]. The matrix $B=[b_{1}\cdots b_{m}]\in\mathbb{R}^{(n+1)\times m}$ is called the incidence matrix of the network. For any connected network, the incidence matrix $B$ has rank $n$ and, additionally, any row of $B$ can be expressed as a linear combination of the remaining $n$ rows, because the sum of all rows is a zero vector. Let $A$ be the submatrix of $B$ obtained by deleting an arbitrary row. Then $A$ satisfies assumption (A.1) and thus, without loss of generality, (A.2). A solution $s$ to $As=b$ can be interpreted as an assignment of flow values to each edge such that the net in-flow at every node $v=1,\ldots,n$ matches the prescribed demand $b_{v}$.

2.2. The dissipation potential

In this section we introduce the dissipation potential, which is the function on which our reformulation of the basis pursuit problem is based.

Definition 2.1**.**

The Laplacian-like matrix relative to a vector $x\in\mathbb{R}^{m}_{\geq 0}$ is the matrix $L(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}AXA^{\top}$, where $X=\mathrm{diag}(x)$.

Remark 2.2*.*

In the network setting described in Remark 2.1, a vector $x\in\mathbb{R}^{m}_{>0}$ can be interpreted as a set of weights, or conductances, on the edges of the network. Then the matrix $BXB^{\top}$ is the weighted Laplacian of the network [12, 34]. The matrix $L(x)=AXA^{\top}$ is sometimes called the reduced Laplacian.

Proposition 2.2**.**

If $x>0$, then $L(x)$ is positive definite.

Proof.

Since $A$ has full rank, so has $AX^{1/2}$; hence $L(x)=(AX^{1/2})(AX^{1/2})^{\top}$ is positive definite. ∎

The following function definition is central to our approach.

Definition 2.2**.**

Let $A\in\mathbb{R}^{n\times m}$, $b\in\mathbb{R}^{n}$ be such that (A.1)–(A.2) hold. Define $f_{0},f:\mathbb{R}^{m}\to(-\infty,+\infty]$ as

[TABLE]

We call $f$ the dissipation potential. An equivalent definition of $f$ is as the convex closure of $f_{0}$, which is the function whose epigraph in $\mathbb{R}^{m+1}$ is the closure of the epigraph of $f_{0}$ [33, Chapter 7]. The effective domain of $f$ is the set

[TABLE]

The functions $f$ and $f_{0}$ differ only on the boundary of the positive orthant. We will show that $f$ always achieves a minimum on $\mathbb{R}^{m}_{\geq 0}$, and hence on $\mathbb{R}^{m}$. One of our main results (Theorem 3.1) is that this minimum equals the minimum of (BP).

Remark 2.3*.*

Consider again the case where the matrix $A$ is derived from a network matrix, as in Remark 2.1. The node of the network corresponding to the row that was removed from the incidence matrix to form $A$ is called the grounded node. Now assume that for some $u=1,\ldots,n$ the vector $b\in\mathbb{R}^{n}$ is such that $b_{v}=0$ if $v\neq u$, $b_{v}=1$ if $v=u$. Then the Laplacian potential $b^{\top}L^{-1}(x)b$ yields the effective resistance between the grounded node and node $u$ when the conductances of the network are specified by the vector $x$. A standard result in network theory is that decreasing the conductance of any edge can only increase the effective resistance between any two nodes (see, for example, [12, 25]). Thus, the minimization of the dissipation potential $f$ involves an equilibrium between two opposing tendencies: decreasing any $x_{j}$ decreases the linear term $\mathbf{1}^{\top}x$, but increases the Laplacian term $b^{\top}L^{-1}(x)b$.

2.3. Basic properties of the dissipation potential

We proceed to show that the dissipation potential attains a minimum. We start with some basic properties of $f_{0}$.

Lemma 2.3**.**

The function $f_{0}$ is positive, convex and differentiable on $\mathbb{R}^{m}_{>0}$.

Proof.

Positivity follows from the positive-definiteness of $L^{-1}(x)$ for $x\in\mathbb{R}^{m}_{>0}$ (implied by Proposition 2.2). For convexity, it suffices to show that the mapping $x\mapsto b^{\top}L^{-1}(x)b$ is convex on $\mathbb{R}^{m}_{>0}$. First observe that $x\mapsto AXA^{\top}$ is a linear matrix-valued function, i.e., each one of the entries of $AXA^{\top}$ is a linear function of $x$, since multiplying $X$ on the left and right with $A$ and $A^{\top}$ yields linear combinations of the elements of $x$. Second, the matrix to scalar function $Y\mapsto b^{\top}Y^{-1}b$ is convex on the cone of positive definite matrices, for any $b\in\mathbb{R}^{n}$ (see for example [15, Section 3.1.7]). By combining the two facts above, it follows that the composition $x\mapsto b^{\top}(AXA^{\top})^{-1}b$ is convex, and hence so is $f_{0}$. Finally, since the entries of $L(x)$ are linear functions of $x$, the function $f_{0}$ is a rational function with no poles in $\mathbb{R}^{m}_{>0}$, hence differentiable. ∎

To argue that $f$ attains a minimum, we first recall some notions from convex analysis [33, 8]. An extended real-valued function $f:\mathbb{R}^{m}\to[-\infty,+\infty]$ is called proper if its domain is nonempty and the function never attains the value $-\infty$. It is called closed if its epigraph is closed. It is called coercive if it is proper and $\lim_{\left\|x\right\|\to\infty}f(x)=+\infty$.

Lemma 2.4**.**

The function $f$ is nonnegative, proper, closed and convex on $\mathbb{R}^{m}$.

Proof.

By Lemma 2.3, $f_{0}$ is convex on $\mathbb{R}^{m}$, since it is convex on its effective domain. Moreover $f_{0}$ is proper, since $L^{-1}(x)$ is positive definite and thus $0<f_{0}(x)<+\infty$ for any $x\in\mathbb{R}^{m}_{>0}$. By construction, $f$ coincides with the closure of $f_{0}$ and thus it is a closed proper convex function [33, Theorem 7.4]. Its nonnegativity follows from the positivity of $f_{0}$ and from (2). ∎

Corollary 2.5**.**

The function $f$ attains a minimum on $\mathbb{R}^{m}_{\geq 0}$.

Proof.

Note that $\lim_{\left\|x\right\|\to\infty}f(x)=\infty$, because $b^{\top}(AXA^{\top})^{-1}b\geq 0$ for any $x\in\operatorname{dom}f_{0}$, and $\mathbf{1}^{\top}x\to\infty$ as $\left\|x\right\|\to\infty$ with $x\in\operatorname{dom}f_{0}$. In other words, $f$ is also a coercive function and therefore, it attains a minimal value over any nonempty closed set intersecting its domain [8, Theorem 2.14]; in particular, it attains its minimal value over $\mathbb{R}^{m}_{\geq 0}$. ∎

Since $f(x)=\liminf_{x^{\prime}\to x}f_{0}(x^{\prime})$, the minimum attained by $f$ over $\mathbb{R}^{m}_{\geq 0}$ equals $\inf_{x>0}f_{0}(x)$. Note also that this minimum may be attained on the boundary of $\operatorname{dom}f$.

2.4. Gradient and Hessian

In this section we derive some formulas for the gradient and Hessian of $f$ on the interior of its domain.

Definition 2.3**.**

Let $x\in\mathbb{R}^{m}_{>0}$. The voltage vector at $x$ is $d(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}A^{\top}L^{-1}(x)b\in\mathbb{R}^{m}$.

Remark 2.4*.*

In the network setting described in Remark 2.1, $d_{j}(x)$ expresses the voltage along edge $j$ when an external current $b_{u}$ enters each node $u=1,\ldots,n$ (and a balancing current $-\sum_{u}b_{u}$ enters the grounded node).

The next lemma relates the gradient $\nabla f(x)$ to the voltage vector at $x$.

Lemma 2.6**.**

Let $x\in\mathbb{R}^{m}_{>0}$. For any $j=1,\ldots,m$, $\frac{\partial f(x)}{\partial x_{j}}=1-(a_{j}^{\top}L^{-1}(x)b)^{2}=1-d^{2}_{j}(x),$ where $a_{j}$ stands for the $j$th column of $A$.

Proof.

First observe that $L(x)=AXA^{\top}=\sum_{j=1}^{m}x_{j}a_{j}a_{j}^{\top}$ and thus $\partial L/\partial x_{j}=a_{j}a_{j}^{\top}$. We apply the formula for the derivative of a matrix inverse:

[TABLE]

We obtain

[TABLE]

The claim follows by the definition of $f$. ∎

To express the Hessian of $f$, in addition to the voltages we need the notion of transfer matrix.

Definition 2.4**.**

Let $x\in\mathbb{R}^{m}_{>0}$. The transfer matrix at $x$ is $T(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}A^{\top}L^{-1}(x)A.$

Remark 2.5*.*

In the network setting described in Remark 2.1, the transfer matrix $T(x)$ expresses the relation between input currents and output voltages, when the conductances are given by the vector $x$. Namely, $T_{ij}(x)$ is the amount of voltage observed along edge $i$ of the network when a unit external current is applied between the endpoints of edge $j$.

Corollary 2.7**.**

For any $x>0$, $\nabla^{2}f(x)=2\cdot(d(x)\cdot d(x)^{\top})\odot T(x),$ where $\odot$ denotes the Schur matrix product defined by $(U\odot V)_{ij}=U_{ij}\cdot V_{ij}$.

Proof.

For any $i,j=1,\ldots,m$, by Lemma 2.6 and applying once more (3), we get

[TABLE]

The claim follows by Definition 2.4. ∎

2.5. Bounds on the norms of gradient and Hessian

In this section we derive some norm bounds for the gradient and Hessian of the dissipation potential $f$; they will be used crucially to derive complexity bounds for the algorithms studied in Section 5.

Two matrices $M$, $M^{\prime}$ are called congruent if there is a nonsingular matrix $S$ such that $M^{\prime}=SMS^{\top}$. For the proofs in this section, the main tool we rely on is the following algebraic fact relating the eigenvalues of congruent matrices; see for example [28, Theorem 4.5.9] for a proof.

Theorem 2.8** (Ostrowski).**

Let $M,S\in\mathbb{R}^{m\times m}$ be two symmetric matrices, with $S$ nonsingular. For $k=1,\ldots,m$, let $\lambda_{k}(M)$, $\lambda_{k}(SMS^{\top})$ denote the $k$-th largest eigenvalue of $M$ and $SMS^{\top}$, respectively. For each $k=1,\ldots,m$ there is a positive real number $\theta_{k}\in[\lambda_{\mathrm{min}}(SS^{\top}),\lambda_{\mathrm{max}}(SS^{\top})]$ such that

[TABLE]

Lemma 2.9**.**

Let $x\in\mathbb{R}^{m}_{>0}$. Each nonzero eigenvalue of $T(x)$ is at least $(\max_{i=1,\ldots,m}x_{i})^{-1}$ and at most $(\min_{i=1,\ldots,m}x_{i})^{-1}$.

Proof.

Consider the matrix $\Pi(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}X^{1/2}T(x)X^{1/2}$. By Definition 2.4,

[TABLE]

Hence, $\Pi(x)$ is the orthogonal projection matrix that projects onto the range of $(AX^{1/2})^{\top}$. In particular, $\Pi(x)^{2}=\Pi(x)$ and each eigenvalue of $\Pi(x)$ equals 0 or 1. Since $T(x)=X^{-1/2}\Pi(x)X^{-1/2}$, the matrices $T(x)$ and $\Pi(x)$ are congruent. By Theorem 2.8, the algebraic multiplicity of the zero eigenvalue of $T(x)$ and $\Pi(x)$ is the same, and each positive eigenvalue of $T(x)$ must lie between the smallest and the largest eigenvalue of $X^{-1}$. These are $(\max_{i}x_{i})^{-1}$ and $(\min_{i}x_{i})^{-1}$, respectively. ∎

Lemma 2.10**.**

Let $x\in\mathbb{R}^{m}_{>0}$. Then $\left\|d(x)\right\|_{\infty}\leq(\min_{i=1,\ldots,m}x_{i})^{-1}\cdot\left\|s\right\|_{2}$, where $s$ is any solution to $As=b$. In particular, for $c_{A,b}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}b^{\top}(AA^{\top})^{-1}b$,

[TABLE]

Additionally, if $s^{*}$ is an optimal solution to (BP),

[TABLE]

Proof.

Note that $d(x)=A^{\top}L^{-1}(x)b=A^{\top}L^{-1}As=T(x)s$. Hence

[TABLE]

Since the largest eigenvalue of $T(x)$ is at most $(\min_{i}x_{i})^{-1}$ by Lemma 2.9, we can bound $\left\|T(x)s\right\|_{2}\leq(\min_{i}x_{i})^{-1}\left\|s\right\|_{2}$, proving the first part of the claim. For the second part, consider the least square solution $u\stackrel{{\scriptstyle\mathrm{def}}}{{=}}A^{\top}(AA^{\top})^{-1}b$. Then $\left\|u\right\|_{2}=c_{A,b}^{1/2}$, and using the optimality of $u$ for the $\ell_{2}$ norm and of $s^{*}$ for the $\ell_{1}$ norm we derive

[TABLE]

Corollary 2.11**.**

If $x\in\mathbb{R}^{m}_{>0}$, then

[TABLE]

Proof.

Combine Lemma 2.10 with Lemma 2.6. ∎

Lemma 2.12**.**

If $x\in\mathbb{R}^{m}_{>0}$, then the largest eigenvalue of $\nabla^{2}f(x)$ satisfies

[TABLE]

Proof.

We can use the matrix identity $M\odot(zz^{\top})=\mathrm{diag}(z)\cdot M\cdot\mathrm{diag}(z)$ to reexpress Corollary 2.7 as

[TABLE]

where $D(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\mathrm{diag}(d(x))$. Hence, by Theorem 2.8, the largest eigenvalue of $\nabla^{2}f(x)$ satisfies

[TABLE]

for some $\theta$ lying between the smallest and largest eigenvalues of $D(x)^{2}$. Since by Lemma 2.10

[TABLE]

combining (10) and (11) with Lemma 2.9 we get $\lambda_{\mathrm{max}}(\nabla^{2}f(x))\leq 2(\min_{i}x_{i})^{-3}c_{A,b}$. ∎

3. Equivalence between basis pursuit and dissipation minimization

In this section we prove the equivalence between basis pursuit and dissipation minimization.

Theorem 3.1**.**

The value of the optimization problem

[TABLE]

is equal to the value of the optimization problem

[TABLE]

We call (DM) the dissipation minimization problem associated to $A$ and $b$. Note that the objective in (DM) is exactly $f_{0}(x)/2$, hence by (2) the minimum of (DM) equals the minimum of $f(x)/2$ over $\mathbb{R}^{m}_{\geq 0}$; the fact that this minimum is achieved is guaranteed by Corollary 2.5.

Definition 3.1**.**

Let $x>0$. The solution induced by $x$ is the vector $q(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}XA^{\top}L^{-1}(x)b$.

The term “solution” is justified by the fact that $Aq(x)=LL^{-1}b=b$. Induced solutions have the following simple characterization.

Lemma 3.2**.**

Let $x\in\mathbb{R}^{m}_{>0}$. The solution induced by $x$, $q(x)$, equals the unique optimal solution to the quadratic optimization problem:

[TABLE]

Proof.

This lemma is a straightforward generalization of Thomson’s principle [12, Chapter 9] from electrical network theory. We adapt an existing proof [13, Lemma 3] to the notation used in this paper. Since the objective function in (QPx) is strictly convex, the problem has a unique optimal solution. Consider any solution $s$, and let $r=s-q(x)$. Then $Ar=b-b=0$ and hence

[TABLE]

since $r^{\top}X^{-1}r\geq 0$ and $r^{\top}X^{-1}q=r^{\top}A^{\top}L^{-1}b=(Ar)^{\top}L^{-1}b=0$. Therefore, the objective function value of any solution $s$ to (QPx) is at least as large as the objective function value of the solution $q(x)$. ∎

The value of (QPx) is, in fact, the Laplacian potential $b^{\top}L^{-1}(x)b$.

Corollary 3.3**.**

The minimum of (QPx) equals $q(x)^{\top}X^{-1}q(x)=b^{\top}L^{-1}(x)b$.

Proof.

We already proved that the minimum of (QPx) is $q(x)^{\top}X^{-1}q(x)$. Substituting the definition of $q(x)$,

[TABLE]

Lemma 3.4**.**

For any $x>0$, $q(x)\in\mathbb{R}^{m}$ is such that $Aq=b$ and $\left\|q(x)\right\|_{1}\leq f(x)/2$. Thus, the value of (BP) is at most that of (DM).

Proof.

For any $x\in\mathbb{R}^{m}_{>0}$, consider its induced solution $q(x)=XA^{\top}L(x)^{-1}b$. We already observed that $q(x)$ is feasible for (BP). Moreover, we can bound:

[TABLE]

where the first upper bound follows from the Cauchy-Schwarz inequality, and the second from the Arithmetic Mean-Geometric Mean inequality. ∎

To prove the converse of Lemma 3.4, we develop an intermediate lemma that relates the value of an optimal solution $s^{*}$ of (BP) to the dissipation value of a vector $x$ such that $x=\left|s\right|$ with $s$ sufficiently close to $s^{*}$.

Lemma 3.5**.**

Let $s\in\mathbb{R}^{m}$, $\epsilon\in(0,1)$ be such that $As=b$, $s_{j}\neq 0$ and $(1-\epsilon)\left|s^{*}_{j}\right|\leq\left|s_{j}\right|\leq\left|s^{*}_{j}\right|+\epsilon/m$ for some $s^{*}$ such that $As^{*}=b$ and each $j=1,\ldots,m$. Then for $x=|s|$,

[TABLE]

Proof.

On one hand, by the assumed upper bound $\left|s_{j}\right|\leq|s^{*}_{j}|+\epsilon/m$, trivially

[TABLE]

On the other hand, consider the solution $q(x)$ induced by $x$ and recall that $q(x)$ is feasible for (BP), since $Aq=b$, and optimal for (QPx). By the assumed lower bound $\left|s_{j}\right|\geq(1-\epsilon)\left|s^{*}_{j}\right|$, and by Lemma 3.2,

[TABLE]

where the first upper bound follows from the fact that $s^{*}$ is a solution to (QPx), and the second follows from the hypothesis. Combining (13) and (14), we get

[TABLE]

Lemma 3.6**.**

The value of (DM) is at most that of (BP).

Proof.

Consider an optimal solution $s^{*}\in\mathbb{R}^{m}$ to (BP). Let $s^{\prime}\in\mathbb{R}^{m}$ be a solution to $As=b$ such that $s^{\prime}_{j}\neq 0$ for all $j=1,\ldots,m$ (such an $s^{\prime}$ exists by assumption (A.2)). For any $\delta\in(0,1)$, let $s(\delta)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(1-\delta)s^{*}+\delta s^{\prime}$ and $x(\delta)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left|s(\delta)\right|>0$. For any $\epsilon\in(0,1)$ we can ensure that the hypotheses of Lemma 3.5 are satisfied by choosing a small enough $\delta>0$. For such a value of $\delta$, Lemma 3.5 yields

[TABLE]

As $\epsilon$ can be chosen arbitrarily small, and the right-hand side of (15) approaches $\left\|s^{*}\right\|_{1}$ as $\epsilon\to 0$, we obtain the claim. ∎

This concludes the proof of Theorem 3.1. Not only are the optimal values of (BP) and (DM) the same, but one can bound the suboptimality of any feasible point of (BP) in terms of the dissipation value of a corresponding vector.

Theorem 3.7**.**

Let $s\in\mathbb{R}^{m}$ be a feasible point of (BP) such that $s_{j}\neq 0$ for all $j=1,\ldots,m$, and let $x=\left|s\right|$, $\rho(x)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left\|d(x)\right\|_{\infty}$. The quantity $\left(1+\rho^{-1}(x)\right)\left\|s\right\|_{1}-\rho^{-1}(x)\cdot f(x)$ is an upper bound on the suboptimality of $s$.

Proof.

Consider the following linear formulation of (BP) (left) and its dual (right):

[TABLE]

[TABLE]

Given any solution $s$ to (BP) such that $x=\left|s\right|>0$, let us take

[TABLE]

Then $\left\|A^{\top}\nu\right\|_{\infty}\leq 1$ by definition of $\rho(x)$; moreover, $\lambda+\mu=\mathbf{1}$, $\lambda-\mu+A^{\top}\nu=0$, and $\lambda,\mu\geq 0$. Thus, $(x,s)$ is a primal feasible solution, $(\lambda,\mu,\nu)$ is a dual feasible solution, and by weak duality

[TABLE]

This implies a duality gap of

[TABLE]

We close this section by observing that a simpler proof of Theorem 3.1 can be obtained by the following quadratic variational formulation of the $\ell_{1}$-norm: for any $s\in\mathbb{R}^{m}$,

[TABLE]

see, for example, Bach et al. [4, Section 1.4.2]. Therefore

[TABLE]

where the last identity follows from Corollary 3.3. However, the full strength of Lemma 3.4 and Lemma 3.6 is crucial to be able to constructively transform feasible points for (DM) into feasible points for (BP) and vice versa.

4. Continuous dynamics for dissipation minimization

Theorem 3.1 readily suggests an approach to the solution of the basis pursuit problem. Namely, the solution of the non-smooth, equality constrained formulation (BP) is reduced to the solution of the differentiable formulation (DM) on the positive orthant.

Mirror descent dynamics. To solve (DM), it is natural to adopt methods for differentiable constrained optimization that are designed for simple constraints. Consider first the following set of ordinary differential equations, aimed at solving $\inf\,\{f(x)\,|\,x>0\}$:

[TABLE]

with initial condition $x(0)=x^{0}$ for some $x^{0}>0$. When $f$ is the dissipation potential, by Lemma 2.6 this yields the explicit dynamics

[TABLE]

The dynamical system (16) is a nonlinear Lotka-Volterra type system of differential equations, of a kind that is common in population dynamics [27]. It is also an example of a Hessian gradient flow [1]: it can be expressed in the form

[TABLE]

where $H(x)=\nabla^{2}h(x)$ is the Hessian of a convex function $h$; namely, here $H(x)=X^{-1}$, and $h:\mathbb{R}^{m}_{>0}\to\mathbb{R}$ is the negative entropy function

[TABLE]

System (18) can also be expressed as $\frac{d}{dt}\frac{\partial h(x)}{\partial x_{j}}=-\frac{\partial f(x)}{\partial x_{j}},j=1,\ldots,m,$ or more succinctly,

[TABLE]

which is known as the mirror descent dynamics or natural gradient flow [30, 2]. The well-posedness of (18) has been considered, for example, in [1]. A dynamics formally similar to (17) is the Physarum dynamics [14, 36, 35, 10], namely,

[TABLE]

Differently from (17), the dynamics (21) is not a gradient flow, that is, there is no function $f$ that allows to write the dynamics in the form (18) or (20) (with $h$ the negative entropy).

Convergence of the dynamics. The fact that the solution of the mirror descent dynamics (18) converges to a minimizer of $f$ with rate $1/t$ is a well-known result; see, for example, [1, 38]. We include a streamlined proof for completeness.

Lemma 4.1**.**

The values $f(x(t))$ with $x(t)$ given by (16) are nonincreasing in $t$.

Proof.

We compute

[TABLE]

A key role in the convergence of the mirror descent dynamics is played by the Bregman divergence of the function $h$.

Definition 4.1**.**

The Bregman divergence of a convex function $h:\mathbb{R}^{m}\to(-\infty,+\infty]$ is defined by $D_{h}(x,y)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}h(x)-h(y)-\langle\nabla h(y),x-y\rangle.$

Convexity of $h$ implies the nonnegativity of $D_{h}(x,y)$. When $h$ is the negative entropy, $D_{h}$ is the relative entropy function (also known as Kullback-Leibler divergence), for which $D_{h}(x,y)=0$ if and only if $x=y$.

Theorem 4.2** ([1, 38]).**

*Let $x^{*}\in\mathbb{R}^{m}_{\geq 0}$ be a minimizer of $f$. As $t\to\infty$, the values $f(x(t))$ with $x(t)$ given by (16) converge to $f(x^{*})$. In particular, *

[TABLE]

Proof.

In the following, to shorten notation we often write $x$ in place of $x(t)$. Since $(d/dt)\nabla h(x)+\nabla f(x)=0$ by (18), for any $y$ we have $\langle(d/dt)\nabla h(x)+\nabla f(x),x-y\rangle=0$. This is equivalent to

[TABLE]

On the other hand, since $(d/dt)h(x)=\langle\nabla h(x),\dot{x}\rangle$, a simple calculation shows

[TABLE]

Combining (22) and (23), and plugging in $y=x^{*}$,

[TABLE]

The proof is concluded by a potential function argument [5, 38]. Consider the function

[TABLE]

Its time derivative is, by (24),

[TABLE]

where the last summand is nonpositive by Lemma 4.1 and the other terms equal, by definition, $-D_{f}(x^{*},x)\leq 0$. Hence, $\mathcal{E}(t)\leq\mathcal{E}(0)$ for all $t\geq 0$, which is equivalent to

[TABLE]

proving the claim. ∎

5. Algorithms for dissipation minimization

We now turn to the problem of designing IRLS-type algorithms for (DM) (and thus (BP)) with provably bounded iteration complexity. Two technical obstacles in the setup of a first-order method for formulation (DM) are: 1) that the positive orthant is not a closed set, and 2) that the gradients of $f$ may not be uniformly bounded on the positive orthant. There is a way to deal with both issues at once: instead of solving $\inf_{x>0}f(x)$, for an appropriately small $\delta>0$ one can minimize $f$ over

[TABLE]

This is established by the next lemma.

Lemma 5.1**.**

Let $x^{*}$ be a minimizer of $f$. Then $f(x^{*})\leq\min_{x\in\Omega_{\delta}}f(x)\leq f(x^{*})+\delta\,m$.

Proof.

The first inequality is trivial. As for the second, recall that $f(x)=\mathbf{1}^{\top}x+b^{\top}L^{-1}(x)b$ for any $x>0$, and that in the latter sum, the second term is non-increasing with $x$ (by Lemma 2.6). Thus, for any $x>0$,

[TABLE]

In other words, for any $x>0$, there is $y\geq\delta\mathbf{1}$ (namely, $y=x+\delta\mathbf{1}$) such that $f(y)\leq f(x)+\delta m$. ∎

In the following, we let $\delta\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\epsilon\,c_{A,b}^{1/2}/(2m)$, where $\epsilon$ is the desired error factor and $c_{A,b}$ is as defined in Lemma 2.10; this, by Lemma 2.10 and Theorem 3.1, ensures that the additional error incurred by restricting solutions to $\Omega_{\delta}$ is at most $(\epsilon/2)\left\|s^{*}\right\|_{1}=(\epsilon/4)f(x^{*})$.

5.1. Primal gradient scheme

Guided by (20), we might consider its forward Euler discretization

[TABLE]

where $x^{k}\in\Omega_{\delta}$ denotes the $k$th iterate, and $\eta\in\mathbb{R}_{>0}$ an appropriate step size. Indeed, the update (25) falls within a well-studied methodology for first-order convex optimization [9, 29]. We adapt this framework to the solution of (DM).

The primal gradient scheme is a first-order method for minimizing a differentiable convex function $f$ over a closed convex set $Q$. This scheme, which is defined with respect to a reference function $h$, proceeds as follows [29, 6]:

1. (1)

Initialize $x^{0}\in Q$. Let $\beta>0$ be a parameter. 2. (2)

At iteration $k=0,1,\ldots$, compute $\nabla f(x^{k})$ and set

[TABLE]

We apply the scheme with $h$ as defined in (19) and with $Q=\Omega_{\delta}$. Then, the minimization in (26) can be carried out analytically; it reduces to

[TABLE]

Update (27) is straightforward to implement as long as one can compute $\nabla f(x^{k})$. This computation is discussed in Section 5.3.

Convergence of the primal gradient scheme. As shown in [29], the primal gradient scheme achieves an absolute error bounded by $O(\beta/k)$ after $k$ iterations provided that the function $f$ is $\beta$-smooth relative to $h$. In our case, where both $f$ and $h$ are twice-differentiable on $Q$, relative $\beta$-smoothness is defined as

[TABLE]

Theorem 5.2** ([29]).**

If $f$ is $\beta$-smooth relative to $h$, then for all $k\geq 1$, the updates (26) satisfy

[TABLE]

where $x^{*}|_{Q}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\operatorname*{argmin}_{x\in Q}f(x)$.

To apply Theorem 5.2 in our setting, we need to bound the smoothness parameter $\beta$. We do this by leveraging the bounds derived in Section 2.5.

Lemma 5.3**.**

Equation (28) holds for $\beta=8m^{2}/\epsilon^{2}$.

Proof.

Condition (28) is equivalent to the condition that the largest eigenvalue of the matrix

[TABLE]

be at most $\beta$ (see [28, Theorem 7.7.3]). The matrix $X\nabla^{2}f(x)$ is similar to $X^{1/2}\nabla^{2}f(x)X^{1/2}$, hence it suffices to bound the eigenvalues of the latter. Since $\nabla^{2}f(x)=2D(x)T(x)D(x)$ with $D(x)=\mathrm{diag}(d(x))$,

[TABLE]

where we used the fact that $X$ and $D(x)$ are diagonal. By the proof of Lemma 2.9, the eigenvalues of $\Pi(x)$ are all 0 or 1. Hence, using again the relation between the eigenvalues of congruent matrices (Theorem 2.8), we conclude that the largest eigenvalue of $X^{1/2}\nabla^{2}f(x)X^{1/2}$ is bounded by that of $2D(x)^{2}$. Since $D(x)=\mathrm{diag}(d(x))$, the latter equals $2\left\|d(x)\right\|_{\infty}^{2}$, which is $2c_{A,b}/\delta^{2}=8m^{2}/\epsilon^{2}$ by Lemma 2.10 and the definitions of $\Omega_{\delta}$ and $\delta$. ∎

Theorem 5.4**.**

The primal gradient scheme (27) applied to the dissipation minimization problem (DM) achieves relative error at most $\epsilon$ after $96m^{2}\log(m/\epsilon)/\epsilon^{3}=\tilde{O}(m^{2}/\epsilon^{3})$ iterations.

Proof.

By Theorem 5.2 and Lemma 5.3, after $k$ iterations it holds that

[TABLE]

where $R\stackrel{{\scriptstyle\mathrm{def}}}{{=}}D_{h}(x^{*}|_{Q},x^{0})$. Since $f(x^{*}|_{Q})\leq(1+\epsilon/4)f(x^{*})$ (by Lemma 5.1, since $Q=\Omega_{\delta}$), this implies

[TABLE]

Thus, $f(x^{k})-f(x^{*})\leq\epsilon f(x^{*})$ if we take $k=\lceil 32Rm^{2}/(3\epsilon^{3}f(x^{*}))\rceil$. We complete the proof by bounding $R/f(x^{*})$ in terms of $\log(m/\epsilon)$. Let

[TABLE]

Observe that since $x^{0}\in Q$,

[TABLE]

with the last inequality following from (6). Thus,

[TABLE]

Hence, $k=\lceil 96m^{2}\log(m/\epsilon)/\epsilon^{3}\rceil$ iterations suffice to achieve relative error $\epsilon$. ∎

5.2. Accelerated gradient scheme

The second optimization scheme that we consider is the accelerated gradient method of Nesterov [31]. This can be summarized as follows:

1. (1)

Initialize $x^{0}\in Q$. Let $\beta>0$ be a parameter. 2. (2)

At iteration $k=0,1,\ldots$, compute $\nabla f(x^{k})$ and set $\alpha_{k}=1/2(k+1)$, $\tau_{k}=2/(k+3)$ and

[TABLE]

In our application of the scheme, $Q=\Omega_{\delta}$ and the minimization in (31) and (32) can be carried out analytically; explicitly, they become

[TABLE]

To implement (34)–(35), it is enough to be able to access the gradient $\nabla f(x^{k})$ and the cumulative gradient $\sum_{i}\alpha_{i}\nabla f(x^{i})$; the latter can be maintained with one additional update at each iteration.

Convergence of the accelerated gradient scheme. The well-known result by Nesterov [31] shows that the accelerated gradient scheme achieves an absolute error bounded by $O(\beta/k^{2})$ after $k$ iterations provided that the gradient of the function $f$ is $\beta$-Lipschitz-continuous over $Q$. In our case, where $f$ is twice-differentiable on $Q$, this means

[TABLE]

Theorem 5.5** ([31]).**

If $\nabla f$ is $\beta$-Lipschitz-continuous over $Q$, then for all $k\geq 1$, the updates (31)–(33) satisfy

[TABLE]

where $x^{*}|_{Q}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\operatorname*{argmin}_{x\in Q}f(x)$.

Again, to apply Theorem 5.5 in our setting, we need to bound the smoothness parameter $\beta$. We do this by exploiting Lemma 2.12.

Lemma 5.6**.**

Equation (36) holds for $\beta=16m^{3}/(\epsilon^{3}c_{A,b}^{1/2})$.

Proof.

Immediate from Lemma 2.12, the fact that $Q=\Omega_{\delta}$ and the definition of $\Omega_{\delta}$. Recall that $\delta=\epsilon c_{A,b}^{1/2}/(2m)$. ∎

Theorem 5.7**.**

If $x^{0}=\left|u\right|$ where $u\stackrel{{\scriptstyle\mathrm{def}}}{{=}}A^{\top}(AA^{\top})^{-1}b$ is the least square solution to $As=b$, the accelerated gradient scheme (31)–(33) applied to the dissipation minimization problem (DM) achieves relative error at most $\epsilon$ after $24m^{2}/\epsilon^{2}$ iterations.

Proof.

By Theorem 5.5 and Lemma 5.6, after $k$ iterations it holds that

[TABLE]

Since $f(x^{*}|_{Q})\leq(1+\epsilon/4)f(x^{*})$ by Lemma 5.1, this implies

[TABLE]

Thus, $f(y^{k})-f(x^{*})\leq(\epsilon/4)f(x^{*})+(\epsilon/2)c_{A,b}^{1/2}<\epsilon f(x^{*})$ if the number of iterations $k$ is at least

[TABLE]

We complete the proof by bounding $(R/c_{A,b})^{1/2}=\left\|x^{*}|_{Q}-x^{0}\right\|_{2}/c_{A,b}^{1/2}$ in terms of $m$. Observe that $R^{1/2}=\left\|x^{*}|_{Q}-x^{0}\right\|_{2}\leq\left\|x^{*}|_{Q}\right\|_{2}+\left\|x^{0}\right\|_{2}$. By the assumption that $x^{0}=|u|$ where $u$ is the least square solution to $As=b$, $\left\|x^{0}\right\|_{2}=\left\|u\right\|_{2}=c_{A,b}^{1/2}$ (recall the definition of $c_{A,b}$ in Lemma 2.10). Moreover,

[TABLE]

Hence $(R/c_{A,b})^{1/2}\leq 3m^{1/2}$ and substitution in (39) yields the theorem. ∎

5.3. Implementing the iterations

We conclude this section by commenting on a few implementations details and in particular on how each iteration of (27) and (31)–(33) could be implemented. A notable point is that each iteration can be reduced to a series of operations that access the matrix $A$ only through the solution of a system of the form $AWA^{\top}p=b$, for some diagonal matrix $W$, or through matrix-vector multiplications of the form $Ax$ or $A^{\top}x$.

Computation of the gradient. By Lemma 2.6, computing the vector $d(x)=A^{\top}L^{-1}(x)b$ is enough to compute the gradient at $x$, since $\nabla f(x)=\mathbf{1}-d^{2}(x)$. To compute $d(x)$, it is enough to solve the linear system $L(x)p=b$ for $p$, then premultiply the solution with $A^{\top}$. Note that since $L(x)=AXA^{\top}$, the system $L(x)p=b$ is a symmetric linear system with a positive definite constraint matrix.

Warm start. Heuristically, the solution of the system $L(x^{k+1})p=b$, which is required to compute the gradient at iteration $k+1$, can be expected to be close to that of the system $L(x^{k})p=b$ when $x^{k+1}$ is close to $x^{k}$. Hence, one possibility in practice is to use the solution obtained at step $k$ to warm-start the linear equation solver at step $k+1$, with a possible substantial reduction in the computational cost of each iteration.

Initial point and exit criterion. We assumed the starting point is the least square solution in Theorem 5.7, but this was only to optimize the worst-case iteration bound. In fact, Theorem 5.4 and Eq. (39) always apply and the schemes we discussed do not require a special initialization apart from membership into $\Omega_{\delta}$; hence, any point that is not too close to the boundary of the positive orthant is a suitable starting point. We can stop the schemes after the number of iterations $k$ is large enough to ensure the error guarantees of Theorems 5.4 and 5.7 (or Eq. (39)). Alternatively, a natural exit criterion in practice can be based on the duality gap provided by Theorem 3.7.

Obtaining feasible iterates for (BP). The algorithms as described above produce iterates in the positive orthant, that is, iterates that are feasible for (DM), but after all, our goal was to obtain feasible iterates of (BP). By using the ideas of Lemma 3.4, we can easily associate with any iterate $x^{k}\in\mathbb{R}^{m}_{>0}$ an iterate $s^{k}$ that is feasible for (BP), and the cost of which is not larger than the dissipation cost of $x^{k}$: namely, take $s^{k}=q(x^{k})=X^{k}A^{\top}L(x^{k})^{-1}b$. By the proof of Lemma 3.4, we know that $\left\|s^{k}\right\|_{1}\leq f(x^{k})/2$. Thus, the error bounds for $f(x^{k})$ can be directly translated into error bounds for $\left\|s^{k}\right\|_{1}$. Note that $s^{k}$ can be computed essentially for free, since $s^{k}=X^{k}d(x^{k})$ and $d(x^{k})$ is a byproduct of the gradient computation at iteration $k$.

6. Numerical comparison with other algorithms for $\ell_{1}$-minimization

We include in this section a numerical comparison of our schemes to other well-known algorithms for $\ell_{1}$-minimization. The results suggest that both the primal scheme and a slightly revised accelerated scheme converge at a geometric rate, that is, much faster than what our theoretical analysis guarantees. This suggests the open problem of improving the quality of our error bounds.

To compare our approaches to other algorithms for $\ell_{1}$-minimization, we implemented them in MATLAB and ran the l1benchmark suite by Yang et al. [39], which includes implementations of many other $\ell_{1}$-minimization solvers. A representative comparison is shown in Figure 1.

The figure plots the relative error of the algorithms as a function of computation time, averaged on 20 randomly generated instances (with $m=1000$, $n=800$, and 20% or 30% nonzeros in the ground truth solution). The implementations based on our approaches are:

• •

the Primal Gradient Scheme of Section 5.1 (PGS, with $\beta=3.5$, $\delta=10^{-15}$),

• •

the Accelerated Gradient Scheme of Section 5.2 (AGS, with $\beta=3.5$, $\delta=10^{-15}$, $\tau_{k}=2/(k+3)$), and

• •

a revised Accelerated Gradient Scheme, which we formulate below (AGS2, with $\beta=1.1$, $\delta=10^{-15}$, $\tau_{k}=10^{-15}$).

Other algorithms measured in the experiment are the Homotopy method, the primal and dual augmented Lagrangian methods (PALM, DALM), the primal-dual interior point method (PDIPA), the truncated Newton interior point method (L1LS), the fast iterative soft-thresholding method (FISTA), and the approximate message passing method (AMP). We refer the reader to Yang et al. [39] for references and discussion of these other methods. We observe, incidentally, that many of these methods construct points that are only approximately feasible for (BP), since they relax the constraint $As=b$ into the objective function, in one form or the other.

In the experiments, PGS clearly exhibits a geometric convergence rate, which is much better than what Theorem 5.4 guarantees, strongly suggesting that an improved theoretical analysis may be possible. Over time, PGS essentially reaches the machine precision barrier ($\approx 10^{-15}$), which is not true for other methods in the benchmark, such as FISTA, L1LS or the interior point method (PDIPA).

AGS, on the other hand, appears to be rather inaccurate in practice and does not exhibit a substantially better behavior than what is guaranteed by Theorem 5.7. This suggests that the entropic form of the updates – used in PGS but not in AGS – might have a high impact in practice. Therefore, we also benchmark a revised algorithm (AGS2) obtained by adopting an entropic form of the AGS updates (34)–(35), as follows (colored terms are new):

[TABLE]

The resulting scheme AGS2 is seen in Figure 1 to exhibit a geometric convergence rate and to be competitive against some of the best results in the benchmark, such as those of the primal augmented Lagrangian method (PALM).

7. Conclusions

We proposed a novel exact reformulation of the basis pursuit problem, which leads to a new family of gradient-based, IRLS-type methods for its solution. We then analyzed the iteration complexity of a natural optimization approach to the reformulation, based on the mirror descent scheme, as well as the iteration complexity of an accelerated gradient method. The first scheme can be seen as the discretization of a Hessian gradient flow and also as a variant on the Physarum dynamics, derived purely from optimization principles. The accelerated method, on the other hand, improves the error dependency for IRLS-type methods for basis pursuit, from $\epsilon^{-8/3}$ to $\epsilon^{-2}$. The experimental convergence rate of the first scheme, as well as that of a simple variant the second scheme, appears to be geometric. We interpret this as evidence that the dissipation minimization perspective may stimulate even more approaches to the design and analysis of efficient and practical IRLS-type methods.