A Field Guide to Forward-Backward Splitting with a FASTA Implementation

Introduction

A large number of non-differentiable and constrained convex optimization problems have the following form:

In many situations, the function $g$ is neither differentiable nor even finite-valued, in which case the problem (1) cannot be minimized using simple gradient-descent methods. However, for a large class of functions $g$ that arise in practice, one can efficiently compute the so-called proximal operator

The proximal operator finds a point close to the minimizer of $g$ without straying too far from a starting point $z$, and is often referred to as a backward (or implicit) gradient-descent step with stepsize $\tau.$ If the proximal operator (2) can be evaluated easily, then one can solve (1) efficiently using the Forward-Backward Splitting (FBS) method (also known as the proximal gradient method). Put simply, FBS can handle non-differentiable objectives and convex constraints while maintaining the simplicity of gradient-descent methods.

Due to the vast applications of FBS and its utility for sparse coding and regression, many variants of FBS have been developed to improve performance and ease of use. In its raw form, FBS requires the user to choose a number of convergence parameters that strongly effect both the performance and reliability of the algorithm. These include stepsizes, stopping condition parameters, acceleration schemes, stability conditions, and initialization. With the right modifications, FBS can be performed without substantial oversight from the user.

This article introduces and reviews FBS from a practical point of view. While this is not the first review article written on FBS, it differs from existing review articles in that it focuses on practical implementation issues. For and excellent theoretical review of FBS, we refer to .

In Section 2, we introduce the forward-backward splitting method and discuss its convergence behavior. Numerous example problems are discussed in Section 3, and for each we detail the formulation and solution via FBS. In Section 4, we discuss practical issues related to the implementation of FBS. The performance of variants of FBS on different test problems is explored in Section 6. Practical issues of FBS are discussed in Section 4 and are incorporated into a new reference implementation of FBS, called FASTA (short for Fast Adaptive Shrinkage/Thresholding Algorithm). FASTA provides a simple interface for applying forward-backward splitting to a broad range of optimization problems.

2. A Word on Notation

Forward-Backward Splitting

Forward-Backward Splitting is a two-stage method that addresses each term in (1) separately. The FBS method is listed in Algorithm 1.

Let’s examine each step of the algorithm in detail. Line (3) performs a simple forward gradient descent step on $f$. This step begins at iterate $x^{k},$ and then moves in the direction of the (negative) gradient of $f,$ which is the direction of steepest descent. The scalar $\tau^{k}$ is the stepsize which controls how far the iterate moves along the gradient direction during iteration $k.$

Equation (4) is called the proximal step, or backward gradient descent step. To understand this terminology, we examine the proximal operator (2). Any $x^{\star}$ that minimizes (2) must satisfy the optimality condition

where $G\in\partial g(x^{\star})$ is some sub-gradient (generalized derivative) of $g.$ Note that when $g$ is differentiable we simply have $G=\nabla g(x^{\star}).$ Equation (5) rearranges to

This shows that $x^{\star}$ is obtained from $z$ by marching down the sub-gradient of $g.$ For this reason, the proximal operator performs a gradient descent step. Because the sub-gradient $G$ is evaluated at the final point $x^{\star}$ rather than the starting point $z,$ this is called backward gradient descent.

Equation (5) is equivalent to the set inclusion $0\in\tau\partial g(x^{\star})+(x^{\star}-z),$ which rearranges to

For this reason, the proximal operator (2) is sometimes written

where $J_{\tau\partial g}=(\tau^{k}\partial g+I)^{-1}$ is the resolvent operator of $\tau\partial g.$ The resolvent is simply another way to express the proximal operator. In plain terms, the proximal/resolvent operator, when applied to $z,$ performs a backward gradient descent step starting at $z$.

Algorithm 1 alternates between forward gradient descent on $f$, and backward gradient descent on $g$. The use of a backward step for $g$ is advantageous in several ways. First, it may be difficult to choose a sub(gradient) of $g$ in cases where the sub-gradient $\partial g$ has a complex form or is not unique. In contrast, it can be shown that problem (4) always has a unique well-defined solution , and (as we will see later) it is often possible to solve this problem in simple closed form. Second, the backward step has an important effect on the convergence of FBS, which is discussed in the next section.

The use of backward (as opposed to forward) descent for the second step of FBS is needed to guarantee convergence. To see this, let $x^{\star}$ denote a fixed point of the FBS iteration. Such a point satisfies

for some $G(x^{\star})\in\partial g(x^{\star}).$ This simplifies to

which is the optimality condition for (1). This simple argument shows that a vector is a fixed-point of the FBS iteration if and only if it is optimal for (1). This equation has a simple interpretation: when FBS is applied to an optimal point $x^{\star},$ the point is moved to a new location by the forward descent step, and the backward descent step puts it back where it started. Both the forward and backward step agree with one another because they both rely on the gradients evaluated at $x^{\star}.$

This simple fixed point property is not enough to guarantee convergence. FBS is only convergent when the stepsize sequence $\{\tau^{k}\}$ satisfies certain stability bounds. An important property of FBS is that this stability condition does not depend on $g,$ but rather on the curvature of $f.$ In the case of a constant stepsize $\tau^{k}=\tau,$ FBS is known to converge for

where $L(\nabla f)$ is a Lipschitz constant of $\nabla f$ (i.e., $\|\nabla(x)-\nabla(y)\|<L\|x-y\|$ for all $x,\,y$). In many applications $f=\frac{1}{2}\|Ax-b\|^{2}$ for some matrix $A$ and vector $b.$ In this case, $L(\nabla f)$ is simply the spectral radius of $A^{T}A.$

For non-constant stepsizes, it is known that convergence is guaranteed if the stepsizes satisfy $0<l<\tau^{k}<u<2/L(\nabla f)$ for some upper bound $u$ and lower bound $l$ (see theorem 3.4 in for this result and its generalization).

In practice, one seldom has accurate knowledge of $L(\nabla f),$ and the best stepsize choice depends on both the problem being solved and the error at each iteration. For this reason, it is better in practice to choose the sequence $\{\tau^{k}\}$ adaptively and enforce convergence using backtracking rules rather than the explicit stepsize restriction (6). These issues will be discussed in depth in Section 4.4.

Applications of Forward-Backward Splitting

In this section we study a variety of problems, and discuss how they are formulated and solved using forward-backward splitting. FBS finds use in a large number of fields, including machine learning, signal and image processing, statistics, and communication systems. Here, we briefly discuss a small subset of potential applications. In Section 6 we present numerical experiments using these test problems.

The projected gradient (PG) method is a special case of FBS involving a convex constraint. Suppose we are interested in solving

for some convex set $\mathcal{C}$. We can rewrite this problem in the form (1) using the (non-smooth) characteristic function $\mathcal{X}_{\mathcal{C}}(x)$ of the set $\mathcal{C}$, which is zero for $x\in\mathcal{C}$ and infinity otherwise. As a consequence, the problem

is then equivalent to (1) with $g(x)=\mathcal{X}_{\mathcal{C}}(x)$. To apply FBS, we must evaluate the proximal operator of $\mathcal{X}_{\mathcal{C}}$:

The solution to (8) is the element of $\mathcal{C}$ closest to $z$; this is simply the orthogonal projection of $z$ onto the set $\mathcal{C}$. The resulting PG method was originally studied by Goldstein, Levitin, and Polyak .

1.2. Lasso Regression

One of the earliest sparse regression tools from high-dimensional statistics is the Lasso regression, which is easily written in the form (7).

An important application of PG methods as in (7) is the Lasso regression problem , which is an important sparse regression tool. The Lasso regression is defined as follows:

In statistics, problem (10) is used to find sparse solutions to under-determined least squares problems. Problem (10) is called basis pursuit denoising (BDPN) in the context of compressive sensing and sparse signal recovery .

Suppose the vector $b\in\{0,1\}^{M}$ contains binary entries representing the outcomes of random Bernoulli trials. The success probability of the $i$th entry is $P(b_{i}=1\,|\,x)=e^{A_{i}x}/(1+e^{A_{i}x})$ where $A_{i}$ denotes the $i$th row of $A.$ One is then interested in solving the so-called sparse logistic regression problem

with the logit penalty function defined as

Both problems, (10) and (11), can be solved using FBS with $g(x)=\|x\|_{1}.$ The proximal operator of $g$ is given by the well-known shrinkage operator $\operatorname*{shrink}(z,\mu\tau),$ whose $i^{\text{th}}$ element is obtained as

1.4. Multiple Measurement Vector (MMV)

In Section (3.1.3) we sought sparse vectors satisfying $Ax\approx b$ for some matrix $A$ and measurement vector $b.$ Suppose now that we have a matrix $B$ containing many measurement vectors, and we want to find a sparse matrix $X$ satisfying $AX\approx B.$ If we further suppose that all columns of $X$ must have the same sparsity pattern, then we arrive at the Multiple Measurement Vector (MMV) problem . To formulate MMV using convex optimization, we need the group sparsity prior

where $X_{i}$ denotes the $i$th row of $X.$

where the scalar $\mu$ controls the strength of the penalty. This is easily solved using FBS. The forward step is simply

The backward step requires the proximal operator of $\mu\operatorname*{MMV}(\cdot),$ which can be evaluated row-by-row. The $i$th row is simply

1.5. Democratic Representations

1.6. Low-Rank (1-bit) Matrix Completion

One prominent formulation of matrix completion takes the form

where $\|X\|_{*}$ is the low-rank inducing nuclear norm of the matrix $X$ and $L(\cdot,Y)$ is a loss function that depends on the model of the observed data $Y$ .

Suppose that $Y$ contains noisy and punctured (or missing) data of the entries of $\hat{X}$, and let $\Omega$ denote the set of indices that have been observed. In this case, one can use

as an appropriate loss function. Another loss function arises in 1-bit matrix completion, which can be viewed as an instance of low-rank logistic regression . Each observation $Y_{\alpha}$ is assumed to be a Bernoulli random variable with success probability $P(Y_{\alpha}=1|\hat{X}_{\alpha})=e^{\hat{X}_{\alpha}}/(1+e^{\hat{X}_{\alpha}}).$ For this case, the appropriate loss function is again the logit function:

For both quadratic and logit link functions, (14) can be solved using FBS. Specifically, one needs to compute the proximal operator of the nuclear norm, defined as

whose solution is given by the matrix $U\operatorname*{shrink}(S,\mu\tau)V^{T}$, where $Z=USV^{T}$ is the singular value decomposition (SVD) of $Z$ .

2. More Complex Applications

Sometimes an objective function does not immediately decompose into a smooth part and a “simple” part for which we can evaluate the proximal operator in closed form. In this case, it is often possible to re-formulate the problem in such a way that FBS is easily applied. In this section, we look at two such problems (total-variation denoising and support vector machines) that can be reformulated using duality and then solved easily and efficiently using FBS. We also look at an example where the convex relaxation of a problem is easily solved using FBS.

Given a 2-dimensional image $u,$ the intensity of the pixel in the $i$th row and $j$th column is denoted $u_{ij}.$ The discrete gradient of $u$ is denoted $\nabla u.$ At each point in the image, the gradient is given by $(\nabla u)_{ij}=(u_{i+1,j}-u_{i,j},u_{i,j+1}-u_{i,j})^{T}$ and contains the discrete derivatives in the row and column directions.

Given a noise-contaminated image $f,$ one can construct a denoised image by solving

where $\mu$ is a parameter to control the level of smoothing, and

denotes the total-variation of $u$ . Rather than imposing sparsity on $u$ itself, the total variation regularizer imposes sparsity on the gradient of $u.$ As a result, problem (16) finds a piecewise-constant approximation to $f.$

Equation (18) follows from the Cauchy-Swartz inequality, which states that $\langle x_{ij},(\nabla u)_{ij}\rangle\leq\|x_{ij}\|\|(\nabla u)_{ij}\|.$ Equality is attained by choosing $x_{ij}$ parallel to $(\nabla u)_{ij}$. If we further choose $x_{ij}$ to have unit norm, then we arrive at (18). If we apply (18) to (17) we get ${|\nabla u|=\max_{\|x\|_{\infty}\leq 1}\langle x,\nabla u\rangle.}$ We now have

The inner minimization in (19) is now differentiable. For a given $x,$ the minimal value of $u$ satisfies $u=f+\mu\nabla\cdot x,$ where $\nabla\cdot x$ is the discrete divergence (which is the negative adjoint of the gradient operator). At pixel $ij,$ this operator takes on the scalar value ${(\nabla\cdot x)_{ij}=x^{r}_{i,j}-x^{r}_{i-1,j}+x^{c}_{i,j}-x^{c}_{i,j-1}.}$ If we plug the optimal value of $u$ into (19) and simplify, we see that the optimal value of $x$ is given by:

This is simply a quadratic minimization with an infinity-norm constraint. Problem (20) is known as the dual form of (16). This problem can be solved using FBS as in Section (3.1.1). The resulting algorithm alternately performs gradient descent steps on (20), and then re-projects the result back into the infinity-norm ball using the formula ${x_{ij}\leftarrow x_{ij}/\max\{\|x_{ij}\|,1\}.}$ Once problem (20) has been solved, the optimal (denoised) image $u^{\star}=f+\mu\nabla\cdot x^{\star}$ is computed.

This approach to total-variation minimization was first taken in using constant stepsize parameters. A similar approach was taken using accelerated variants of FBS in .

2.2. Support Vector Machines

The hinge loss function penalizes points that either lie on the “wrong” side of the hyper-plane $w^{T}x=0$ or that lie too close to this hyperplane.

where the constraints on the vector $x$ are interpreted element-wise, $L$ is a diagonal matrix with $L_{i,i}=l_{i},$ and $1_{n}$ is a vector of 1’s. This objective is quadratic in $w$, and the optimal value of $w$ is given by $D^{T}Lx.$ If we plug this value in for $w$ in (22) and multiply by -1 to convert the maximization for $x$ into a minimization, we arrive at the dual problem

This is a simple quadratic minimization with constraints of the form $0\leq x_{i}\leq C.$ The proximal operator corresponding to these constraints is the projection $\operatorname{prox}(z,t)=\min\{\max\{z,0\},C\}.$ Once the dual problem (23) is solved using FBS, the primal solution is recovered from the formula $w=D^{T}Lx.$ SVM’s and the dual minimization are reviewed in .

2.3. Phase Retrieval and Rank Minimization

A variety of non-convex problems can be relaxed into convex problems involving symmetric positive semi-definite (SPDP) matrices. This idea was popularized by , in which it was shown that the NP-complete MaxCut problem has a convex relaxation with the form of a semi-definite program. Here, we consider the more recent PhaseLift algorithm, which can be used for phase retrieval applications .

Recovery of $x$ from the phase-less measurements $b$ requires the solution of a system of non-linear equations. However, one can relax the equations into a convex problem by defining $A_{i}=a_{i}a_{i}^{T},$ and observing that $|\langle a_{i},x\rangle|^{2}=\langle A_{i},xx^{T}\rangle=b_{i}.$ Letting $\mathcal{A}(xx^{T})_{i}=\langle A_{i},xx^{T}\rangle,$ the set of phase-less measurements can be written compactly as $\mathcal{A}(xx^{T})=b.$ Finally, note that any SPSD rank-1 matrix $X$ can be factorized as $X=xx^{T}.$ This observation allows us to pose the recovery problem in the following form:

which can be solved using FBS. The relevant proximal operator corresponds to

whose solution is given by the matrix $U\operatorname*{shrink}(\Lambda,\mu\tau)U^{T}$, where $Z=U\Lambda U^{T}$ is the eigenvalue decomposition of $Z$.

3. Non-convex Problems

All applications discussed so far involved convex optimization. In practice, FBS also works quite well when applied to non-convex problems. However, in this case there are no theoretical guarantees that the method with converge to a global minimizer, or that the algorithm will even converge at all. That being said, it is often the case for non-convex problems that no known method can guarantee optimality in polynomial time, and so this lack of theoretical guarantees makes FBS no weaker than other options.

When solving non-convex problems the user should be cautious of several things. First, unlike in the convex case, the final result will be sensitive to the initial iterate. For this reason, it is important to either initialize the methods with a good approximate solution, or else generate numerous solutions from different random initializations and choose the solution with the lowest objective value. Second, because the method is no longer guaranteed to converge, the user must put an intelligent limit on the number of iterations. In some situations, it may even be best to run the algorithm for a pre-determined number of iterations rather than using one of the precision-based termination conditions to be discussed in Section 4.6.

where the constraints on $W$ and $C$ are interpreted element-wise.

This problem is known as non-negative matrix factorization (NMF) . While various methods have been proposed to solve this problem , forward backward splitting remains one of the simplest. The gradient of (26) is easily obtained using the chain rule, resulting in the forward step

The backward step is a trivial projection onto the set of non-negative matrices – i.e., negative entries in the matrices are replaced by zeros.

3.2. Max-Norm Optimization and the Max Cut Problem

where $X_{i}$ is the $i$th row of $X.$ It can shown that the max-norm of a matrix is a good approximation to the nuclear norm (up to a constant factor), and can be used to promote low-rank solutions when the nuclear norm is not available because of the expense of computing singular value decompositions .

An important application of the max-norm is for solving the max-cut problem on a graph . Consider a graph with $N$ vertices and an edge of weight $w_{ij}$ between vertex $i$ and $j.$ An assignment of a binary label $x_{i}\in\{-1,1\}$ to each vertex is called a “cut.” The “value” of a cut is the sum of all edge weights connecting vertices with different labels.

Finding the graph cut with maximum value is NP-complete , however a good approximation can be found by solving

The objective in (27) differentiable, and the forward step of FBS is given by

The backward step is simply a projection onto the constraint set. In this case, we have

Bells and Whistles

FBS in its raw form suffers from numerous problems including: (i) The convergence speed of FBS depends strongly on the choice of the stepsize parameters $\{\tau_{k}\}.$ The best stepsize choices may not be intuitively obvious. (ii) Convergence is only guaranteed when the stepsizes satisfy the stability condition (6) which depends on the Lipschitz constant for $\nabla f.$ In real applications one often has no knowledge of this Lipschitz constant and no practical way to estimate it. (iii) To automatically stop the FBS iteration, a good measure of convergence is needed.

In this section, we discuss practical methods for overcoming these problems. We begin with adaptive schemes for automatic stepsize selection. We then discuss backtracking methods that can guarantee stability even when the user has no explicit knowledge of $L(\nabla f).$ Finally, we discuss stopping conditions and measures of convergence.

The efficiency of FBS (and gradient methods in general) is very sensitive to the choice of the stepsize $\tau_{k}.$ For this reason, much work has been devoted to studying adaptive stepsize strategies. Adaptive methods automatically tune stepsize parameters in real time (as the algorithm runs) to achieve fast convergence. In this section, we discuss spectral (also called Barzilai-Borwein) stepsize methods, and how they can be adapted to FBS.

Before considering the full-scale FBS method for (1), we begin by considering the case $g=0.$ In this case $h=f$ and FBS reduces to simple gradient descent of the form

Spectral schemes for gradient descent were proposed by Barzilai and Borwein , who model the function $f$ as the simple quadratic function

It can be shows that the optimal stepsize choice for (29) is $\tau=1/\alpha.$ With this choice, the gradient descent method achieves a perfect minimizer of the simple quadratic (29) in one iteration. This motivates the following stepsize scheme for gradient descent: Before applying the descent step (28), approximate $f$ with a quadratic of the form (29). Then, take a step of length $\tau^{k}=1/\alpha.$

Spectral stepsizes for general FBS with arbitrary $f$ and $g$ were first studied in for use with the solver SpaRSA. The method presented here is a hybrid of the spectral scheme with the “adaptive” stepsize rule presented in .

It was observed in the introduction that the stepsize restriction for FBS depends only on $f$ and not on $g.$ Spectral methods for FBS exploit this property. The idea is to build a quadratic approximation for $f$ of the form (29) at each iteration, and then choose the optimal gradient descent stepsize $\tau^{k}=1/a.$ Let

If we assume a quadratic model for $f$ of the form (29), then we have

The function (29) is fit to $f$ using a least squares method which chooses $a$ to minimize either $\|\Delta F^{k}-a\Delta x^{k}\|^{2}$ or $\|a^{-1}\Delta F^{k}-\Delta x^{k}\|^{2}.$ Just like in the case $g=0,$ we then select a stepsize of length $\tau^{k}=1/a$. The resulting stepsize choice is given by

respectively for each least squares problem. The value $\tau_{s}^{k}$ is known as the “steepest descent” stepsize, and $\tau_{m}^{k}$ is called the “minimum residual” stepsize .

A number of variations of spectral descent are reviewed by Fletcher , several of which perform better in practice than the “classic” stepsize rules (32). We recommend the “adaptive” BB method , which uses the rule

Note that (particularly for non-convex problems), the stepsize $\tau_{m}^{k}$ or $\tau_{s}^{k}$ may be negative. If this happens, the stepsize should be discarded and replaced with its value from the previous iteration. When complex-valued problems are considered, it is important to use only the real part of the inner product in (33).

2. Preconditioning

FBS can be “preconditioned” by multiplying the gradient directions of the forward step with a symmetric positive definite matrix $\Gamma.$ The backward step must be modified as well to ensure convergence. The resulting approach, sometimes called the proximal newton method , is described in Algorithm 2.

The preconditioned backward step (35) now requires the generalized proximal operator, which involves the following norm weighted by the matrix $\Gamma^{-1}$

For general $\Gamma$ the minimization (35) may not have a closed form solution, even when $g$ is “simple.” However, when the function $g$ acts separately on each element of $x,$ adding a diagonal preconditioner is a trivial modification. For example, when $g=\mu|\cdot|$ (i.e., for the sparse least squares problems in section 3.1.3) and $\Gamma$ is diagonal, the $i$th entry of $\operatorname{prox}_{g}(z,\tau,\Gamma)$ is simply given by $\operatorname*{shrink}(z_{i},\mu\tau\Gamma_{ii}).$ Preconditioning with a diagonal $\Gamma$ is similarly easy in the case of the Support Vector Machine (Section 3.2.2).

The preconditioner has another interpretation. One can derive the preconditioned iteration (Algorithm 2) by making the change of variables $x\to Py$ for some invertible matrix $P$ and solving the problem

Each iterate $y^{k}$ can then be converted to an approximate solution to (1) by computing $x^{k}=Py^{k}.$ It can be shown that the iterates $\{x^{k}\}$ obtained through this procedure are equivalent to those obtained by Algorithm 2 with preconditioner $\Gamma=P^{T}P.$

The latter interpretation of the preconditioner is important because it suggests how to choose $\Gamma.$ When $f$ has the form

for some matrix $A$ and function $\hat{f},$ the preconditioned problem (36) involves the objective $\hat{f}(APy).$ For example, when $f(x)=\|Ax-b\|^{2},$ the problem (36) contains the term $\|APy-b\|^{2}.$ If $A$ is poorly conditioned, then one should choose $P$ to be a diagonal matrix such that $AP$ is more well conditioned. One effective method chooses the entries in $P$ so that all columns of $AP$ have unit norm. The corresponding choice of $\Gamma=P^{T}P$ is then

where $a_{i}$ denotes the $i$th column of $A.$

Note that, because of the equivalence between preconditioned FBS and the original FBS on problem (36), all the tricks described in this section (including adaptivity and backtracking) can be applied to Algorithm 2 as long as $\Gamma$ remains unchanged between iterations. In the event that $\Gamma$ changes on every iteration, the converge theory becomes somewhat more complex (see ).

3. Acceleration

Adaptive stepsize selection helps to speed the convergence of the potentially slow FBS method. Another approach to dealing with slow convergence is to use predictor-corrector schemes that “accelerate” the convergence of FBS.

While several predictor-corrector variants of FBS have been proposed (see for example ), the algorithm FISTA has become quite popular because of its lack of tuning parameters and good worst-case performance. FISTA relies on a one-size-fits-all sequence of acceleration parameters that works for any objective. FISTA is listed in Algorithm 3. Step 2 of FISTA simply performs an FBS step. Step 4 performs prediction, in which the current iterate is advanced further in the direction it moved during the previous iteration. The aggressiveness of this prediction step is controlled by the scalar parameter $\alpha^{k}.$ This acceleration parameter is updated in step 3, and increases on each iteration causing the algorithm to become progressively more aggressive. The FISTA method alternates between “predicting” an aggressive estimate of the solution, and “correcting” this estimate using FBS to attain greater accuracy.

One notable advantage of this approach is the worst-case convergence rate. It has been shown that FISTA decreases the optimality gap (the difference between the objective at iteration $k$ and the optimal objective value) with rate $O(\frac{1}{k^{2}})$ . In contrast, the worst-case performance of conventional FBS is known to be $O(\frac{1}{k}).$

Rather than requiring the user to choose a stepsize, FISTA is generally coupled with a backtracking line search such as those described in the next section.

In theory, FISTA achieves superior worst-case performance by increasing the prediction parameter $\alpha^{k}$ on every iteration. In practice, however, FISTA does not always perform well when the prediction parameter $\alpha^{k}$ becomes large. This is particularly true when a large/aggressive stepsize is used, in which case the objective values may oscillate rapidly before convergence is reached. For this reason FISTA (and other Nesterov-type methods) generally perform better when a “restart” method is used . Such methods allow $\alpha^{k}$ to increase while the performance of Algorithm 3 is good, but reset the value to $\alpha^{k}=1$ if oscillations develop at iteration $k.$

One restart strategy simply sets $\alpha^{k}=1$ on iterations where the objective value increases. However, this method requires the objective value to be calculated at each iteration, which is expensive and in general unnecessary. The authors of suggest a restart method that does not require the objective. The method restarts when the difference between iterates $x^{k}-x^{k-1}$ points in an ascent direction for the objective. It can be shown that this occurs whenever

in which case the method should be restarted . We will see in Section 4.6 that the vector $y^{k}-x^{k}$ is parallel to the gradient of the objective function. For this reason, condition (37) restarts the method whenever the change between iterates forms an acute angle with the gradient (which is the direction of steepest ascent).

4. Backtracking Line Search

Convergence of FBS can be guaranteed by enforcing the stability condition (6). In practice, though, the user generally has no knowledge of the global properties of $\nabla f,$ and so the actual value of this stepsize restriction is unknown. In Section 4.1 we discuss adaptive methods that automatically choose stepsizes for us. This does not free us from the need for stability conditions: spectral methods are, in general, not guaranteed to converge, even for convex problems .

Even without knowing the stepsize restriction (6), convergence can be guaranteed by incorporating a backtracking line search. Such methods proceed by checking a line search condition after each iteration of FBS. The search condition usually enforces that the objective has decreased sufficiently. If this condition fails to hold, then backtracking is performed – the stepsize is decreased and the FBS iteration repeated until the backtracking condition (i.e., sufficient decrease of the objective) is satisfied.

Line search methods were originally proposed for smooth problems , and later for projected gradient schemes . A backtracking scheme for general FBS was proposed by Beck and Teboulle for use with FISTA . Many authors consider these line searches to be overly conservative, particularly for poorly conditioned problems , and prefer non-monotone line search conditions . Rather than insisting upon an objective decrease on every iteration, non-monotone line search conditions allow the objective to increase within limitations.

Non-monotone line search methods are advantageous for two reasons. First, for poorly conditioned problems (i.e., objective functions with long, narrow valleys) it may be the case that iterate $x^{k+1}$ lies much closer to the minimizer than $x^{k},$ despite $x^{k+1}$ having the larger objective value. A non-monotone line search prevents such iterates from being rejected.

Second, the objective function must be evaluated every time the line-search condition is tested. For complex problems where evaluating the objective is costly and multiple backtracking steps are necessary for each iteration, the line search procedure may dominate the runtime of the method. Non-monotone line search conditions are less likely to be violated and thus backtracking terminates faster, which alleviates this computational burden.

The method proposed here is inspired by the monotone search proposed in for general FBS, and generalizes the non-monotone strategy of for SPG.

Let $M>0$ be an integer line search parameter, and define

After each step of FBS, the following line search condition is checked.

If the backtracking condition (38) fails, the stepsize is decreased until (38) is satisfied. This process is formalized in Algorithm 4. Note that Algorithm 4 always terminates because condition (38) is guaranteed to hold whenever $\tau^{k}$ is less than the reciprocal of the Lipschitz constant of $\nabla f.$

A convergence proof for FBS with the line search described in Algorithm 4 is given in the appendix.

5. Continuation

Some types of regression problems are solved very quickly when solutions are highly sparse, but become numerically difficult when solutions lack sparsity. This is sometimes the case for problem (10). When $\mu$ is large, solutions to (10) are highly sparse. For small $\mu,$ solutions may lack sparsity, resulting in slow convergence if $A$ is poorly conditioned.

Continuation methods for (10) exploit this observation by choosing a large initial value of $\mu$ and decreasing this parameter over time. This way, the solver can use the results of “easy” problems as a warm start for more difficult problems with less sparsity. Furthermore, using a warm start keeps the support of each iterate small, whereas the support might “blow up” for many iterations if difficult problems are attacked with a bad initializer.

Continuation techniques have been proposed by a number of authors (see for example ) and are a keystone component in methods such as fixed-point continuation that aim to approximate solutions to under-determined problems of the form

When choosing values for $\mu,$ it helps to observe that the solution to (10) is zero when $\mu\geq\|A^{T}b\|_{\infty}.$ For this reason it is suggested to choose an initial parameter of $\mu=\eta\|A^{T}b\|_{\infty}$ for some $\eta<1.$ After we obtain a solution using the current value of $\mu,$ we replace $\mu\leftarrow\eta\mu$ and continue iterating until the desired value of $\mu$ is reached. We suggest to choose $\eta=\frac{1}{5},$ although the practical performance of continuation is not highly sensitive to this parameter.

Continuation is mostly effective for problem (10) when $A$ is poorly conditioned, or when extremely large regularization parameters are needed. Note that continuation does not always enhance performance, and may sometimes even make performance dramatically worse.

6. Stopping Conditions

While the accuracy of FBS becomes arbitrarily good as the number of iterations approaches infinity, we must of course stop after a finite number of iterations. A good stopping criteria should be strict enough to guarantee an acceptable degree of accuracy without requiring an excessive number of iterations.

Our stopping conditions will be based on the residual, which is simply the derivative of the objective function (or a sub-gradient in the case that $g$ is non-differentiable). Because $f$ is assumed to be smooth, we can differentiate this term in the objective directly. While we may not be able to differentiate $g$, we see from equation (5) that a sub-gradient is given by $(\hat{x}^{k+1}-x^{k+1})/\tau^{k}\in\partial g(x^{k+1}).$ We now have the following formula for the residual $r^{k+1}$ at iterate $x^{k+1}:$

A simple termination rule would stop the algorithm when $\|r^{k+1}\|<tol$ for some small tolerance $tol>0$. However, this rule is problematic because it is not scale invariant. To understand what this means, consider the minimization of some objective function $h(\cdot).$ This function can be re-scaled by a factor of 1000 to obtain $\hat{h}=1000h.$ The new rescaled objective has the same minimizer as the original, however the sub-gradient $\partial\hat{h}(x^{k})$ is 1000 times larger than $\partial h(x^{k}).$ A stopping parameter $tol$ may be reasonable for minimizing $h$ but overly strict for minimizing $\hat{h},$ even though the solutions to the problems are identical. Ideally, we would like scale invariant stopping rules that treat both of these problems equally.

One way to achieve scale invariance is by replacing the residual with the relative residual. To define the relative residual, we observe that (40) is small when

In plain words, the residual (40) measures the difference between the gradient of $f$ and the negative sub-gradient of $g$. The relative residual $r^{k+1}_{r}$ measures the relative difference between these two quantities, which is given by

where $\epsilon^{r}$ is some small positive constant to avoid dividing by zero.

Another more general scale invariant stopping condition uses the normalized residual, which is given by

where the small constant $\epsilon_{n}$ prevents division by zero. Rather than being an absolute measure of accuracy, the normalized residual measures how much the approximate solution has improved relative to $x^{1}.$

Both scale invariant conditions have advantages and disadvantages. The relative residual works well for a wide range of problems and is insensitive to the initial choice of $x^{0}.$ However, the relative residual looses scale invariance when $\nabla f(x^{\star})=0$ (in which case the denominator of (42) nearly vanishes). This happens, for example, when $g$ is the characteristic function of a convex set and the constraint is inactive at the optimal point (i.e., the optimal point lies in the interior of the constraint set). In contrast, the normalized residual can be effective even if $\nabla f(x^{\star})=0.$ However, this measure of convergence is potentially sensitive to the choice of the initial iterate, and so it requires the algorithm to be initialized in some consistent way. The strictness of this condition can also depend on the problem being solved.

For general applications, we suggest a combined stopping condition that terminates the algorithm when either $r^{k+1}_{r}$ or $r^{k+1}_{n}$ gets small. In this case, $r^{k+1}_{r}$ terminates the iteration if a high degree of accuracy is attained, and $r^{k+1}_{n}$ terminates the residual appropriately in cases where $0\in\partial g(x^{\star}).$

FASTA: A Handy Forward-Backward Solver

To create a common interface for testing different FBS variants, we have created the solver FASTA (Fast Adaptive Shrinkage/Thresholding), which implements forward-backward splitting for arbitrary problems. Many improvement to FBS are implemented in FASTA including adaptively, acceleration, backtracking, a variety of stopping conditions, and more. FASTA enables different FBS variants to be compared objectively while controlling for the effects of stepsize rules, programming language, and other implementation details.

Rather than addressing the problem (1) directly, FASTA solves general problems of the form

Numerical Experiments

We compare several variants of FBS and compare their performance using the test problems described in Section 3. The variants considered here are the original FBS with constant stepsize, and the accelerated variant FISTA described in Section 4.3. We also consider the adaptive method described in Section 4.1, which is implemented in the solver FASTA. This method uses the spectral stepsize rules proposed in and .

At each iteration of the algorithms considered, the relative residual (42) was computed and used as a measure of convergence. Time trials for all methods were terminated on the first iteration that satisfied $r^{k+1}_{r}<10^{-4}.$ If this stopping condition was not reached, iterations were terminated after 1000 iterations (except for the SVM problem, which was allowed up to 5000 iterations). Time trial results are averaged over 100 random trials.

We generated problems of the form (10) using a similar procedure as for Lasso. This time, however, with the noise scaled to achieve $\text{SNR}=20$ dB. Recovery was performed with $\mu=0.1.$

We generated frames of dimension $500\times 1000$ by randomly selecting subset of rows from a unitary discrete Fourier transform matrix. The signal $b$ was a complex-valued random Gaussian vector of length 500. Equation (13) was solved with $\mu=300.$

We generated a $200\times 1000$ matrix $X$ with i.i.d. random Gaussian entries of standard deviation $10.$ The SVD of $X$ was then computed, and all but the largest 5 singular values were set to zero to obtain a rank-$5$ matrix. The sigmoidal logistic function was then applied element-wise to $X$ to obtain a matrix of probabilities; a random Bernoulli matrix $Y$ was drawn from the resulting distribution. Problem (14) was then solved using the logistic penalty function with $\mu=25.$

The $256\times 256$ Shepp-Logan phantom was constructed with pixel intensities that spanned the unit interval. Test images were then contaminated with Gaussian noise ($\sigma=0.05$). Images were denoised by solving (16) with $\mu=0.1.$

Feature vectors were generated from two classes: one class containing random Gaussian variables with mean $-1$, and one with mean $+1.$ Each trial contained 1000 feature vectors of dimension 15. A support vector machine was trained to separate the two classes with regularization parameter $C=10^{-2}.$

We generated random vectors of length $N=200$ with random Gaussian real and imaginary parts. The set $\{a_{i}\}$ was created by randomly drawing complex Gaussian vectors. The measurement vector $b$ of length $M=600$ was creating by setting $b_{i}=|\langle a_{i},x\rangle|^{2},$ and then, contaminating the resulting measurements with real-valued Gaussian noise to have $\text{SNR}=13$ dB. Problem (3.2.3) was then solved with parameter $\mu=15,$ which was chosen to be large enough that the solution matrix had rank 1.

2. Results and Discussion

For all problems considered, both the accelerated and adaptive FBS dramatically out-performed “vanilla” FBS without these modifications. However FBS does not reliably converge when adaptivity and acceleration are used simultaneously, and so the user must pick one.

To explore the efficiency of different approaches, we applied three variants of FBS to each test problem: Plain FBS, FBS with acceleration (FISTA), and FBS with adaptive stepsizes (SpaRSA). The accelerated FBS was accompanied by the restart rule (37), and all methods used backtracking line search. For each algorithm, the number of iterations (and time in seconds) needed to solve each test problem is reported in Table 1.

For most problems considered, the adaptive method out-performed the accelerated scheme by a factor of 3-to-5. For the SVM problem, the performance gap was somewhat larger. We can take a closer look at the behavior of each method with the convergence curves in Figure 1. Figure 1a shows the convergence for matrix completion, which looks typical for most problems including Lasso, penalized least-squares, logistic regression, and MMV. For such problems, we see smooth, exponential decay of the error until machine precision is reached. Note this empirical behavior is much better than the $O(1/k)$ worse-case global convergence bounds .

We also show some less-typical convergence curves, including SVM for which the performance gap between methods was extremely large. The curves for non-negative matrix completion show more irregular behavior (including oscillations) because of non-convexity. Finally, we show convergence curves for the total variation minimization in Figure 1d. For this problem, adaptive and accelerated methods were competitive. The adaptive method was superior in the low-precision regime, with the accelerated variant winning out in the high-precision regime. This is largely because total variation involves the gradient operator, which has a large condition number. Nesterov-type acceleration tends to be most effective for these types of poorly-conditioned problems, however the advantages over adaptivity are still slim in the examples shown here.

We note that the advantage of adaptivity is most pronounced for problems where $f$ is non-quadratic, i.e., for problems with a logistic data term. In this case, the Hessian of $f$ varies over the problem domain, and the optimal stepsize for FBS varies with it. For such problems, an adaptive scheme is able to effectively match the stepsize to the local structure of the objective, resulting in fast convergence.

Conclusion

The forward-backward splitting method is a surprisingly simple way to solve a wide range of optimization problems. Even seemingly complex problems involving non-differentiable objectives (total-variation, support vector machine, sparse regression), and complex constraint sets (semi-definite programing and max-norm regularization, etc…) can be reduced to a sequence of extremely simple steps.

Historically, FBS it is most commonly used for simple sparse regression problems despite its much wider applicability. In many common domains, more complex splitting methods involving Lagrange multipliers (such as ADMM and its variants ) are more commonly used. However, these method are often more complex, memory intensive, and computationally burdensome than FBS. Furthermore, FBS has a major advantage over other splitting methods – algorithm parameters like stepsizes and stopping conditions are easily automated. For other splitting methods, it is considerably more difficult to guarantee convergence for adaptive methods . For this reason it is fair to say that FBS is under-utilized for complex problems.

Appendix A Convergence Proof for Non-Monotone Line Search

We now consider the convergence of the backtracking line search discussed in section 4.4.

Suppose that FBS is applied to (1) with convex $g$ and differentiable $f$. Suppose further that $h=f+g$ is proper, lower semi-continuous, and has bounded level sets. If $\{\tau^{k}\}$ is bounded below by a positive constant and

then $\lim_{k\to\infty}h(x^{k})=h^{\star}$, where $h^{\star}$ denotes the minimum value of $h$.

From the optimality condition (4), we have $0\in\tau^{k}\partial g(x^{k+1})+x^{k+1}-\bar{x}^{k+1}$ and so

for some $G^{k+1}\in\partial g(x^{k+1})$ and $F^{k}=\nabla f(x^{k}).$ From this we arrive at

Subtracting (48) from (45) and applying (47) yields

where $\hat{h}^{k}=\max\{h^{k-1},h^{k-2},\dots,h^{k-\min\{M,k\}}\}.$ Note that $\{\hat{h}^{k}\}$ is a monotonically decreasing bounded sequence, and thus has a limit $\hat{h}^{*}.$

We aim to show that $\hat{h}^{\star}$ is the minimal value of $h$. Observe that $\hat{h}^{k}=h^{k^{\prime}}$ for some $k^{\prime}$ with $k-M\leq k^{\prime}\leq k.$ It is clear from (A) that there must exist a sub-sequence $\{x^{k(i)}\}$ with $h(x^{k(i)})=\hat{h}^{k}$ such that

(otherwise, equation (A) would imply $\hat{h}^{k}\to-\infty$ for large $k$). Note that the level sets of $h$ are assumed to be bounded, and by equation (A) $\{h(x^{k})\}$ is bounded as well. By compactness, we may assume without loss of generality that $\{x^{k(i)}\}$ is a convergent sub-sequence of iterates with limit point $x^{\star}$.

Equation (51), together with (50) and the fact that $\tau^{k}$ is bounded away from zero, implies that

Because $\nabla f$ is Lipschitz continuous and $\|x^{k(i)+1}-x^{k(i)}\|\to 0$, we also conclude that $x^{k(i)+1}\to x^{\star}$ and

Note that $G^{k(i)+1}+F^{k(i+1)}\in\partial h(x^{k(i)+1})$. Because the sub-differential of a convex function is continuousMore formally, the convex sub-differential is a multi-valued function which is upper semicontinous in a topological sense. See . , (52) implies that $0\in\partial h(x^{\star}),$ and so $x^{\star}$ is a minimizer of $h$.

We have shown that $\lim_{k\to\infty}\hat{h}^{k}=h(x^{\star})=h^{\star}.$ Because $h^{\star}\leq h(x^{k})\leq\hat{h}^{k},$ we arrive at the conclusion

References