iPiano: Inertial Proximal Algorithm for Non-Convex Optimization

Introduction

The gradient method is certainly one of the most fundamental but also one of the most simple algorithms to solve smooth convex optimization problems. In the last decades, the gradient method has been modified in many ways. One of those improvements is to consider so-called multi-step schemes . It has been shown that such schemes significantly boost the performance of the plain gradient method. Triggered by practical problems in signal processing, image processing and machine learning, there has been an increased interest in so-called composite objective functions, where the objective function is given by the sum of a smooth function and a non-smooth function with an easy to compute proximal map. This initiated the development of the so-called proximal gradient or forward-backward method , that combines explicit (forward) gradient steps w.r.t. the smooth part with proximal (backward) steps w.r.t. the non-smooth part.

In this paper, we combine the concepts of multi-step schemes and the proximal gradient method to efficiently solve a certain class of non-convex, non-smooth optimization problems. Although, the transfer of knowledge from convex optimization to non-convex problems is very challenging, it aspires to find efficient algorithms for certain non-convex problems. Therefore, we consider the subclass of non-convex problems

where $g$ is a convex (possibly non-smooth) and $f$ is a smooth (possibly non-convex) function. The sum $f+g$ comprises non-smooth, non-convex functions. Despite the non-convexity, the structure of $f$ being smooth and $g$ being convex makes the forward-backward splitting algorithm well-defined. Additionally, an inertial force is incorporated into the design of our algorithm, which we termed iPiano. Informally, the update scheme of the algorithm that will be analyzed is

where $\alpha$ and $\beta$ are the step size parameters. The term $x^{n}-\alpha\nabla f(x^{n})$ is referred as forward step, $\beta(x^{n}-x^{n-1})$ as inertial term, and $(I+\alpha\partial g)^{-1}$ as backward or proximal step.

For $g\equiv 0$ the proximal step is the identity and the update scheme is usually referred as Heavy-ball method. This reduced iterative scheme is an explicit finite differences discretization of the so-called Heavy-ball with friction dynamical system

It arises when Newton’s law is applied to a point subject to a constant friction $\gamma>0$ (of the velocity $\dot{x}(t)$) and a gravity potential $f$. This explains the naming “Heavy-ball method” and the interpretation of $\beta(x^{n}-x^{n-1})$ as inertial force.

Setting $\beta=0$ results in the forward-backward splitting algorithm, which has the nice property that in each iteration the function value decreases. Our convergence analysis reveals that the additional inertial term prevents our algorithm from monotonically decreasing the function values. Although this may look like a limitation on first glance, demanding monotonically decreasing function values anyway is too strict as it does not allow for provably optimal schemes. We refer to a statement of Nesterov : “In convex optimization the optimal methods never rely on relaxation. Firstly, for some problem classes this property is too expensive. Secondly, the schemes and efficiency estimates of optimal methods are derived from some global topological properties of convex functions”Relaxation is to be interpreted as the property of monotonically decreasing function values in this context. Topological properties should be associated with geometrical properties.. The negative side of better efficiency estimates of an algorithm is usually the convergence analysis. This is even true for convex functions. In case of non-convex and non-smooth functions, this problem becomes even more severe.

Despite this problem, we can establish convergence of the sequence of function values for the general case, where the objective function is only required to be a composition of a convex and a differentiable function. Regarding the sequence of arguments generated by the algorithm, existence of a converging subsequence is shown. Furthermore, we show that each limit point is a critical point of the objective function.

To establish convergence of the whole sequence in the non-convex case is very hard. However, with slightly more assumptions to the objective, namely that it satisfies the Kurdyka-Łojasiewicz inequality , several algorithms have been shown to converge . In an abstract convergence theorem for descent methods with certain properties is proved. It applies to many algorithms. However, it can not be used for our algorithm. Based on their analysis, we prove an abstract convergence theorem for a different class of descent methods, which applies to iPiano. By verifying the requirements of this abstract convergence theorem, we manage to also show such a strong convergence result. From a practical point of view of image processing, computer vision, or machine learning, the Kurdyka-Łojasiewicz inequality is almost always satisfied. For more details about properties of Kurdyka-Łojasiewicz functions and a taxonomy of functions that have this property, we refer to .

The last part of the paper is devoted to experiments. We exemplarily present results on computer vision tasks, such as denoising and image compression, and show that entering the staggering world of non-convex functions pays off in practice.

Related Work

In convex optimization, splitting algorithms usually originate from the proximal point algorithm . It is a very general algorithm, and results on its convergence affect many other algorithms. Practically, however, computing one iteration of the algorithm can be as hard as the original problem. Among the strategies to tackle this problem are splitting approaches like Douglas-Rachford , several primal-dual algorithms , and forward-backward splitting ; see for a survey.

Especially the forward-backward splitting schemes seem to be appealing to generalize to non-convex problems. This is due to their simplicity and the existence of simpler formulations in some special cases like, for example, the gradient projection method, where the backward-step is the projection onto a set . In the classical forward-backward algorithm, where the backward step is the solution of a proximal term involving a convex function, is studied for a non-convex problem. In fact, the same class of objective functions as in the present paper is analyzed. The algorithm presented here comprises the algorithm from as a special case. Also Nesterov briefly accounts this algorithm in a general setting. Even the reverse setting is generalized in the non-convex setting , namely where the backward-step is performed on a non-smooth non-convex function.

As the amount of data to be processed is growing and algorithms are supposed to exploit all the data in each iteration, inexact methods become interesting, though we do not consider erroneous estimates in this paper. Forward-backward splitting schemes also seem to work for non-convex problems with erroneous estimates . A mathematical analysis of inexact methods can be found, e.g., in , but with the restriction that the method is explicitly required to decrease the function values in each iteration. The restriction comes with significantly improved results with regard of the convergence of the algorithm. The algorithm proposed in this paper provides strong convergence results, although it does not require the function values to decrease.

In his seminal work , Polyak investigates multi-step schemes to accelerate the gradient method. It turns out that a particularly interesting case is given by a two-step algorithm, which has been coined the Heavy-ball method. The name of the method is because it can be interpreted as an explicit finite differences discretization of the so-called Heavy-ball with friction dynamical system. It differs from the usual gradient method by adding an inertial term that is computed by the difference of the two preceding iterations. Polyak showed that this method can speed up convergence in comparison to the standard gradient method, while the cost of each iteration stays basically unchanged.

The popular accelerated gradient method of Nesterov obviously shares some similarities with the Heavy-ball method, but it differs from it in one regard: while the Heavy-ball method uses gradients based on the current iterate, Nesterov’s accelerated gradient method evaluates the gradient at points that are extrapolated by the inertial force. On strongly convex functions, both methods are equally fast (up to constants), but Nesterov’s accelerated gradient method converges much faster on weakly convex functions .

The Heavy-ball method requires knowledge about the function parameters (Lipschitz constant of the gradient and the modulus of strong convexity) to achieve the optimal convergence rate, which can be seen as a disadvantage. Interestingly, the conjugate gradient method for minimizing strictly convex quadratic problems can be expressed as Heavy-ball method. Hence, it can be seen as a special case of the Heavy-ball method for quadratic problems. In this special case, no additional knowledge is required about the function parameters, as the algorithm parameters are computed online.

The Heavy-ball method was originally proposed for minimizing differentiable convex functions, but it has been generalized in different ways. In , it has been generalized to the case of smooth non-convex functions. It is shown that, by considering an appropriate Lyapunov objective function, the iterations are attracted by the connected components of stationary points. In Section 4 it will become evident that the non-convex Heavy-ball method is a special case of our algorithm, and also the convergence analysis of shows some similarities to ours.

In , the Heavy-ball method has been extended to maximal monotone operators, e.g., the subdifferential of a convex function. In a subsequent work , it has been applied to a forward-backward splitting algorithm, again in the general framework of maximal monotone operators.

An abstract convergence result

In order to give a sound description of the first order optimality condition for a non-convex non-smooth optimization problem, we have to introduce the generalization of the subdifferential for convex functions.

The limiting-subdifferential (or simply subdifferential) is defined by (see [40, Def. 8.3])

which makes use of the Fréchet subdifferential defined by

when $x\in\operatorname{dom}F$ and by $\widehat{\partial}F(x)=\varnothing$ else.

In what follows, we will consider the problem of finding a critical point $x^{*}\in\operatorname{dom}F$ of $F$, which is characterized by the necessary first-order optimality condition $0\in\partial F(x^{*})$.

We state the definition of the Kurdyka-Łojasiewicz property from .

If the function $F$ satisfies the Kurdyka-Łojasiewicz inequality at each point of $\operatorname{dom}\partial F$, it is called KL function.

Roughly speaking, this condition says that we can bound the subgradient of a function from below by a reparametrization of its function values. In the smooth case, we can also say that up to a reparametrization the function $h$ is sharp, meaning that any non-zero gradient can be bounded away from . This is sometimes called a desingularization. It has been shown in that a proper lower semi-continuous extended valued function $h$ always satisfies this inequality at each non-stationary point. For more details and other interpretations of this property, also for different formulations, we refer to .

A big class of functions that have the KL-property is given by real semi-algebraic functions . Real semi-algebraic functions are defined as functions whose graph is a real semi-algebraic set.

2 Inexact descent convergence result for KL functions

Based on these conditions, we derive the same convergence result as in . The statements and proofs of the subsequent results follow the same ideas as . We modified the involved calculations according to our conditions H1, H2, and H3.

These conditions are very similar to the ones in , however, they are not identical. The difference comes from the fact that does not consider a two-step algorithm.

In the corresponding condition to H1 (sufficient decrease condition) is $F(x^{n+1})+a\Delta_{n+1}^{2}\leq F(x^{n})$.

The corresponding condition to H2 (relative error condition) is $\|{w^{n+1}}\|_{2}\leq b\Delta_{n+1}$. In some sense, our condition H2 accepts a larger relative error.

Our proof and the proof in mainly differ in the calculations that are involved, the outline is the same. There is hope to find an even more general convergence result, which comprises ours and .

Moreover, the initial point $z^{0}=(x^{0},x^{-1})$ is such that $F(z^{*})\leq F(z^{0})<F(z^{*})+\eta$ and

The key points of the proof are the facts that for all $j\geq 1$:

As for $n\geq 1$ the set $\partial F(z^{n})$ is nonempty (see Condition H2) every $z^{n}$ belongs to $\operatorname{dom}F$. For notational convenience, we define

Now, we want to show that for $n\geq 1$ holds: if $F(z^{n})<F(z^{*})+\eta$ and $z^{n}\in B(z^{*},\rho)$, then

Obviously, we can assume that $\Delta_{n}\neq 0$ (otherwise it is trivial), and therefore H1 and (2) imply $F(z^{n})>F(z^{n+1})\geq F(z^{*})$. The KL inequality shows $w^{n}\neq 0$ and H2 shows $\Delta_{n}+\Delta_{n-1}>0$. Since $w^{n}\in\partial F(z^{n})$, using KL inequality and H2, we obtain

As $\varphi$ is concave and increasing ($\varphi^{\prime}>0$), Condition H1 and (2) yield

which by applying $2\sqrt{uv}\leq u+v$ establishes (7).

As (2) does only imply $z^{n+1}\in B(z^{*},\sigma)$, $\sigma>\rho$, we can not use (7) directly for the whole sequence. However, (5) and (6) can be shown by induction on $j$. For $j=0$, (2) yields $z^{1}\in B(z^{*},\sigma)$ and $F(z^{1}),F(z^{2})\geq F(z^{*})$. From Condition H1 with $n=1$, $F(z^{2})\geq F(z^{*})$ and $F(z^{1})\leq F(z^{0})$, we infer

and therefore $z^{1}\in B(z^{*},\rho)$. Direct use of (7) with $n=1$ shows that (6) holds with $j=1$.

Suppose (5) and (6) are satisfied for $j\geq 1$. Then, using the triangle inequality and (6), we have

which shows, using $\Delta_{j+1}\leq\sqrt{\frac{1}{a}(F(z^{j+1})-F(z^{j+2}))}\leq\sqrt{\frac{1}{a}(F(z^{0})-F(z^{*}))}$ and (3), that $z^{j+1}\in B(z^{*},\rho)$. As a consequence (7), with $n=j+1$, can be added to (6) and we can conclude (6) with $j+1$. This shows the desired induction on $j$.

Therefore, $x^{n}$ converges to some $\bar{x}$ as $n\to\infty$, and $z^{n}$ converges to $\bar{z}=(\bar{x},\bar{x})$. As $\varphi$ is concave, $\varphi^{\prime}$ is decreasing. Using this and Condition H2 yields $w^{n}\to 0$ and $F(z^{n})\to\zeta\geq F(z^{*})$. Suppose we have $\zeta>F(z^{*})$, then KL-inequality reads $\varphi^{\prime}(\zeta-F(z^{*}))\|{w^{n}}\|_{2}\geq 1$ for all $n\geq 1$, which contradicts $w^{n}\to 0$.

The next corollary and the subsequent theorem follow as in by replacing the calculation with our conditions.

By Condition H1, for $z^{n}\in B(z^{*},\rho)$, we have

Using the triangle inequality on $\|{z^{n+1}-z^{*}}\|_{2}$ shows that $z^{n+1}\in B(z^{*},\sigma)$, which implies (2) and concludes the proof. ∎

The work that is done in Lemma 5 and Corollary 6 allows us to formulate an abstract convergence theorem for sequences satisfying the Conditions H1, H2, and H3. It follows, with a few modifications, as in .

The proposed algorithm - iPiano

The proposed algorithm, which is stated in Subsection 4.3, seeks for a critical point $x^{*}\in\operatorname{dom}h$ of $h$, which is characterized by the necessary first-order optimality condition $0\in\partial h(x^{*})$. In our case, this is equivalent to

This equivalence is explicitly verified in the next subsection, where we collect some details and state some basic properties, which are used in the convergence analysis in Subsection 4.5.

2 Preliminaries

Consider the function $f$ first. It is required to be $C^{1}$-smooth with $L$-Lipschitz continuous gradient on $\operatorname{dom}g$, i.e., there exists a constant $L>0$ such that

This directly implies that $\operatorname{dom}h=\operatorname{dom}g$ is a non-empty convex set, as $\operatorname{dom}g\subset\operatorname{dom}f$. This property of $f$ plays a crucial role in our convergence analysis due to the following lemma (stated as in ).

We assume that the function $g$ is a proper lower semi-continuous convex function with an efficient to compute proximal map.

Let $g$ be a proper lower semi-continuous convex function. Then, we define the proximal map

An important (basic) property that the convex function $g$ contributes to the convergence analysis is the following:

Let $g$ be a proper lower semi-continuous convex function, then it holds for any $x,y\in\operatorname{dom}g$, $s\in\partial g(x)$ that

This result follows directly from the convexity of $g$. ∎

Finally, consider the optimality condition $0\in\partial h(x^{*})$ more in detail. The following proposition proves the equivalence to $-\nabla f(x^{*})\in\partial g(x^{*})$. The proof is mainly based on Definition 2 of the limiting-subdifferential.

Let $h$, $f$, and $g$ be like before, i.e., let $h=f+g$ with $f$ continuously differentiable and $g$ convex. Sometimes, $h$ is then called a $C^{1}$-perturbation of a convex function. Then, for $x\in\operatorname{dom}h$ holds

We first prove “$\subset$”. Let $\xi^{h}\in\partial h(x)$, i.e., there is a sequence $(y_{k})_{k=0}^{\infty}$ such that $y_{k}\to x$, $h(y_{k})\to h(x)$, and $\xi_{k}^{h}\to\xi^{h}$, where $\xi_{k}^{h}\in\widehat{\partial}h(y_{k})$. We want to show that $\xi^{g}:=\xi^{h}-\nabla f(x)\in\partial g(x)$. As $f\in C^{1}$ and $\xi^{h}\in\partial h(x)$, we have

where $\liminf$ and $\lim$ are over $y_{k}^{\prime}\to y_{k},y_{k}^{\prime}\neq y_{k}$. Therefore, $\xi^{g}_{k}\in\widehat{\partial}g(y_{k})$. The other inclusion “$\supset$” is trivial. ∎

As a consequence, a critical point can also be characterized by the following definition.

Let $f$ and $g$ be as afore. Then, we define the proximal residual

It can be easily seen that $r(x)=0$ is equivalent to $x=(I+\partial g)^{-1}(x-\nabla f(x))$ and $(I+\partial g)(x)=(I-\nabla f)(x)$, which is the first-order optimality condition. The proximal residual is defined with respect to a fixed step size of $1$. The rationale behind this becomes obvious when $g$ is the indicator function of a convex set. In this case, a small residual could be caused by small step sizes as the reprojection onto the convex set is independent of the step size.

3 The generic algorithm

In this paper, we propose an algorithm, iPiano, with the generic formulation in Algorithm 1. It is a forward-backward splitting algorithm incorporating an inertial force. In the forward step, $\alpha_{n}$ determines the step size in the direction of the gradient of the differentiable function $f$. The step in gradient direction is aggregated with the inertial force from the previous iteration weighted by $\beta_{n}$. Then, the backward step is the solution of the proximity operator for the function $g$ with the weight $\alpha_{n}$.

In order to make the algorithm specific and convergent, the step size parameters must be chosen appropriately. What “appropriately” means, will be specified in Subsection 4.4 and proved in Subsection 4.5.

4 Rules for choosing the step size

In this subsection, we propose several strategies for choosing the step sizes. This will make it easier to implement the algorithm. One may choose among the following variants of step size rules depending on the knowledge about the objective function.

The most simple one, which requires most knowledge about the objective function, is outlined in Algorithm 2. All step size parameters are chosen a priori and are constant.

Observe that our law on $\alpha,\beta$ is equivalent to the law found in for minimizing a smooth non-convex function. Hence, our result can be seen as an extension of their work to the presence of an additional non-smooth convex function.

The case where we have only limited knowledge about the objective function occurs more frequently. It can be very challenging to estimate the Lipschitz constant of $\nabla f$ beforehand. Using backtracking the Lipschitz constant can be estimated automatically. A sufficient condition that the Lipschitz constant at iteration $n$ to $n+1$ must satisfy is

Although, there are different strategies to determine $L_{n}$, the most common one is by defining an increment variable $\eta>1$ and looking for $L_{n}\in\{L_{n-1},\eta L_{n-1},\eta^{2}L_{n-1},\ldots\}$ minimal satisfying (15). Sometimes, it is also feasible to decrease the estimated Lipschitz constant after a few iterations. A possible strategy is as follows: if $L_{n}=L_{n-1}$, then search for the minimal $L_{n}\in\{\eta^{-1}L_{n-1},\eta^{-2}L_{n-1},\ldots\}$ satisfying (15).

In Algorithm 3 we propose an algorithm with variable step sizes. Any strategy for estimating the Lipschitz constant may be used. When changing the Lipschitz constant from one iteration to another, all step size parameters must be adapted. The rules for adapting the step sizes will be justified during the convergence analysis in Subsection 4.5.

Algorithm 5 defines the general rules that the step size parameters must satisfy.

It contains the Algorithms 2, 3, and 4 as special instances. This is easily verified for Algorithms 2 and 4. For Algorithm 3 the step size rules are derived from the proof of Lemma 13.

As Algorithm 5 is the most general one, now, let us analyze the behavior of this algorithm.

5 Convergence analysis

Let us first verify that the algorithm makes sense. We have to show that the requirements to the parameters are not contradictory, i.e., that it is possible to choose a feasible set of parameters. In the following Lemma, we will only show existence of such a parameter set, however, the proof helps us to formulate specific step size rules.

For all $n\geq 0$, there are $\delta_{n}\geq\gamma_{n}$, $\beta_{n}\in[0,1)$, and $\alpha_{n}<{2(1-\beta_{n})}/{L_{n}}$. Furthermore, given $L_{n}>0$, there exists a choice of parameter $\alpha_{n}$ and $\beta_{n}$ such that additionally $(\delta_{n})_{n=0}^{\infty}$ is monotonically decreasing.

The upper bound for $\beta_{n}$ and $\alpha_{n}$ come from rearranging $\gamma_{n}\geq c_{2}$ to $\beta_{n}\leq 1-\alpha_{n}L_{n}/2-c_{2}\alpha_{n}$ and $\alpha_{n}\leq 2(1-\beta_{n})/(L_{n}+2c_{2})$, respectively. The last statement follows by incorporating the descent property of $\delta_{n}$. Let $\delta_{-1}\geq c_{2}$ be chosen initially. Then, the decent property of $(\delta_{n})_{n=0}^{\infty}$ requires one of the equivalent statements

to be true. An upper bound on $\alpha_{n}$ is obtained by

The only thing that remains to show is that there exists $\alpha_{n}>c_{1}$ and $\beta_{n}\in[0,1)$ such that these two relations are fulfilled. Consider the condition for a non-negative gap between the upper and lower bound for $\alpha_{n}$

Defining $b:=(\delta_{n-1}+\frac{L_{n}}{2})/(c_{2}+\frac{L_{n}}{2})\geq 1$, it is easily verified that there exists $\beta_{n}\in[0,1)$ satisfying the equivalent condition

As a consequence, the existence of a feasible $\alpha_{n}$ follows, and the decent property for $\delta_{n}$ holds. ∎

In the following proposition, we state a result which will be very useful. Although, iPiano does not imply a descent property of the function values, we construct a majorizing function that enjoys a monotonically descent property. This function reveals the connection to the Lyapunov direct method for convergence analysis as used in .

The sequence $(H_{\delta_{n}}(x^{n},x^{n-1}))_{n=0}^{\infty}$ is monotonically decreasing and thus converging. In particular, it holds

It holds $\sum_{n=0}^{\infty}\Delta_{n}^{2}<\infty$ and, thus, $\lim_{n\to\infty}\Delta_{n}=0$.

Now using $x=x^{n+1}$ and $y=x^{n}$ in (11) and (12) and summing both inequalities it follows that

which establishes (23) as $\delta_{n}$ is monotonically decreasing. Obviously, the sequence $(H_{\delta_{n}}(x^{n},x^{n-1}))_{n=0}^{\infty}$ is monotonically decreasing if and only if $\gamma_{n}\geq 0$, which is true by the algorithm’s requirements. By assumption, $h$ is bounded from below by some constant $\underline{h}>-\infty$, hence $(H_{\delta_{n}}(x^{n},x^{n-1}))_{n=0}^{\infty}$ converges.

Summing up (23) from $n=0,\ldots,N$ yields (note that $H_{\delta_{n}}(x^{0},x^{-1})=h(x^{0})$)

Letting $N$ tend to $\infty$ and remembering that $\gamma_{N}\geq c_{2}>0$ holds implies the statement.

The function $H_{\delta}$ is a Lyapunov function for the dynamical system of described by the Heavy-ball method. It corresponds to a discretized version of the kinetic energy of the Heavy-ball with friction.

In the following theorem, we state our general convergence results about Algorithm 5.

The sequence $(h(x^{n}))_{n=0}^{\infty}$ converges.

There exists a converging subsequence $(x^{n_{k}})_{k=0}^{\infty}$.

Any limit point ${x^{*}:=\lim_{k\to\infty}x^{n_{k}}}$ is a critical point of (9) and $h(x^{n_{k}})\to h(x^{*})$ as $k\to\infty$.

This follows from the Squeeze theorem as for all $n\geq 0$ holds

and thanks to Proposition 14(a) and (b) holds

To show that each limit point $x^{*}:=\lim_{j\to\infty}x^{n_{j}}$ is a critical point of (9) recall that the subdifferential (1) is closed . Define

it holds $\lim_{j\to\infty}\xi^{j}=0$. It remains to show that $\lim_{j\to\infty}h(x^{n_{j}})=h(x^{*})$. By the closure property of the subdifferential $\partial h$ it is $(x^{*},0)\in\operatorname{Graph}(\partial h)$, which means that $x^{*}$ is a critical point of $h$.

The continuity statement about the limiting process as $j\to\infty$ follows by the lower semi-continuity of $g$, the existence $\lim_{j\to\infty}\xi^{j}=0$, and the convexity property in Lemma 10

The first equality holds because the subadditivity of $\limsup$ becomes an equality when the limit exists for one of the two summed sequencesIn general, the existence of $(\xi^{j})_{j=0}^{\infty}$ is not guaranteed. Compared to the general case, additionally $\lim_{j\to\infty}\xi^{j}=0$ is known here., here it exists $\lim_{j\to\infty}\left\langle\xi^{j},x^{*}-x^{n_{j}}\right\rangle=0$. Moreover, as $f$ is differentiable it is also continuous, thus $\lim_{j\to\infty}f(x^{n_{j}})=f(x^{*})$. This implies $\lim_{j\to\infty}h(x^{n_{j}})=h(x^{*})$.

The convergence properties shown in Theorem 15 should be the basic requirement of any algorithm. Very loosely speaking, it states that the algorithm ends up in a meaningful solution. It allows us to formulate stopping conditions, e.g., the residual between successive function values.

Condition H1 is proved in Proposition 14(a) with $a=c_{2}\leq\gamma_{n}$.

To proof Condition H2, consider $w^{n+1}:=(w_{x}^{n+1},w_{y}^{n+1})^{\top}\in\partial H_{\delta}(x^{n+1},x^{n})$ with $w_{x}^{n+1}\in\partial g(x^{n+1})+\nabla f(x^{n+1})+2\delta(x^{n+1}-x^{n})$ and $w_{y}^{n+1}=-2\delta(x^{n+1}-x^{n})$. The Lipschitz continuity of $\nabla f$ and using (19) to specify an element from $\partial g(x^{n+1})$ imply

As $\alpha_{n}L_{n}\leq 2(1-\beta_{n})\leq 2$ and $\delta\alpha_{n}=1-\frac{1}{2}\alpha_{n}L_{n}-\frac{1}{2}\beta_{n}\leq 1$, setting $b=\frac{7}{c_{1}}$ verifies condition H2, i.e., $\|{w^{n+1}}\|_{2}\leq b(\Delta_{n}+\Delta_{n+1})$.

Now, the abstract convergence Theorem 7 concludes the proof. ∎

The next corollary makes use of the fact that semi-algebraic functions (Definition 4) have the Kurdyka-Łojasiewicz property.

As $h$ and $\delta\|{x-y}\|_{2}$ are semi-algebraic, $H_{\delta}(x,y)$ is semi-algebraic and has the KL property. Then, Theorem 16 concludes the proof. ∎

6 Convergence rate

In the following, we are interested in determining a convergence rate with respect to the proximal residual from Definition 12. Since all preceding estimations are according to $\|{x^{n+1}-x^{n}}\|_{2}$ we establish the relation to $\|{r(x)}\|_{2}$ first. The following lemmas about the monotonicity and the non-expansiveness of the proximity operator turn out to be very useful for that. Coarsely speaking, Lemma 18 states that the residual is sub-linearly increasing. Lemma 19 formulates a standard property of the proximal operator.

Then, $p_{g}(\alpha)$ is a decreasing function of $\alpha$, and $q_{g}(\alpha)$ increasing in $\alpha$.

See e.g. [36, Lemma 1] or [44, Lemma 4]. ∎

Let $g$ be a convex function and $\alpha>0$, then, for all $x,y\in\operatorname{dom}g$ we obtain the non-expansiveness of the proximity operator

It is a well-known fact. See for example . ∎

The two preceding lemmas allow us to establish the following relation.

First, we observe the relations $1\leq\alpha\Rightarrow q_{g}(1)\leq q_{g}(\alpha)$ and $1\geq\alpha\Rightarrow p_{g}(1)\leq p_{g}(\alpha)=\frac{1}{\alpha}q_{g}(\alpha)$, which are based on Lemma 18. Then, invoking the non-expansiveness of the proximity operator (Lemma 19) we obtain

This allows us to compute the following lower bound

where the first inequality arises from adding zero and using (26), the second uses the triangle inequality, the next one applies Lemma 18 and $\beta_{n}<1$. Now, summing both sides from $n=0,\ldots,N$ and using $x^{-1}=x^{0}$ the statement easily follows. ∎

Algorithm 5 guarantees that for all $N\geq 0$

In view of Proposition 14(a), and the definition of $\gamma_{N}$ in (21), summing up both sides of (23) for $n=0,\ldots,N$ and using that $\delta_{N}>0$ from (21) we obtain

As it is $\gamma_{n}>c_{2}$, a simple rearrangement invoking Lemma 20 concludes the proof. ∎

A similar result can be found in for the case $\beta=0$.

Numerical experiments

Let us present some of the qualitative properties of the proposed algorithm. For this, we consider to minimize the following simple problem

where $x$ is the unknown vector, $u^{0}$ is some given vector, and $\lambda,\mu>0$ are some free parameters. A contour plot and the energy landscape of $h$ in the case of $N=2$, $\lambda=1$, $\mu=100$, and $u^{0}=(1,1)^{\top}$ is depicted in Figure 1. It turns out that the function $h$ has four stationary points, i.e. points $\bar{x}$, such that $0\in\nabla f(\bar{x})+\partial g(\bar{x})$. These points are marked by small black diamonds.

Clearly the function $f$ is non-convex but has a Lipschitz continuous gradient with components

The Lipschitz constant of $\nabla f$ is easily computed as $L=\mu$. The function $g$ is non-smooth but convex and the proximal operator with respect to $g$ is given by the well-known shrinkage operator

where all operations are understood component-wise. Let us test the performance of the proposed algorithm on the example shown in Figure 1. We set $\alpha=2(1-\beta)/L$. Figure 2 shows the results of using the iPiano algorithm for different settings of the extrapolation factor $\beta$. We observe that iPiano with $\beta=0$ is strongly attracted by the closest stationary points while switching on the inertial term can help to overcome the spurious stationary points. The reason for this desired property is that while the gradient might vanish at some points, the inertial term $\beta(x^{n}-x^{n-1})$ is still strong enough to drive the sequence out of the stationary region. Clearly, there is no guarantee that iPiano always avoids spurious stationary points. iPiano is not designed to find the global optimum. However, our numerical experiments suggest that in many cases, iPiano finds lower energies than the respective algorithm without inertial term. A similar observation about the Heavy-ball method is described in .

2 Image processing applications

It is well-known that non-convex regularizers are better models for many image processing and computer vision problems, see e.g. . However, convex models are still preferred over non-convex ones, since they can be efficiently optimized using convex optimization algorithms. In this section, we demonstrate the applicability of the proposed algorithm to solve a class of non-convex regularized variational models. We present examples for natural image denoising, and linear diffusion based image compression. We show that iPiano can be easily adapted to all these problems and yields state-of-the-art results.

In this subsection, we investigate the task of natural image denoising. For this we exploit an optimized MRF (Markov random field) model, which is learned in following , and make use of the iPiano algorithm to solve it. In order to evaluate the performance of iPiano, we compare it to the well-known bound constrained limited memory quasi Newton method (L-BFGS) We make use of the implementation distributed at http://www.cs.toronto.edu/~liam/software.shtml.. As an error measure, we use the energy difference

where $h^{n}$ is the energy of the current iteration $n$ and $h^{*}$ is the energy of the true solution. Clearly, this error measure makes sense only when different algorithms can achieve the same true energy $h^{*}$ which is in general wrong for non-convex problems. In our image denoising experiments, however, we find, that all tested algorithms find the same solution, independent of the initialization. This can be explained by the fact that the learning procedure also delivers models that are relatively easy to optimize, since otherwise they would have resulted in a bad training error. In order to compute a true energy $h^{*}$, we run the iPiano algorithm with a proper $\beta$ (e.g., $\beta=0.8$) for enough iterations ($\sim$1000 iterations). We run all the experiments in Matlab on a 64-bit Linux server with 2.53GHz CPUs.

The MRF image denoising model based on learned filters is formulated as

and for the impulse noise (e.g., salt & pepper noise), $g_{1,2}$ is given as

The parameter $\lambda>0$ is used to define the tradeoff between regularization and data fitting.

In this paper, we consider the following non-convex penalty function, which is derived from the Student-t distribution:

Let us now explain how to solve (30) using the iPiano algorithm. Casting (30) in the form of (9), we see that $f(u)=\sum_{i=1}^{N_{f}}\vartheta_{i}\Phi(K_{i}u)$ and $g(u)=g_{1,2}(u,u^{0})$. Thus, we have

where $\Phi^{\prime}(K_{i}u)=[\varphi^{\prime}((K_{i}u)_{1})\,,\varphi^{\prime}((K_{i}u)_{2}),\dots,\varphi^{\prime}((K_{i}u)_{p})]^{\top}$ and $\varphi^{\prime}(t)={2t}/{(1+t^{2})}$. The proximal map with respect to $g$ simply poses point-wise operations. For the case of $g_{2}$, it is given by

where $\beta$ is a free parameter to be evaluated in the experiment. In order to make use of possible larger step sizes in practice, we use a following trick: when the inequality (15) is fulfilled, we decrease the evaluated Lipschitz constant $L_{n}$ slightly by setting $L_{n}=L_{n}/1.05$.

The iPiano algorithm has an additional advantage of simplicity. The iPiano version without backtracking basically relies on matrix vector products (filter operations in the denoising examples) and simple pointwise operations. Therefore, the iPiano algorithm is well suited for a parallel implementation on GPUs which an lead to speedup factors of 20-30.

2.2 Linear diffusion based image compression

In this example we apply the iPiano algorithm to linear diffusion based image compression. Recent works have shown that image compression based on linear and non-linear diffusion can outperform the standard JPEG standard and even the more advanced JPEG 2000 standard, when the interpolation points are carefully chosen. Therefore, finding optimal data for interpolation is a key problem in the context of PDE-based image compression. There exist only few prior works for this topic, see e.g. , and the very recent approach presented in defines the state-of-the-art.

The problem of finding optimal data for homogeneous diffusion-based interpolation is formulated as the following constrained minimization problem:

Observe that if $c\in[0,1)^{N}$, we can multiply the constraint equation in (33) from the left by $(I-C)^{-1}$ such that it becomes

where $E(c)=\mathop{\operator@font diag}\nolimits(c_{1}/(1-c_{1}),...,c_{N}/(1-c_{N}))$. This shows that problem (33) is in fact a reduced formulation of the bilevel optimization problem

where $D$ is the nabla operator and hence $-L=D^{\top}D$.

Problem (33) is non-convex due to the non-convexity of the equality constraint. In , the above problem is solved by a successive primal-dual (SPD) algorithm, which successively linearizes the non-convex constraint and solves the resulting convex problem with the first-order primal-dual algorithm . The main drawback of SPD is, that it requires tens of thousands inner iterations and thousands of outer iterations to reach a reasonable solution. However, as we now demonstrate, iPiano can solve this problem with higher accuracy in 1000 iterations.

Observe that we can rewrite the problem (33) by solving $u$ from the constraints equation, which gives

where $A=C+(C-I)L$. In , it is shown that the $A$ is invertible as long as at least one element of $c$ is non-zero, which is the case for non-degenerate problems. Substituting back the above equation into (33), we arrive at the following optimization problem, which now only depends on the inpainting mask $c$:

Casting (35) in the form of (9), we have $f(c)=\frac{1}{2}\|A^{-1}Cu^{0}-u^{0}\|_{2}^{2}$, and $g(c)=\lambda\|c\|_{1}$. In order to minimize the above problem using iPiano, we need to calculate the gradient of $f$ with respect to $c$. This is shown by the following lemma.

By substituting (38) into (37), we obtain

Finally, we need to compute the proximal map with respect to $g(c)$ which is again given by a pointwise application of the shrinkage operator (28).

Now, we can make use of the iPiano algorithm to solve the problem (35). We set $\beta=0.8$, which generally performs very well in practice. We additionally accelerate the SPD algorithm used in the previous work by applying the diagonal preconditioning technique , which significantly reduces the required iterations for the primal-dual algorithm in the inner loop.

Figure 7 shows examples of finding optimal interpolation data for the three test images. Table 3 summarizes the results of two different algorithms. Regarding the reconstruction quality, we make use of the mean squared error (MSE) as an error measurement to keep consistent with previous work, which is computed by

From Table 3, one can see that the Successive PD algorithm requires $200\times 4000$ iterations to converge. iPiano only needs 1000 iterations to reach already a lower energy. Note that in each iteration of the iPiano algorithm, two linear systems have to be solved. In our implementation we use the Matlab “backslash” operator which effectively exploits the strong sparseness of the systems. A lower energy basically implies that iPiano can solve the minimization problem (33) better. Regarding the final compression result, usually the result of iPiano has slightly less density, but slightly worse MSE. Following the work , we also consider the so-called gray value optimization (GVO) as a post-processing step to further improve the MSE of the reconstructed images.

Conclusions

In this paper, we have proposed a new optimization algorithm, which we call iPiano. It is applicable to a broad class of non-convex problems. More specifically, it addresses objective functions, which are composed as a sum of a differentiable (possibly non-convex) and a convex (possibly non-differentiable) function. The basic methodologies have been derived from the forward-backward splitting algorithm and the Heavy-ball method.

Our theoretical convergence analysis is divided into two steps. First, we have proved an abstract convergence result about inexact descent methods. Then, we analyze the convergence of iPiano. For iPiano, we have proved that the sequence of function values converges, that the subsequence of arguments generated by the algorithm is bounded, and that every limit point is a critical point of the problem. Requiring the Kurdyka-Łojasiewicz property for the objective function establishes deeper insights into the convergence behavior of the algorithm. Using the abstract convergence result, we have shown that the whole sequence converges and the unique limit point is a stationary point.

The analysis includes an examination of the convergence rate. A rough upper bound of $O(1/n)$ has been found for the squared proximal residual. Experimentally, iPiano has been shown to have a much faster convergence rate.

Finally, the applicability of the algorithm has been demonstrated and iPiano achieved state-of-the-art performance. The experiments comprised denoising and image compression. In the first two experiments, iPiano helped learning a good prior for the problem. In the case of image compression, iPiano has demonstrated its use in a huge optimization problem for computing an optimal mask for a Laplacian PDE-based image compression method.

In summary, iPiano has many favorable theoretical properties, is simple and efficient. Hence, we recommend it as a standard solver for the considered class of problems.

Acknowledgements

We are grateful to Joachim Weickert for discussions about the image compression by diffusion problem.

References