A forward-backward view of some primal-dual optimization methods in image recovery

Introduction

Many image recovery problems can be formulated in Hilbert spaces $\mathcal{H}$ and $(\mathcal{G}_{i})_{1\leqslant i\leqslant m}$ as structured optimization problems of the form

where, for every $i\in\{1,\ldots,m\}$, $g_{i}$ is a proper lower semicontinuous convex function from $\mathcal{G}_{i}$ to $\,\left]-\infty,+\infty\right]$ and $L_{i}$ is a bounded linear operator from $\mathcal{H}$ to $\mathcal{G}_{i}$. For example, the functions $(g_{i}\circ L_{i})_{1\leqslant i\leqslant m}$ may model data fidelity terms, smooth or nonsmooth measures of regularity, or hard constraints on the solution. In recent years, many algorithms have been developed to solve such a problem by taking advantage of recent advances in convex optimization, especially in the development of proximal tools (see and the references therein). In image processing, however, solving such a problem still poses a number of conceptual and numerical challenges. First of all, one often looks for methods which have the ability to split the problem by activating each of the functions through elementary processing steps which can be computed in parallel. This makes it possible to reduce the complexity of the original problem and to benefit from existing parallel computing architectures. Secondly, it is often useful to design algorithms which can exploit, in a flexible manner, the structure of the problem. In particular, some of the functions may be Lipschitz differentiable in which case they should be exploited through their gradient rather than through their proximity operator, which is usually harder to implement (examples of proximity operators with closed-form expression can be found in ). In some problems, the functions $(g_{i})_{1\leqslant i\leqslant m}$ can be expressed as the infimal convolution of simpler functions (see and the references therein). Last but not least, in image recovery, the operators $(L_{i})_{1\leqslant i\leqslant m}$ may be of very large size so that their inversions are costly (e.g., in reconstruction problems). Finding algorithms which do not require to perform inversions of these operators is thus of paramount importance.

Note that all the existing convex optimization algorithms do not have these desirable properties. For example, the Alternating Direction Method of Multipliers (ADMM) requires a stringent assumption of invertibility of the involved linear operator. Parallel versions of ADMM and related Parallel Proximal Algorithm (PPXA) usually necessitate a linear inversion to be performed at each iteration. Also, early primal-dual algorithms did not make it possible to handle smooth functions through their gradients. Only recently, have primal-dual methods been proposed with this feature. Such work was initiated in in the line of and subsequent developments can be found in . As will be seen in the present paper, another advantage of these approaches is that they can be coupled with variable metric strategies which can potentially accelerate their convergence.

In Section 2, we provide some background on convex analysis and monotone operator theory. In Section 3, we introduce a general form of the forward-backward algorithm which uses a variable metric. This algorithm is employed in Section 4 to develop a versatile family of primal-dual proximal methods. Several particular instances of this framework are discussed. Finally, we provide illustrating numerical results in Section 5.

Notation and background

Monotone operator theory provides a both insightful and elegant framework for dealing with convex optimization problems and developing new solution algorithms that could not be devised using purely variational tools. We summarize a number of related concepts that will be needed.

Throughout, $\mathcal{H}$, $\mathcal{G}$, and $(\mathcal{G}_{i})_{1\leqslant i\leqslant m}$ are real Hilbert spaces. We denote the scalar product of a Hilbert space by $\left\langle{\cdot}\mid{\cdot}\right\rangle$ and the associated norm by $\|\cdot\|$. The symbol $\rightharpoonup$ denotes weak convergence, In a finite dimensional space, weak convergence is equivalent to strong convergence. and $\operatorname{Id}$ denotes the identity operator. We denote by $\mathcal{B}\,(\mathcal{H},\mathcal{G})$ the space of bounded linear operators from $\mathcal{H}$ to $\mathcal{G}$, we set $\mathcal{S}\,(\mathcal{H})=\big\{{L\in\mathcal{B}\,(\mathcal{H},\mathcal{H})}\mid{L=L^{*}}\big\}$, where $L^{*}$ denotes the adjoint of $L$. The Loewner partial ordering on $\mathcal{S}\,(\mathcal{H})$ is denoted by $\succcurlyeq$. For every $\alpha\in\left[0,+\infty\right[$, we set $\mathcal{P}_{\alpha}(\mathcal{H})=\big\{{U\in\mathcal{S}\,(\mathcal{H})}\mid{U\succcurlyeq\alpha\operatorname{Id}}\big\},$ and we denote by $\sqrt{U}$ the square root of $U\in\mathcal{P}_{\alpha}(\mathcal{H})$. Moreover, for every $U\in\mathcal{P}_{\alpha}(\mathcal{H})$ and $\alpha>0$, we define the norm $\|x\|_{U}=\sqrt{\left\langle{Ux}\mid{x}\right\rangle}$.

We denote by $\boldsymbol{\mathcal{G}}=\mathcal{G}_{1}\oplus\cdots\oplus\mathcal{G}_{m}$ the Hilbert direct sum of the Hilbert spaces $(\mathcal{G}_{i})_{1\leqslant i\leqslant m}$, i.e., their product space equipped with the scalar product $\colon(\boldsymbol{x},\boldsymbol{y})\mapsto\sum_{i=1}^{m}\left\langle{x_{i}}\mid{y_{i}}\right\rangle$ where ${\boldsymbol{x}}=(x_{i})_{1\leqslant i\leqslant m}$ and ${\boldsymbol{y}}=(y_{i})_{1\leqslant i\leqslant m}$ denote generic elements in $\boldsymbol{\mathcal{G}}$.

Let $A\colon\mathcal{H}\to 2^{\mathcal{H}}$ be a set-valued operator. We denote by $\operatorname{gra}A=\big\{{(x,u)\in\mathcal{H}\times\mathcal{H}}\mid{u\in Ax}\big\}$ the graph of $A$, by $\operatorname{zer}A=\big\{{x\in\mathcal{H}}\mid{0\in Ax}\big\}$ the set of zeros of $A$, and by $\operatorname{ran}A=\big\{{u\in\mathcal{H}}\mid{(\exists\;x\in\mathcal{H})\;u\in Ax}\big\}$ its range. The inverse of $A$ is $A^{-1}\colon\mathcal{H}\mapsto 2^{\mathcal{H}}\colon u\mapsto\big\{{x\in\mathcal{H}}\mid{u\in Ax}\big\}$, and the resolvent of $A$ is $J_{A}=(\operatorname{Id}+A)^{-1}$. Moreover, $A$ is monotone if

and maximally monotone if it is monotone and there exists no monotone operator $B\colon\mathcal{H}\to 2^{\mathcal{H}}$ such that $\operatorname{gra}A\subset\operatorname{gra}B$ and $A\neq B$. An operator $B\colon\mathcal{H}\to\mathcal{H}$ is $\beta$-cocoercive for some $\beta\in\,\left]0,+\infty\right[$ if

The conjugate of a function $f\colon\mathcal{H}\to\,\left]-\infty,+\infty\right]$ is

and the infimal convolution of $f$ with $g\colon\mathcal{H}\to\,\left]-\infty,+\infty\right]$ is

The class of lower semicontinuous convex functions $f\colon\mathcal{H}\to\,\left]-\infty,+\infty\right]$ such that $\operatorname{dom}f=\big\{{x\in\mathcal{H}}\mid{f(x)<+\infty}\big\}\neq{\varnothing}$ is denoted by $\Gamma_{0}(\mathcal{H})$. If $f\in\Gamma_{0}(\mathcal{H})$, then $f^{*}\in\Gamma_{0}(\mathcal{H})$ and the subdifferential of $f$ is the maximally monotone operator

Let $U\in\mathcal{P}_{\alpha}(\mathcal{H})$ for some $\alpha\in\,\left]0,+\infty\right[$. The proximity operator of $f\in\Gamma_{0}(\mathcal{H})$ relative to the metric induced by $U$ is [22, Section XV.4]

When $U=\operatorname{Id}$, we retrieve the standard definition of the proximity operator . Let $C$ be a nonempty subset of $\mathcal{H}$. The indicator function of $C$ is defined on $\mathcal{H}$ as

A general form of Forward-Backward algorithm

Optimization problems can often be reduced to finding a zero of a sum of two maximally monotone operators $A$ and $B$ acting on $\mathcal{H}$. When $B$ is cocoercive (see (3)), a useful algorithm to solve this problem is the forward-backward algorithm, which can be formulated in a general form involving a variable metric as shown in the next result.

Then $x_{n}\rightharpoonup\overline{x}$ for some $\overline{x}\in Z$.

A variable metric primal-dual method

A wide array of optimization problems encountered in image processing are instances of the following one, which was first investigated in and can be viewed as a more structured version of the minimization problem in (1):

Let us now introduce the product space $\boldsymbol{\mathcal{K}}=\mathcal{H}\oplus\boldsymbol{\mathcal{G}}$ and the operators

The operator ${\boldsymbol{A}}$ can be shown to be maximally monotone,whereas ${\boldsymbol{B}}$ is cocoercive. A key observation in this context is that, if there exists $(\overline{x},\overline{\boldsymbol{v}})\in\boldsymbol{\mathcal{K}}$ such that $(\overline{x},\overline{\boldsymbol{v}})\in\operatorname{zer}({\boldsymbol{A}}+{\boldsymbol{B}})$, then $(\overline{x},\overline{\boldsymbol{v}})$ is a pair of primal-dual solutions to Problem 4.1 . This connection with the construction for a zero of ${\boldsymbol{A}}+{\boldsymbol{B}}$ makes it possible to apply a forward-backward algorithm as discussed in Section 3, by using a linear operator $\boldsymbol{V}_{n}\in\mathcal{B}\,(\boldsymbol{\mathcal{K}},\boldsymbol{\mathcal{K}})$ to change the metric at each iteration $n$. Depending on the form of this operator various algorithms can be obtained.

2 A first class of primal-dual algorithms

The following result constitutes a direct extension of [14, Example 6.4]:

3 A second class of primal-dual algorithms

where $\widetilde{\boldsymbol{U}}_{n}$ is given by (18).

The following result can then be deduced from Theorem 3.1. Its proof is skipped due to the lack of space.

Application to image restoration

An estimate $\widetilde{x}$ of $\overline{x}$ is computed as a solution to (12) where $m=2$, $z=0$, $r_{1}=0$, $r_{2}=0$,

The restored image is displayed in Fig. 1(d). Fig. 2 shows the convergence profile of the algorithm. We plot the evolution of the normalized Euclidean distance (in log scale) between the iterates and $\widetilde{x}$ in terms of computational time (Matlab R2011b codes running on a single-core Intel i7-2620M CPU@2.7 GHz with 8 GB of RAM). An approximation of $\widetilde{x}$ obtained after 5000 iterations is used. This result illustrates the fact that an appropriate choice of the metric may be beneficial in terms of speed of convergence.

References