algorithms | regularize

June 6, 2018

Best stepsizes for Douglas-Rachford for Lipschitz + strongly montone

Posted by Dirk under Math, Optimization, Uncategorized | Tags: algorithms, Douglas-Rachford, monotone inclusion, stepsize |
Leave a Comment

This is a short follow up on my last post where I wrote about the sweet spot of the stepsize of the Douglas-Rachford iteration. For the case $\beta$ -Lipschitz + $\mu$ -strongly monotone, the iteration with stepsize $t$ converges linear with rate

$\displaystyle r(t) = \tfrac{1}{2(1+t\mu)}\left(\sqrt{2t^{2}\mu^{2}+2t\mu + 1 +2(1 - \tfrac{1}{(1+t\beta)^{2}} - \tfrac1{1+t^{2}\beta^{2}})t\mu(1+t\mu)} + 1\right)$

Here is animated plot of this contraction factor depending on $\beta$ and $\mu$ and $t$ acts as time variable:

DR_contraction

What is interesting is, that this factor has increasing or decreasing in $t$ depending on the values of $\beta$ and $\mu$ .

For each pair $(\beta,\mu)$ there is a best $t^*$ and also a smallest contraction factor $r(t^*)$ . Here are plots of these quantities:
DR_opt_stepsize

Comparing the plot of te optimal contraction factor to the animated plot above, you see that the right choice of the stepsize matters a lot.

June 6, 2018

The sweet spot of the Douglas-Rachford iteration

Posted by Dirk under Math, Optimization | Tags: algorithms, Douglas-Rachford, monotone inclusion, stepsize |
1 Comment

I blogged about the Douglas-Rachford method before here and here. It’s a method to solve monotone inclusions in the form

$\displaystyle 0 \in Ax + Bx$

with monotone multivalued operators ${A,B}$ from a Hilbert space into itself. Using the resolvent ${J_{A} = (I+A)^{-1}}$ and the reflector ${R_{A} = 2J_{A} - I}$ , the Douglas-Rachford iteration is concisely written as

$\displaystyle u^{n+1} = \tfrac12(I + R_{B}R_{A})u_{n}.$

The convergence of the method has been clarified is a number of papers, see, e.g.

Lions, Pierre-Louis, and Bertrand Mercier. “Splitting algorithms for the sum of two nonlinear operators.” SIAM Journal on Numerical Analysis 16.6 (1979): 964-979.

for the first treatment in the context of monotone operators and

Svaiter, Benar Fux. “On weak convergence of the Douglas–Rachford method.” SIAM Journal on Control and Optimization 49.1 (2011): 280-287.

for a recent very general convergence result.

Since ${tA}$ is monotone if ${A}$ is monotone and ${t>0}$ , we can introduce a stepsize for the Douglas-Rachford iteration

$\displaystyle u^{n+1} = \tfrac12(I + R_{tB}R_{tA})u^{n}.$

It turns out, that this stepsize matters a lot in practice; too small and too large stepsizes lead to slow convergence. It is a kind of folk wisdom, that there is “sweet spot” for the stepsize. In a recent preprint Quoc Tran-Dinh and I investigated this sweet spot in the simple case of linear operators ${A}$ and ${B}$ and this tweet has a visualization.

A few days ago Walaa Moursi and Lieven Vandenberghe published the preprint “Douglas-Rachford splitting for a Lipschitz continuous and a strongly monotone operator” and derived some linear convergence rates in the special case they mention in the title. One result (Theorem 4.3) goes as follows: If ${A}$ is monotone and Lipschitz continuous with constant ${\beta}$ and ${B}$ is maximally monotone and ${\mu}$ -strongly monotone, than the Douglas-Rachford iterates converge strongly to a solution with a linear rate

$\displaystyle r = \tfrac{1}{2(1+\mu)}\left(\sqrt{2\mu^{2}+2\mu + 1 +2(1 - \tfrac{1}{(1+\beta)^{2}} - \tfrac1{1+\beta^{2}})\mu(1+\mu)} + 1\right).$

This is a surprisingly complicated expression, but there is a nice thing about it: It allows to optimize for the stepsize! The rate depends on the stepsize as

$\displaystyle r(t) = \tfrac{1}{2(1+t\mu)}\left(\sqrt{2t^{2}\mu^{2}+2t\mu + 1 +2(1 - \tfrac{1}{(1+t\beta)^{2}} - \tfrac1{1+t^{2}\beta^{2}})t\mu(1+t\mu)} + 1\right)$

and the two plots of this function below show the sweet spot clearly.

If one knows the Lipschitz constant of ${A}$ and the constant of strong monotonicity of ${B}$ , one can minimize ${r(t)}$ to get on optimal stepsize (in the sense that the guaranteed contraction factor is as small as possible). As Moursi and Vandenberghe explain in their Remark 5.4, this optimization involves finding the root of a polynomial of degree 5, so it is possible but cumbersome.

Now I wonder if there is any hope to show that the adaptive stepsize Quoc and I proposed here (which basically amounts to ${t_{n} = \|u^{n}\|/\|Au^{n}\|}$ in the case of single valued ${A}$ – note that the role of ${A}$ and ${B}$ is swapped in our paper) is able to find the sweet spot (experimentally it does).
<p

February 24, 2017

The origin of the Douglas-Rachford iteration

Posted by Dirk under Math, Optimization | Tags: algorithms, Douglas-Rachford, monotone inclusion, optimization, publishing |
1 Comment

I blogged about the Douglas-Rachford method before and in this post I’d like to dig a bit into the history of the method.

As the name suggests, the method has its roots in a paper by Douglas and Rachford and the paper is

Douglas, Jim, Jr., and Henry H. Rachford Jr., “On the numerical solution of heat conduction problems in two and three space variables.” Transactions of the American mathematical Society 82.2 (1956): 421-439.

At first glance, the title does not suggest that the paper may be related to monotone inclusions and if you read the paper you’ll not find any monotone operator mentioned. So let’s start and look at Douglas and Rachford’s paper.

1. Solving the heat equation numerically

So let us see, what they were after and how this is related to what is known as Douglas-Rachford splitting method today.

Indeed, Douglas and Rachford wanted to solve the instationary heat equation

$\displaystyle \begin{array}{rcl} \partial_{t}u &=& \partial_{xx}u + \partial_{yy}u \\ u(x,y,0) &=& f(x,y) \end{array}$

with Dirichlet boundary conditions (they also considered three dimensions, but let us skip that here). They considered a rectangular grid and a very simple finite difference approximation of the second derivatives, i.e.

$\displaystyle \begin{array}{rcl} \partial_{xx}u(x,y,t)&\approx& (u^{n}_{i+1,j}-2u^{n}_{i,j}+u^{n}_{i-1,j})/h^{2}\\ \partial_{yy}u(x,y,t)&\approx& (u^{n}_{i,j+1}-2u^{n}_{i,j}+u^{n}_{i,j-1})/h^{2} \end{array}$

(with modifications at the boundary to accomodate the boundary conditions). To ease notation, we abbreviate the difference quotients as operators (actually, also matrices) that act for a fixed time step

$\displaystyle \begin{array}{rcl} (Au^{n})_{i,j} &=& (u^{n}_{i+1,j}-2u^{n}_{i,j}+u^{n}_{i-1,j})/h^{2}\\ (Bu^{n})_{i,j} &=& (u^{n}_{i,j+1}-2u^{n}_{i,j}+u^{n}_{i,j+1})/h^{2}. \end{array}$

With this notation, our problem is to solve

$\displaystyle \begin{array}{rcl} \partial_{t}u &=& (A+B)u \end{array}$

in time.

Then they give the following iteration:

$\displaystyle Av^{n+1}+Bw^{n} = \frac{v^{n+1}-w^{n}}{\tau} \ \ \ \ \ (1)$

$\displaystyle Bw^{n+1} = Bw^{n} + \frac{w^{n+1}-v^{n+1}}{\tau} \ \ \ \ \ (2)$

(plus boundary conditions which I’d like to swipe under the rug here). If we eliminate ${v^{n+1}}$ from the first equation using the second we get

$\displaystyle (A+B)w^{n+1} = \frac{w^{n+1}-w^{n}}{\tau} + \tau AB(w^{n+1}-w^{n}). \ \ \ \ \ (3)$

This is a kind of implicit Euler method with an additional small term ${\tau AB(w^{n+1}-w^{n})}$ . From a numerical point of it has one advantage over the implicit Euler method: As equations (1) and (2) show, one does not need to invert ${I-\tau(A+B)}$ in every iteration, but only ${I-\tau A}$ and ${I-\tau B}$ . Remember, this was in 1950s, and solving large linear equations was a much bigger problem than it is today. In this specific case of the heat equation, the operators ${A}$ and ${B}$ are in fact tridiagonal, and hence, solving with ${I-\tau A}$ and ${I-\tau B}$ can be done by Gaussian elimination without any fill-in in linear time (read Thomas algorithm). This is a huge time saver when compared to solving with ${I-\tau(A+B)}$ which has a fairly large bandwidth (no matter how you reorder).

How do they prove convergence of the method? They don’t since they wanted to solve a parabolic PDE. They were after stability of the scheme, and this can be done by analyzing the eigenvalues of the iteration. Since the matrices ${A}$ and ${B}$ are well understood, they were able to write down the eigenfunctions of the operator associated to iteration (3) explicitly and since the finite difference approximation is well understood, they were able to prove approximation properties. Note that the method can also be seen, as a means to calculate the steady state of the heat equation.

We reformulate the iteration (3) further to see how ${w^{n+1}}$ is actually derived from ${w^{n}}$ : We obtain

$\displaystyle (-I + \tau(A+B) - \tau^{2}AB)w^{n+1} = (-I-\tau^{2}AB)w^{n} \ \ \ \ \ (4)$

2. What about monotone inclusions?

What has the previous section to do with solving monotone inclusions? A monotone inclusion is

$\displaystyle \begin{array}{rcl} 0\in Tx \end{array}$

with a monotone operator, that is, a multivalued mapping ${T}$ from a Hilbert space ${X}$ to (subsets of) itself such that for all ${x,y\in X}$ and ${u\in Tx}$ and ${v\in Ty}$ it holds that

$\displaystyle \begin{array}{rcl} \langle u-v,x-y\rangle\geq 0. \end{array}$

We are going to restrict ourselves to real Hilbert spaces here. Note that linear operators are monotone if they are positive semi-definite and further note that monotone linear operators need not to be symmetric. A general approach to the solution of monotone inclusions are so-called splitting methods. There one splits ${T}$ additively ${T=A+B}$ as a sum of two other monotone operators. Then one tries to use the so-called resolvents of ${A}$ and ${B}$ , namely

$\displaystyle \begin{array}{rcl} R_{A} = (I+A)^{-1},\qquad R_{B} = (I+B)^{-1} \end{array}$

to obtain a numerical method. By the way, the resolvent of a monotone operator always exists and is single valued (to be honest, one needs a regularity assumption here, namely one need maximal monotone operators, but we will not deal with this issue here).

The two operators ${A = \partial_{xx}}$ and ${B = \partial_{yy}}$ from the previous section are not monotone, but ${-A}$ and ${-B}$ are, so the equation ${-Au - Bu = 0}$ is a special case of a montone inclusion. To work with monotone operators we rename

$\displaystyle \begin{array}{rcl} A \leftarrow -A,\qquad B\leftarrow -B \end{array}$

and write the iteration~(4) in terms of monotone operators as

$\displaystyle \begin{array}{rcl} (I + \tau(A+B) + \tau^{2}AB)w^{n+1} = (I+\tau^{2}AB)w^{n}, \end{array}$

i.e.

$\displaystyle \begin{array}{rcl} w^{n+1} = (I+\tau A+\tau B+\tau^{2}AB)^{-1}(I+\tau AB)w^{n}. \end{array}$

Using ${I+\tau A+\tau B + \tau^{2}A = (I+\tau A)(I+\tau B)}$ and ${(I+\tau^{2}AB) = (I-\tau B) + (I + \tau A)\tau B}$ we rewrite this in terms of resolvents as

$\displaystyle \begin{array}{rcl} w^{n+1} & = &(I+\tau B)^{-1}[(I+\tau A)^{-1}(I-\tau B) + \tau B]w^{n}\\ & =& R_{\tau B}(R_{\tau A}(w^{n}-\tau Bw^{n}) + \tau Bw^{n}). \end{array}$

This is not really applicable to a general monotone inclusion since there ${A}$ and ${B}$ may be multi-valued, i.e. the term ${Bw^{n}}$ is not well defined (the iteration may be used as is for splittings where ${B}$ is monotone and single valued, though).

But what to do, when both and ${A}$ and ${B}$ are multivaled? The trick is, to introduce a new variable ${w^{n} = R_{\tau B}(u^{n})}$ . Plugging this in throughout leads to

$\displaystyle \begin{array}{rcl} R_{\tau B} u^{n+1} & = & R_{\tau B}(R_{\tau A}(R_{\tau B}u^{n}-\tau B R_{\tau B}u^{n}) + \tau B R_{\tau B}u^{n}). \end{array}$

We cancel the outer ${R_{\tau B}}$ and use ${\tau B R_{\tau B}u^{n} = u^{n} - R_{\tau B}u^{n}}$ to get

$\displaystyle \begin{array}{rcl} u^{n+1} & = & R_{\tau A}(2R_{\tau B}u^{n} - u^{n}) + u^{n} - R_{\tau B}u^{n} \end{array}$

and here we go: This is exactly what is known as Douglas-Rachford method (see the last version of the iteration in my previous post). Note that it is not ${u^{n}}$ that converges to a solution, but ${w^{n} = R_{\tau B}u^{n}}$ , so it is convenient to write the iteration in the two variables

$\displaystyle \begin{array}{rcl} w^{n} & = & R_{\tau B}u^{n}\\ u^{n+1} & = & R_{\tau A}(2w^{n} - u^{n}) + u^{n} - w^{n}. \end{array}$

The observation, that these splitting method that Douglas and Rachford devised for linear problems has a kind of much wider applicability is due to Lions and Mercier and the paper is

Lions, Pierre-Louis, and Bertrand Mercier. “Splitting algorithms for the sum of two nonlinear operators.” SIAM Journal on Numerical Analysis 16.6 (1979): 964-979.

Other, much older, splitting methods for linear systems, such as the Jacobi method, the Gauss-Seidel method used different properties of the matrices such as the diagonal of the matrix or the upper and lower triangluar parts and as such, do not generalize easily to the case of operators on a Hilbert space.

February 10, 2017

How to obtain a finite primal dual gap

Posted by Dirk under Math, Optimization, Uncategorized | Tags: algorithms, duality, image processing, optimization, TV |
Leave a Comment

Consider a convex optimization problem of the form

$\displaystyle \begin{array}{rcl} \min_{x}F(x) + G(Ax) \end{array}$

with convex ${F}$ and ${G}$ and matrix ${A}$ . (We formulate everything quite loosely, skipping over details like continuity and such, as they are irrelevant for the subject matter). Optimization problems of this type have a specific type of dual problem, namely the Fenchel-Rockafellar dual, which is

$\displaystyle \begin{array}{rcl} \max_{y}-F^{*}(-A^{T}y) - G^{*}(y) \end{array}$

and under certain regularity conditions it holds that the optimal value of the dual equals the the objective value of the primal and, moreover, that a pair ${(x^{*},y^{*})}$ is both primal and dual optimal if and only if the primal dual gap is zero, i.e. if and only if

$\displaystyle \begin{array}{rcl} F(x^{*})+G(Ax^{*}) + F^{*}(-A^{T}y^{*})+G^{*}(y^{*}) = 0. \end{array}$

Hence, it is quite handy to use the primal dual gap as a stopping criteria for iterative methods to solve these problems. So, if one runs an algorithm which produces primal iterates ${x^{k}}$ and dual iterates ${y^{k}}$ one can monitor

$\displaystyle \begin{array}{rcl} \mathcal{G}(x^{k},y^{k}) = F(x^{k})+G(Ax^{k}) + F^{*}(-A^{T}y^{k})+G^{*}(y^{k}). \end{array}$

and stop if the value falls below a desired tolerance.

There is some problem with this approach which appears if the method produces infeasible iterates in the sense that one of the four terms in ${\mathcal{G}}$ is actually ${+\infty}$ . This may be the case if ${F}$ or ${G}$ are not everywhere finite or, loosely speaking, have linear growth in some directions (since then the respective conjugate will not be finite everywhere). In the rest of the post, I’ll sketch a general method that can often solve this particular problem.

For the sake of simplicity, consider the following primal dual algorithm

$\displaystyle \begin{array}{rcl} x^{k+1} & = &\mathrm{prox}_{\tau F}(x^{k}-A^{T}y^{k})\\ y^{k+1} & = &\mathrm{prox}_{\sigma G^{*}}(y^{k}+\sigma A(2x^{k+1}-x^{k})) \end{array}$

(also know as primal dual hybrid gradient method or Chambolle-Pock’s algorithm). It converges as soon as ${\sigma\tau\leq \|A\|^{-2}}$ .

While the structure of the algorithm ensures that ${F(x^{k})}$ and ${G^{*}(y^{k})}$ are always finite (since always ${\mathrm{prox}_{F}(x)\in\mathrm{dom}(F)}$ ), is may be that ${F^{*}(-A^{T}y^{k})}$ or ${G(Ax^{k})}$ are indeed infinite, rendering the primal dual gap useless.

Let us assume that the problematic term is ${F^{*}(-A^{T}y^{k})}$ . Here is a way out in the case where one can deduce some a-priori bounds on ${x^{*}}$ , i.e. a bounded and convex set ${C}$ with ${x^{*}\in C}$ . In fact, this is often the case (e.g. one may know a-priori that there exist lower bounds ${l_{i}}$ and upper bounds ${u_{i}}$ , i.e. it holds that ${l_{i}\leq x^{*}_{i}\leq u_{i}}$ ). Then, adding these constraints to the problem will not change the solution.

Let us see, how this changes the primal dual gap: We set ${\tilde F(x) = F(x) + I_{C}(x)}$ where ${C}$ is the set which models the bound constraints. Since ${C}$ is a bounded convex set and ${F}$ is finite on ${C}$ , it is clear that

$\displaystyle \begin{array}{rcl} \tilde F^{*}(\xi) = \sup_{x\in C}\,\langle \xi,x\rangle - F(x) \end{array}$

is finite for every ${\xi}$ . This leads to a finite duality gap. However, one should also adapt the prox operator. But this is also simple in the case where the constraint ${C}$ and the function ${F}$ are separable, i.e. ${C}$ encodes bound constraints as above (in other words ${C = [l_{1},u_{1}]\times\cdots\times [l_{n},u_{n}]}$ ) and

$\displaystyle \begin{array}{rcl} F(x) = \sum_{i} f_{i}(x_{i}). \end{array}$

Here it holds that

$\displaystyle \begin{array}{rcl} \mathrm{prox}_{\sigma \tilde F}(x)_{i} = \mathrm{prox}_{\sigma f_{i} + I_{[l_{i},u_{i}]}}(x_{i}) \end{array}$

and it is simple to see that

$\displaystyle \begin{array}{rcl} \mathrm{prox}_{\sigma f_{i} + I_{[l_{i},u_{i}]}}(x_{i}) = \mathrm{proj}_{[l_{i},u_{i}]}\mathrm{prox}_{\tau f_{i}}(x_{i}), \end{array}$

i.e., only uses the proximal operator of ${F}$ and project onto the constraints. For general ${C}$ , this step may be more complicated.

One example, where this makes sense is ${L^{1}-TV}$ denoising which can be written as

$\displaystyle \begin{array}{rcl} \min_{u}\|u-u^{0}\|_{1} + \lambda TV(u). \end{array}$

Here we have

$\displaystyle \begin{array}{rcl} F(u) = \|u-u^{0}\|_{1},\quad A = \nabla,\quad G(\phi) = I_{|\phi_{ij}|\leq 1}(\phi). \end{array}$

The guy that causes problems here is ${F^{*}}$ which is an indicator functional and indeed ${A^{T}\phi^{k}}$ will usually be dual infeasible. But since ${u}$ is an image with a know range of gray values one can simple add the constraints ${0\leq u\leq 1}$ to the problem and obtains a finite dual while still keeping a simple proximal operator. It is quite instructive to compute ${\tilde F}$ in this case.

November 7, 2016

Sparse reconstruction with ℓ∞ constraints

Posted by Dirk under Optimization, Sparsity, Uncategorized | Tags: algorithms, arxiv, convexity, duality, optimization, papers, sparsity |
Leave a Comment

Last week Christoph Brauer, Andreas Tillmann and myself uploaded the paper A Primal-Dual Homotopy Algorithm for ${\ell_{1}}$ -Minimization with ${\ell_{\infty}}$ -Constraints to the arXiv (and we missed being the first ever arXiv-paper with a non-trivial five-digit identifier by twenty-something papers…). This paper specifically deals with the optimization problem

$\displaystyle \begin{array}{rcl} \min_{x}\|x\|_{1}\quad\text{s.t.}\quad \|Ax-b\|_{\infty}\leq \delta \end{array}$

where ${A}$ and ${b}$ are a real matrix and vector, respecively, of compatible size. While the related problem with ${\ell_{2}}$ constraint has been addressed quite often (and the penalized problem ${\min_{x}\|x\|_{1} + \tfrac1{2\lambda}\|Ax-b\|_{2}^{2}}$ is even more popular) there is not much code around to solve this specific problem. Obvious candidates are

Linear optimization: The problem can be recast as a linear program: The constraint is basically linear already (rewriting it with help of the ones vector ${\mathbf{1}}$ as ${Ax\leq \delta\mathbf{1}+b}$ , ${-Ax\leq \delta\mathbf{1}-b}$ ) and for the objective one can, for example, perform a variable split ${x = x^{+}-x^{-}}$ , ${x^{-},x^{+}\geq 0}$ and then write ${\|x\|_{1} = \mathbf{1}^{T}x^{+}+ \mathbf{1}^{T}x^{-}}$ .
Splitting methods: The problem is convex problem of the form ${\min_{x}F(x) + G(Ax)}$ with
$\displaystyle \begin{array}{rcl} F(x) & = & \|x\|_{1}\\ G(y) & = & \begin{cases} 0 & \|y-b\|\leq\delta\\ \infty & \text{else.} \end{cases} \end{array}$

and hence, several methods for these kind of problem are available, such as the alternating direction method of multipliers or the Chambolle-Pock algorithm.

The formulation as linear program has the advantage that one can choose among a lot of highly sophisticated software tools and the second has the advantage that the methods are very easy to code, usually in just a few lines. Drawbacks are, that both methods do exploit too much specific structure of the problem.

Application of the problem with ${\ell_{\infty}}$ constraint are, for example:

The Dantzig selector, a statistical estimation technique, were the problem is
$\displaystyle \begin{array}{rcl} \min_{x}\|x\|_{1}\quad\text{s.t.}\quad \|A^{T}(Ax-b)\|_{\infty}\leq\delta. \end{array}$
Sparse dequantization as elaborated here by Jacques, Hammond and Fadili and applied here to de-quantizaton of speech signals by Christoph, Timo Gerkmann and myself.

We wanted to see if one of the most efficient methods known for sparse reconstruction with ${\ell_{2}}$ penalty, namely the homotopy method, can be adapted to this case. The homotopy method for ${\min_{x}\|x\|_{1} + \tfrac1{2\lambda}\|Ax-b\|_{2}^{2}}$ builds on the observation that the solution for ${\lambda\geq \|A^{T}b\|_{\infty}}$ is zero and that the set of solutions ${x_{\lambda}}$ , parameterized by the parameter ${\lambda}$ , is piecewise linear. Hence, on can start from ${\lambda_{0}= \|A^{T}b\|_{\infty}}$ , calculate which direction to go, how far the breakpoint is away, go there and start over. I’ve blogged on the homotopy method here already and there you’ll find some links to great software packages, but also the fact that there can be up to exponentially many breakpoints. However, in practice the homotopy method is usually blazingly fast and with some care, can be made numerically stable and accurate, see, e.g. our extensive study here (and here is the optimization online preprint).

The problem with ${\ell_{\infty}}$ constraint seems similar, especially it is clear that for ${\delta = \|b\|_{\infty}}$ , ${x=0}$ is a solution. It is also not so difficult to see that there is a piecewise linear path of solutions ${x_{\delta}}$ . What is not so clear is, how it can be computed. It turned out, that in this case the whole truth can be seen when the problem is viewed from a primal-dual viewpoint. The associated dual problem is

$\displaystyle \begin{array}{rcl} \max_{y}\ -b^{T}y - \delta\|y\|_{1}\quad\text{s.t.}\quad\|A^{T}y\|_{\infty}\leq\infty \end{array}$

and a pair ${(x^{*},y^{*})}$ is primal and dual optimal if and only if

$\displaystyle \begin{array}{rcl} -A^{T}y^{*}&\in&\text{Sign}(x^{*})\\ Ax^{*}-b & \in & \delta\text{Sign}(y^{*}) \end{array}$

(where ${\text{Sign}}$ denotes the sign function, multivalued at zero, giving ${[-1,1]}$ there). One can note some things from the primal-dual optimality system:

You may find a direction ${d}$ such that ${(x^{*}+td,y^{*})}$ stays primal-dual optimal for the constraint ${\leq\delta-t}$ for small ${t}$ ,
for a fixed primal optimal ${x_{\delta}}$ there may be several dual optimal ${y_{\delta}}$ .

On the other hand, it is not that clear which of the probably many dual optimal ${y_{\delta}}$ allows to find a new direction ${d}$ such that ${x_{\delta}+td}$ with stay primal optimal when reducing ${\delta}$ . In fact, it turned out that, at a breakpoint, a new dual variable needs to be found to allow for the next jump in the primal variable. So, the solution path is piecewise linear in the primal variable, but piecewise constant in the dual variable (a situation similar to the adaptive inverse scale space method).

What we found is, that some adapted theorem of the alternative (a.k.a. Farkas’ Lemma) allows to calculate the next dual optimal ${y}$ such that a jump in ${x}$ will be possible.

What is more, is that the calculation of a new primal or dual optimal point amounts to solving a linear program (in contrast to a linear system for ${\ell_{2}}$ homotopy). Hence, the trick of up- and downdating a suitable factorization of a suitable matrix to speed up computation does not work. However, one can somehow leverage the special structure of the problem and use a tailored active set method to progress through the path. Our numerical tests indicated is that the resulting method, which we termed ${\ell_{1}}$ -Houdini, is able to solve moderately large problems faster than a commercial LP-solver (while also not only solving the given problem, but calculating the whole solution path on the fly) as can be seen from this table from the paper:

085_runtime_l1houdini

The code of $\ell_1$ -Houdini is on Christoph’s homepage, you may also reproduce the data in the above table with your own hardware.

October 12, 2016

The randomized sparse Kaczmarz method also converges linearly

Posted by Dirk under Math, Optimization, Uncategorized | Tags: algorithms, Kaczmarz, optimization, sparsity |
Leave a Comment

Yesterday I uploaded the paper “Linear convergence of the Randomized Sparse Kaczmarz Method” by Frank Schöpfer and myself to the arXiv.

Recall that the Kaczmarz method for linear systems

$\displaystyle \begin{array}{rcl} Ax&=&b \end{array}$

iterates

$\displaystyle \begin{array}{rcl} x^{k+1} &=& x^{k} - \frac{\langle a_{i},x_{k}\rangle-b_{i}}{\|a_{i}\|^{2}}a_{i} \end{array}$

where ${a_{i}}$ is the ${i}$ -th row of ${A}$ , ${b_{i}}$ is the ${i}$ th entry of ${b}$ and the index ${i}$ is chosen according to some rule. We could, e.g. choose the rows in a cyclic manner, i.e. starting with the first row, proceeding to the next row and once we came to the last row we would start over from the first again. It is known (and probably proved by Kaczmarz himself) that the method converges to a solution of ${Ax=b}$ whenever the system has a solution. Moreover, it easy to see that we converge to the minimum norm solution in case of underdetermined systems when the method is initialized with zero. This is due to the fact that the whole iteration takes place in the range space of ${A^{T}}$ .

In this and this paper we proposed a simple modification of the Kaczmarz method, that makes it converge to sparse solutions. The modification is simply

$\displaystyle \begin{array}{rcl} z^{k+1} & = & z^{k} - \frac{\langle a_{i},x_{k}\rangle-b_{i}}{\|a_{i}\|^{2}}a_{i}\\ x^{k+1}& = & S_{\lambda}(z^{k+1}) \end{array}$

where ${S_{\lambda}(x) = \max(|x|-\lambda,0)\text{sign}(x)}$ is the soft thresholding function. In this paper we proved that this variant converges, when initialized with zero and for a consistent system, to the solution of

$\displaystyle \begin{array}{rcl} \min_{x}\lambda\|x\|_{1} + \tfrac12\|x\|_{2}^{2},\quad\text{s.t.}\quad Ax=b. \end{array}$

For not too small values of ${\lambda}$ this is indeed a sparse solution of ${Ax=b}$ and Frank also proved that there is a threshold such that for ${\lambda}$ larger than this threshold the solution is also the minimum ${\ell^{1}}$ solution.

In general, convergence rates for the Kaczmarz method (and its sparse variant) are hard to prove. To convince oneself of this fact note that the convergence speed can drastically change if the rows of the system are reordered. The situation changes if one uses a randomization of the sparse Kaczmarz method where, in each iteration, a row is chose at random. Strohmer and Vershynin proved that this randomization leads to linear convergence rate. In the above mentioned paper we were able to prove the same result for the randomized sparse Kaczmarz method. While this sounds like an obvious generalization, the methods we use are totally different. While the linear convergence of the randomized Kaczmarz method can be proven in a few lines(well, probably one page) using only very basic tools, we need, among other things, quite intricate error bounds for Bregman projections w.r.t. piecewise linear-quadratic convex functions.

In fact, the linear convergence can be observed in practice and we illustrate the usefulness of the randomization and also the “sparsification” on some examples in the paper. For example the following figure shows the decay of the residual for the the randomized Kaczmarz method (black), the randomized sparse Kaczmarz method (red) and the randomized sparse Kaczmarz method with exact steps (green) which I did not explain here.

083_randomizedkaczmarz-figure0

More details are in the paper…

July 1, 2011

Solving Basis Pursuit itself with ISAL1

Posted by Dirk under Math, Optimization, Signal and image processing, Sparsity | Tags: algorithms, Basis pursuit, sparsity, subgradient method |
[4] Comments

Recently Andreas Tillmann presented the poster “An Infeasible-Point Subgradient Algorithm and a Computational Solver Comparison for l1-Minimization” at SPARS11. This poster summarized some results of the project SPEAR on sparse exact and approximate recovery of Marc Pfetsch an myself. We used this as an opportunity to release a draft of the accompanying paper with the same title. Although this draft is not totally ready to be submitted yet, I already summarize its content here.

Is this paper we considered the Basis Pursuit problem (beware: the linked Wikipedia page is stub at this time) from a purely optimization point of view. The Basis Pursuit problem is: For given matrix ${A\in{\mathbb R}^{m\times n}}$ (with ${m<n}$ ) and a vector ${b\in{\mathbb R}^m}$ , find the solution to

$\displaystyle \min_{x} \|x\|_1\quad\text{s.t.}\quad Ax = b. \ \ \ \ \ (1)$

Hence, we mainly neglected all its interesting features of reproducing the sparsest solution of an underdetermined linear system and so on and solely concentrated on its solution as an optimization problem.

The paper has three somehow separated contributions:

The new algorithm ISAL1: The problem (1) is a convex nonsmooth constrained optimization problem. Marc and Andreas are optimizers and they wondered how the most basic method for this class of problems would perform: The projected subgradient method: For solving
$\displaystyle \min_x f(x)\quad\text{s.t.}\quad x\in C$

take steps along some negative subgradient and project back to ${C}$ : ${x^{k+1} = P_C(x^k - \alpha_k h^k)}$ . For (1) subgradients are readily available, e.g. ${h^k = \text{sgn}(x^k)}$ (taken coordinate-wise). However, projecting onto the constraint ${Ax=b}$ is not too easy. Denoting the projection simply by ${P}$ , we can give a closed form expression (assuming that ${A}$ has full rank) as

$\displaystyle P(z) = (I - A^T (AA^T)^{-1} A) z + A^T(AA^T)^{-1}b,$

this has the drawback that one needs to explicitly invert a matrix (which, however, is just ${m\times m}$ and hence, is usually not too large since we assume ${m<<n}$ ). However, we proposed replace the exact projection by an approximate one: In each step we solve for the projection by a truncated conjugate gradient method. While we expected that one should increase the accuracy of the approximate projection by increasing the number of CG-steps during iteration, surprisingly that is not true: Throughout the iteration, a fixed small number of iterations (say ${5}$ for matrices of size ${1000\times 4000}$ but mainly independently of the size) suffices to obtain convergence (and especially feasibility of the iterates). In this paper we give a proof of convergence of the methods under several assumptions on the step-sizes and projection accuracies building on our previous paper in which we analyzed this method in the general case. Moreover, we described several ways to speed up and stabilize the subgradient method. Finally, we called this method “Infeasible point subgradient algorithm for ${\ell^1}$ ”: ISAL1. A Matlab implementation can be found and the SPEAR website.
HSE, the heuristic support evaluation: That’s a pretty neat device which can be integrated in any Basis Pursuit solver (beware: not Basis Pursuit denoising; we want the equality constraint). The idea is based on the following small lemma:

Lemma 1 A feasible vector ${\bar x}$ (i.e. ${A\bar x = b}$ ) is optimal for (1) if and only if there is ${w\in{\mathbb R}^m}$ such that ${A^Tw \in\partial\|\bar x\|_1}$ .

The proof basically consists of noting that the normal cone on the constraint ${\{Ax=b\}}$ is the image space of ${A^T}$ and hence, the condition is equivalent to saying that this normal cone intersects the subgradient ${\partial\| \bar x\|_1}$ which is necessary and sufficient for ${\bar x}$ being optimal. In practice the HSE does the following:
- deduce candidate (approx.) support ${S}$ from a given ${x}$
- compute approximate solution ${\hat{w}}$ to ${A_{S}^T w = \text{sgn}(x_{S})}$ by ${w = (A_S^T)^\dagger\text{sgn}(x_S)}$ with the help of CG
- if ${\|A^T \hat{w}\|_\infty \approx 1}$ check existence of a ${\hat{x}}$ with ${A_{S} \hat{x}_{S} = b}$ and ${\hat{x}_i = 0}$ ${\forall\, i \notin S}$
- if that ${\hat x}$ exists, check if the relative duality gap ${(\|\hat{x}\|_1 + b^T (-\hat{w}))/\|\hat{x}\|_1}$ is small and return “success” if so, i.e. take ${\hat x}$ as an optimal solution
Again, CG usually performs great here and only a very few iterations (say ${5}$ ) are needed. In practice this methods did never return any vector ${\hat x}$ marked as optimal which was wrong.
Computational comparison: We faced the challenge of a computational comparison for Basis Pursuit solvers.
The first step was, to design a testset. We constructed 100 matrices (74 of which are dense, 26 are sparse) by several constructions and concatenations (see Section 5 in the draft). More complicated was the construction of appropriate right hand sides. Why? Because we wanted to have unique solutions! That is, because we wanted to have the norm difference between optimal and computed solution as a measure for both optimality and feasibility. In the first place we used the ERC due to Joel Tropp (e.g. described in this blog post of Bob’s blog). However, this does not only guarantee uniqueness of solutions but also that the minimum -norm solution is also the sparsest. Since that is probably too much to have for solutions (i.e. they have to be very sparse) we constructed some more right hand sides using L1TestPack: Construct an such that there is such that and use . This also leads to unique solutions for Basis Pursuit if is injective when restricted to the columns which related to the entries in which but allows for much larger supports.For the results of the comparison of ISAL1, SPGL1, YALL1, ${\ell^1}$ -MAGIC, SparseLAB, the homotopy solves of Salman Asif and CPLEX check the paper. However, some things are interesting:
1. homotopy is the overall winner (which is somehow clear for the instances constructed with ERC but not for others). Great work Salman!
2. ISAL1 is quite good (although it is the simplest among all methods).
3. HSE works great: Including it e.g. in SPGL1 produces “better” solution in less time.
4. CPLEX is remarkably good (we used the dual simplex). So: How does it come that so many people keep saying that standard LP-solves do not work well for Basis Pursuit? That is simply not true for the dual simplex! (However, the interior point methods in CPLEX was not competitive at all.)
We plan to make a somehow deeper evaluation of our computational results before submitting the paper to have some more detailed conclusions on the performance of the solvers an different instances.