Playing with Duality: An Overview of Recent Primal-Dual Approaches for Solving Large-Scale Optimization Problems

I Motivation and importance of the topic

Optimization is an extremely popular paradigm which constitutes the backbone of many branches of applied mathematics and engineeering, such as signal processing, computer vision, machine learning, inverse problems, and network communications, to mention just a few. The popularity of optimization approaches often stems from the fact that many problems from the above fields are typically characterized by a lack of closed form solutions and by uncertainties. In signal and image processing, for instance, uncertainties can be introduced due to noise, sensor imperfectness, or ambiguities that are often inherent in the visual interpretation. As a result, perfect or exact solutions hardly exist, whereas inexact but optimal (in a statistical or an application-specific sense) solutions and their efficient computation is what one aims at. At the same time, one important characteristic that is nowadays shared by increasingly many optimization problems encountered in the above areas is the fact that these problems are often of very large scale. A good example is the field of computer vision where one often needs to solve low level problems that require associating at least one (and typically more than one) variable to each pixel of an image (or even worse of an image sequence as in the case of video) . This leads to problems that easily can contain millions of variables, which are therefore the norm rather than the exception in this context.

Similarly, in fields like machine learning , due to the great ease with which data can now be collected and stored, quite often one has to cope with truly massive datasets and to train very large models, which thus naturally lead to optimization problems of very high dimensionality . Of course, a similar situation arises in many other scientific domains, including application areas such as inverse problems (e.g., medical image reconstruction or satellite image restoration) or telecommunications (e.g., network design, network provisioning) and industrial engineering. Due to this fact, computational efficiency constitutes a major issue that needs to be thoroughly addressed. This, therefore, makes mandatory the use of tractable optimization techniques that are able to properly exploit the problem structures, but which at the same time remain applicable to a class of problems as wide as possible.

A bunch of important advances that took place in this regard over the last years concerns a particular class of optimization approaches known as primal-dual methods. As their name implies, these approaches proceed by concurrently solving a primal problem (corresponding to the original optimization task) as well as a dual formulation of this problem. As it turns out, in doing so they are able to exploit more efficiently the problem specific properties, thus offering in many cases important computational advantages, some of which are briefly mentioned next for two very broad classes of problems.

Primal-dual methods have been primarily employed in convex optimization problems where strong duality holds. They have been successfully applied to various types of nonlinear and nonsmooth cost functions that are prevalent in the above-mentioned application fields.

Many such applied problems can essentially be expressed under the form of a minimization of a sum of terms, where each term is given by the composition of a convex function with a linear operator. One first advantage of primal-dual methods pertains to the fact that they can yield very efficient splitting optimization schemes, according to which a solution to the original problem is iteratively computed through solving a sequence of easier subproblems, each one involving only one of the terms appearing in the objective function.

The resulting primal-dual splitting schemes can also handle both differentiable and nondifferentiable terms, the former by use of gradient operators (i.e., through explicit steps) and the latter by use of proximity operators (i.e., through implicit steps) . Depending on the target functions, either explicit or implicit steps may be easier to implement. Therefore, the derived optimization schemes exploit the properties of the input problem, in a flexible manner, thus leading to very efficient first-order algorithms.

Even more importantly, primal-dual techniques are able to achieve what is known as full splitting in the optimization literature, meaning that each of the operators involved in the problem (i.e., not only the gradient or proximity operators but also the involved linear operators) is used separately . As a result, no call to the inversion of a linear operator, which is an expensive operation for large scale problems, is required during the optimization process. This is an important feature which gives these methods a significant computational advantage compared with all other splitting-based approaches.

Last but not least, primal-dual methods lead to algorithms that are easily parallelizable, which is nowadays becoming increasingly important for efficiently handling high-dimensional problems.

I-2 Discrete optimization

Besides convex optimization, another important area where primal-dual methods play a prominent role is discrete optimization. This is of particular significance given that a large variety of tasks from signal processing, computer vision, and pattern recognition are formulated as discrete labeling problems, where one seeks to optimize some measure related to the quality of the labeling . This includes, for instance, tasks such as image segmentation, optical flow estimation, image denoising, stereo matching, to mention a few examples from image analysis. The resulting discrete optimization problems not only are of very large size, but also typically exhibit highly nonconvex objective functions, which are generally intricate to optimize.

Similarly to the case of convex optimization, primal-dual methods again offer many computational advantages, leading often to very fast graph-cut or message-passing-based algorithms, which are also easily parallelizable, thus providing in many cases a very efficient way for handling discrete optimization problems that are encountered in practice . Besides being efficient, they are also successful in making little compromises regarding the quality of the estimated solutions. Techniques like the so-called primal-dual schema are known to provide a principled way for deriving powerful approximation algorithms to difficult combinatorial problems, thus allowing primal-dual methods to often exhibit theoretical (i.e., worst-case) approximation properties. Furthermore, apart from the aforementioned worst-case guaranties, primal-dual algorithms can also provide (for free) per-instance approximation guaranties. This is essentially made possible by the fact that these methods are estimating not only primal but also dual solutions.

Convex optimization and discrete optimization have different background theory originally. Convex optimization may appear as the most tractable topic in optimization, for which many efficient algorithms have been developed allowing a broad class of problems to be solved. By contrast, combinatorial optimization problems are generally NP-hard. However, many convex relaxations of certain discrete problems can provide good approximate solutions to the original ones . The problems encountered in discrete optimization therefore constitute a source of inspiration for developing novel convex optimization techniques.

Goals of this tutorial paper. Based on the above observations, our objectives will be the following:

To provide a thorough introduction that intuitively explains the basic principles and ideas behind primal-dual approaches.

To describe how these methods can be employed both in the context of continuous optimization and in the context of discrete optimization.

To explain some of the recent advances that have taken place concerning primal-dual algorithms for solving large-scale optimization problems.

To detail useful connections between primal-dual methods and some widely used optimization techniques like the alternating direction method of multipliers (ADMM) .

Finally, to provide examples of useful applications in the context of image analysis and signal processing.

The remainder of the paper is structured as follows. In Section II, we introduce the necessary methodological background on optimization. Our presentation is grounded on the powerful notion of duality known as Fenchel’s duality, from which duality properties in linear programming can be deduced. We also introduce useful tools from functional analysis and convex optimization, including the notions of subgradient and subdifferential, conjugate function, and proximity operator. The following two sections explain and describe various primal-dual methods. Section III is devoted to convex optimization problems. We discuss the merits of various algorithms and explain their connections with ADMM, that we show to be a special case of primal-dual proximal method. Section IV deals with primal-dual methods for discrete optimization. We explain how to derive algorithms of this type based on the primal-dual schema which is a well-known approximation technique in combinatorial optimization, and we also present primal-dual methods based on LP relaxations and dual decomposition. In Section V, we present applications from the domains of signal processing and image analysis, including inverse problems and computer vision tasks related to Markov Random Field energy minimization. In Section VI, we finally conclude the tutorial with a brief summary and discussion.

II Optimization background

In this section, we introduce the necessary mathematical definitions and concepts used for introducing primal-dual algorithms in later sections. Although the following framework holds for general Hilbert spaces, for simplicity we will focus on the finite dimensional case.

Any vector $u$ in $\partial f(x)$ is called a subgradient of $f$ at $x$ (see Fig. 2).

Fermat’s rule states that is a subgradient of $f$ at $x$ if and only if $x$ belongs to the set of global minimizers of $f$. If $f$ is a proper convex function which is differentiable at $x$, then its subdifferential at $x$ reduces to the singleton consisting of its gradient, i.e. $\partial f(x)=\{\nabla f(x)\}$. Note that, in the nonconvex case, extended definitions of the subdifferential may be useful such as the limiting subdifferential , but this one reduces to the Moreau subdifferential when the function is convex.

II-B Proximity operator

II-C Conjugate function

A geometrical interpretation of this result is that the epigraph of any proper lower-semicontinuous convex function always is an intersection of closed half-spaces.

II-D Duality results

A wide array of problems in signal and image processing can be expressed under the following variational form:

The latter problem may be easier to solve than the former one, especially when $K$ is much smaller than $N$.

Another useful result follows from the fact that, by using the definition of the conjugate function of $g$, Problem (10) can be reexpressed as the following saddle-point problem:

II-E Duality in linear programming

In linear programming (LP) , we are interested in convex optimization problems of the form:

By using the properties of the conjugate function and by setting $y=-v$, it is readily shown that the dual problem (11) can be reexpressed as

Since $f$ is a convex function, strong duality holds in LP. If $\widehat{x}=(\widehat{x}^{(j)})_{1\leq j\leq N}$ is a solution to Primal-LP, a solution $\widehat{y}=(\widehat{y}^{(i)})_{1\leq i\leq K}$ to Dual-LP can be obtained by the primal complementary slackness condition:

whereas, if $\widehat{y}$ is a solution to Dual-LP, a solution $\widehat{x}$ to Primal-LP can be obtained by the dual complementary slackness condition:

III Convex optimization algorithms

In this section, we present several primal-dual splitting methods for solving convex optimization problems, starting from the basic forms to the more sophisticated highly parallelized ones.

A wide range of convex optimization problems can be formulated as follows:

For examples, the functions $f$, $g\circ L$, and $h$ may model various data fidelity terms and regularization functions encountered in the solution of inverse problems. In particular, the Lipschitz differentiability property is satisfied for least squares criteria.

With respect to Problem (10), we have introduced an additional smooth term $h$. This may be useful in offering more flexibility for taking into account the structure of the problem of interest and the properties of the involved objective function. We will however see that not all algorithms are able to possibly take advantage of the fact that $h$ is a smooth term.

Based on the results in Section II-D and Property (I) in Table I, the dual optimization problem reads:

Note that, in the particular case when $h=0$, the inf-convolution f^{*}\mbox{\footnotesize\,\square\,}h^{*} (see the definition in Table I(I)) of the conjugate functions of $f$ and $h$ reduces to $f^{*}$ and we recover the basic form (11) of the dual problem.

It has to be mentioned that some specific forms of Problem (19) (e.g. when $g=0$) can be solved in a quite efficient manner by simpler proximal algorithms (see ) than those described in the following.

III-B ADMM

The celebrated ADMM (Alternating Direction Method of Multipliers) can be viewed as a primal-dual algorithm. This algorithm belongs to the class of augmented Lagrangian methods since a possible way of deriving this algorithm consists of looking for a saddle point of an augmented version of the classical Lagrange function . This augmented Lagrangian is defined as

where $\gamma\in\,\left]0,+\infty\right[$ and $\gamma z$ corresponds to a Lagrange multiplier. ADMM simply splits the step of minimizing the augmented Lagrangian with respect to $(x,y)$ by alternating between the two variables, while a gradient ascent is performed with respect to the variable $z$. The resulting iterations are given in Algorithm 1.

This algorithm has been known for a long time although it has attracted recently much interest in the signal and image processing community (see e.g. ). A condition for the convergence of ADMM is as follows:

A convergence rate analysis is conducted in .

It must be emphasized that ADMM is equivalent to the application of the Douglas-Rachford algorithm , another famous algorithm in convex optimization, to the dual problem. Other primal-dual algorithms can be deduced from the Douglas-Rachford iteration or an augmented Lagrangian approach .

III-C Methods based on a Forward-Backward approach

Note that when $L=0$ and $g^{*}=0$ the basic form of the forward-backward algorithm (also called the proximal gradient algorithm) is recovered, a popular example of which is the iterative soft-thresholding algorithm .

Convergence guarantees were established in , as well as for a more general version of this algorithm in :

Algorithm 2 also constitutes a generalization of (designated by some authors as PDHG, Primal-Dual Hybrid Gradient). Preconditioned or adaptive versions of this algorithm were proposed in which may accelerate its convergence. Convergence rate results were also recently derived in .

Another primal-dual method (see Algorithm 5) was proposed in which also results from a forward-backward approach . This algorithm is restricted to the case when $f=0$ in Problem (19).

As shown by the next convergence result, the conditions on the step-sizes $\tau$ and $\sigma$ are less restrictive than for Algorithm 2.

It must be emphasized that Algorithms 2-5 present two interesting features which are very useful in practice. At first, they allow to deal with the functions involved in the optimization problem at hand either through their proximity operator or through their gradient. Indeed, for some functions, especially non differentiable or non finite ones, the proximity operator can be a very powerful tool but, for some smooth functions (e.g. the Poisson-Gauss neg-log-likelihood ) the gradient may be easier to handle. Secondly, these algorithms do not require to invert any matrix, but only to apply $L$ and its adjoint. This advantage is of main interest when large-size problems have to be solved for which the inverse of $L$ (or $L^{\top}L$) does not exist or it has a no tractable expression.

III-D Methods based on a Forward-Backward-Forward approach

Primal-dual methods based on a forward-backward-forward approach were among the first primal-dual proximal methods proposed in the optimization literature, inspired from the seminal work in . They were first developed in the case when $h=0$ , then extended to more general scenarios in (see also for further refinements).

The convergence of the algorithm is guaranteed by the following result:

Algorithm 6 is often refered to as the M+LFBF (Monotone+Lipschitz Forward Backward Forward) algorithm. It enjoys the same advantages as FB-based primal-dual algorithms we have seen before. It however makes it possible to compute the proximity operators of scaled versions of functions $f$ and $g^{*}$ in parallel. In addition, the choice of its parameters in order to satisfy convergence conditions may appear more intuitive than for Algorithms 2-4. With respect to FB-based algorithms, an extra forward step however needs to be performed. This may lead to a slower convergence if, for example, the computational cost of the gradient is high and an iteration of a FB-based algorithm is at least as efficient as an iteration of Algorithm 6.

III-E A projection-based primal-dual algorithm

Another primal-dual algorithm was recently proposed in which relies on iterative projections onto half-spaces including the set of Kuhn-Tucker points (see Algorithm 7).

We have then the following convergence result:

Although few numerical experiments have been performed with this algorithm, one of its potential advantages is that it introduces few constraints on the choice of the parameters $\gamma_{n}$, $\mu_{n}$ and $\lambda_{n}$ at iteration $n$ and that it does not require any knowledge on the norm of the matrix $L$. Nonetheless, the use of this algorithm does not allow us to exploit the fact that $h$ is a differentiable function.

III-F Extensions

More generally, one may be interested in more challenging convex optimization problems of the form:

IV Discrete optimization algorithms

As already mentioned in the introduction, another common class of problems in signal processing and image analysis are discrete optimization problems, for which primal-dual algorithms also play an important role. Problems of this type are often stated as integer linear programs (ILPs), which can be expressed under the following form:

where $L=(L^{(i,j)})_{1\leq i\leq K,1\leq j\leq N}$ represents a matrix of size $K\times N$, and $b=(b^{(i)})_{1\leq i\leq K}$, $c=(c^{(j)})_{1\leq j\leq N}$ are column vectors of size $K$ and $N$, respectively. Note that integer linear programming provides a very general formulation suitable for modeling a very broad range of problems, and will thus form the setting that we will consider hereafter. Among the problems encountered in practice, many of them lead to a Primal-ILP that is NP-hard to solve. In such cases, a principled approach for finding an approximate solution is through the use of convex relaxations (see framebox), where the original NP-hard problem is approximated with a surrogate one (the so-called relaxed problem), which is convex and thus much easier to solve. The premise is the following: to the extent that the surrogate problem provides a reasonably good approximation to the original optimization task, one can expect to obtain an approximately optimal solution for the latter by essentially making use of or solving the former.

The type of relaxations that are typically preferred in large scale discrete optimization are based on linear programming, involving the minimization of a linear function subject to linear inequality constraints. These can be naturally obtained by simply relaxing the integrality constraints of Primal-ILP, thus leading to the relaxed primal problem (14) as well as its dual (16). It should be noted that the use of LP-relaxations is often dictated by the need of maintaining a reasonable computational cost. Although more powerful convex relaxations do exist in many cases, these may become intractable as the number of variables grows larger, especially for Semidefinite Programming (SDP) or Second-Order Cone Programming (SOCP) relaxations.

Based on the above observations, in the following we aim to present some very general primal-dual optimization strategies that can be used in this context, focusing a lot on their underlying principles, which are based on two powerful techniques, the so-called primal-dual schema and dual decomposition. As we shall see, in order to estimate an approximate solution to Primal-ILP, both approaches make heavy use of the dual of the underlying LP relaxation, i.e., Problem (16). But their strategies for doing so is quite different: the second one essentially aims at solving this dual LP (and then converting the fractional solution into an integral one, trying not to increase the cost too much in the process), whereas the first one simply uses it in the design of the algorithm.

IV-B The primal-dual schema for integer linear programming

The primal-dual schema is a well-known technique in the combinatorial optimization community that has its origins in LP duality theory. It is worth noting that it started as an exact method for solving linear programs. As such, it had initially been used in deriving exact polynomial-time algorithms for many cornerstone problems in combinatorial optimization that have a tight LP relaxation. Its first use probably goes back to Edmond’s famous Blossom algorithm for constructing maximum matchings on graphs, but it had been also applied to many other combinatorial problems including max-flow (e.g., Ford and Fulkerson’s augmenting path-based techniques for max-flow can essentially be understood in terms of this schema), shortest path, minimum branching, and minimum spanning tree . In all of these cases, the primal-dual schema is driven by the fact that optimal LP solutions should satisfy the complementary slackness conditions (see (17) and (18)). Starting with an initial primal-dual pair of feasible solutions, it therefore iteratively steers them towards satisfying these complementary slackness conditions (by trying at each step to minimize their total violation). Once this is achieved, both solutions (the primal and the dual) are guaranteed to be optimal. Moreover, since the primal is always chosen to be updated integrally during the iterations, it is ensured that an integral optimal solution is obtained at the end. A notable feature of the primal-dual method is that it often reduces the original LP, which is a weighted optimization problem, to a series of purely combinatorial unweighted ones (related to minimizing the violation of complementary slackness conditions at each step).

Interestingly, today the primal-dual schema is no longer used for providing exact algorithms. Instead, its main use concerns deriving approximation algorithms to NP-hard discrete problems that admit an ILP formulation, for which it has proved to be a very powerful and widely applicable tool. As such, it has been applied to many NP-hard combinatorial problems up to now, including set-cover, Steiner-network, scheduling, Steiner tree, feedback vertex set, facility location, to mention only a few . With regard to problems from the domains of computer vision and image analysis, the primal-dual schema has been introduced recently in , and has been used for modeling a broad class of tasks from these fields.

Then, $x$ can be shown to be a $\nu$-approximation to an unknown optimal integral solution $\widehat{x}$, i.e.

PRIMAL-DUAL PRINCIPLE IN THE DISCRETE CASE Essentially, the proof of this principle relies on the fact that the sequence of optimal costs of problems Dual-LP, Primal-LP, and Primal-ILP is increasing. Specifically, by weak LP duality, the optimal cost of Dual-LP is known to not exceed the optimal cost of Primal-LP. As a result of this fact, the cost $c^{\top}\widehat{x}$ (of an unknown optimal integral solution $\widehat{x}$) is guaranteed to be at least as large as the cost $b^{\top}y$ of any dual feasible solution $y$. On the other hand, by definition, $c^{\top}\widehat{x}$ cannot exceed the cost $c^{\top}x$ of an integral-primal feasible solution $x$. Therefore, if the gap $\Delta(y,x)$ between the costs of $y$ and $x$ is small (e.g., it holds $c^{\top}x\leq\nu\,b^{\top}y$), the same will be true for the gap $\Delta(\widehat{x},x)$ between the costs of $\widehat{x}$ and $x$ (i.e., $c^{\top}x\leq\nu\,c^{\top}\widehat{x}$), thus proving that $x$ is a $\nu$-approximation to optimal solution $\widehat{x}$.

Although the above principle lies at the heart of many primal-dual techniques (i.e., in one way or another, primal-dual methods often try to fulfill the assumptions imposed by this principle), it does not directly specify how to estimate a primal-dual pair of solutions $(x,y)$ that satisfies these assumptions. This is where the so-called relaxed complementary slackness conditions come into play, as they typically provide an alternative and more convenient (from an algorithmic viewpoint) way for generating such a pair of solutions. These conditions generalize the complementary slackness conditions associated with an arbitrary pair of primal-dual linear programs (see Section II-E). The latter conditions apply only in cases when there is no duality gap, like between Primal-LP and Dual-LP, but they are not applicable to cases like Primal-ILP and Dual-LP, when a duality gap exists as a result of the integrality constraint imposed on variable $x$. As in the exact case, two types of relaxed complementary slackness conditions exist, depending on whether the primal or dual variables are checked for being zero.

where J_{x}=\big{\{}{j\in\{1,\ldots,N\}}~{}\big{|}~{}{x^{(j)}>0}\big{\}}.

where I_{y}=\big{\{}{i\in\{1,\ldots,K\}}~{}\big{|}~{}{y^{(i)}>0}\big{\}}.

This result simply follows from the inequalities

Based on the above result, iterative schemes can be devised yielding a primal-dual $\nu$-approximation algorithm. For example, we can employ the following algorithm:

Note that, in this scheme, primal solutions are always updated integrally. Also, note that, when applying the primal-dual schema, different implementation strategies are possible. The strategy described in Algorithm 8 is to maintain feasible primal-dual solutions $(x_{n},y_{n})$ at iteration $n$, and iteratively improve how tightly the (primal or dual) complementary slackness conditions get satisfied. This is performed through the introduction of slackness variables $(q^{(i)})_{i\in I_{y_{n}}}$ and $(r^{(j)})_{j\in J_{x_{n}}}$ the sums of which measure the degrees of violation of each relaxed slackness condition and have thus to be minimized. Alternatively, for example, we can opt to maintain solutions $(x_{n},y_{n})$ that satisfy the relaxed complementary slackness conditions, but may be infeasible, and iteratively improve the feasibility of the generated solutions. For instance, if we start with a feasible dual solution but with an infeasible primal solution, such a scheme would result into improving the feasibility of the primal solution, as well as the optimality of the dual solution at each iteration, ensuring that a feasible primal solution is obtained at the end. No matter which one of the above two strategies we choose to follow, the end result will be to gradually bring the primal and dual costs $c^{\top}x_{n}$ and $b^{\top}y_{n}$ closer and closer together so that asymptotically the primal-dual principle gets satisfied with the desired approximation factor. Essentially, at each iteration, through the coupling by the complementary slackness conditions the current primal solution is used to improve the dual, and vice versa.

The above problem can be expressed as the following ILP:

where indicator variables $(x^{(j)})_{1\leq j\leq N}$ are used for determining if a set in $\mathcal{S}$ has been included in the set cover or not, and (37) ensures that each one of the elements of $\mathcal{V}$ is contained in at least one of the sets that were chosen for participating to the set cover.

An LP-relaxation for this problem is obtained by simply replacing the Boolean constraint with the constraint $x\in[0,+\infty[^{N}$. The dual of this LP relaxation is given by the following linear program:

Primal Complementary Slackness Conditions

Relaxed Dual Complementary Slackness Conditions (with relaxation factor $F_{\max}$)

A set $S_{j}$ with $j\in\{1,\ldots,N\}$ for which $\sum_{\begin{subarray}{c}i\in\{1,\ldots,K\}\\ \upsilon^{(i)}\in S_{j}\end{subarray}}y^{(i)}=\varphi(S_{j})$ will be called packed. Based on this definition, and given that the primal variables $(x^{(j)})_{1\leq j\leq N}$ are always kept integral (i.e., either or $1$) during the primal-dual schema, Conditions (40) basically say that only packed sets can be included in the set cover (note that overpacked sets are already forbidden by feasibility constraints (39)). Similarly, Conditions (41) require that an element $\upsilon^{(i)}$ with $i\in\{1,\ldots,K\}$ associated with a nonzero dual variable $y^{(i)}$ should not be covered more than $F_{\max}$ times, which is, of course, trivially satisfied given that $F_{\max}$ represents the maximum frequency of any element in $\mathcal{V}$.

IV-C Dual decomposition

We will next examine a different approach for discrete optimization, which is based on the principle of dual decomposition . The core idea behind this principle essentially follows a divide and conquer strategy: that is, given a difficult or high-dimensional optimization problem, we decompose it into smaller easy-to-handle subproblems and then extract an overall solution by cleverly combining the solutions from these subproblems.

To explain this technique, we will consider the general problem of minimizing the energy of a discrete Markov Random Field (MRF), which is a ubiquitous problem in the fields of computer vision and image analysis (applied with great success on a wide variety of tasks from these domains such as stereo-matching, image segmentation, optical flow estimation, image restoration and inpainting, or object detection) . This problem involves a graph $G$ with vertex set $\mathcal{V}$ and edge set $\mathcal{E}$ (i.e., $G=(\mathcal{V},\mathcal{E})$) plus a finite label set $\mathcal{L}$. The goal is to find a labeling $z=(z^{(p)})_{p\in\mathcal{V}}\in\mathcal{L}^{|\mathcal{V}|}$ for the graph vertices that has minimum cost, that is

where, for every $p\in\mathcal{V}$ and $e\in\mathcal{E}$, $\varphi_{p}\colon\mathcal{L}\to\,\left]-\infty,+\infty\right[$ and $\sansmath\varphi_{e}\colon\mathcal{L}^{2}\to\,\left]-\infty,+\infty\right[$ represent the unary and pairwise costs (also known connectively as MRF potentials $\varphi=\left\{\{\varphi_{p}\}_{p\in\mathcal{V}},\{\sansmath\varphi_{e}\}_{e\in\mathcal{E}}\right\}$), and $\mathsf{z}^{(e)}$ denotes the pair of components of $z$ defined by the variables corresponding to vertices connected by $e$ (i.e., $\mathsf{z}^{(e)}=(z^{(p)},z^{(q)})$ for $e=(p,q)\in\mathcal{E}$).

The above problem is NP-hard, and much of the recent work on MRF optimization revolves around the following equivalent ILP formulation of (42) , which is the one that we will also use here:

where the set $C_{G}$ is defined for any graph $G=(\mathcal{V},\mathcal{E})$ as

In the above formulation, for every $p\in\mathcal{V}$ and $e\in\mathcal{E}$, the unary binary function $x_{p}(\cdot)$ and the pairwise binary function $\mathsf{x}_{e}(\cdot)$ indicate the labels assigned to vertex $p$ and to the pair of vertices connected by edge $e=(p^{\prime},q^{\prime})$ respectively, i.e.,

Minimizing with respect to the vector $x$ regrouping all these binary functions is equivalent to searching for an optimal binary vector of dimension $N=|\mathcal{V}||\mathcal{L}|+|\mathcal{E}||\mathcal{L}|^{2}$. The first constraints in (44) simply encode the fact that each vertex must be assigned exactly one label, whereas the rest of the constraints enforces consistency between unary functions $x_{p}(\cdot)$, $x_{q}(\cdot)$ and the pairwise function $\mathsf{x}_{e}(\cdot)$ for edge $e=(p,q)$, ensuring essentially that if $x_{p}(z^{(p)})=x_{q}(z^{(q)})=1$, then $\mathsf{x}_{e}(z^{(p)},z^{(q)})=1$.

As mentioned above, our goal will be to decompose the MRF problem (43) into easier subproblems (called slaves), which, in this case, involve optimizing MRFs defined on subgraphs of $G$. More specifically, let $\{G_{m}=(\mathcal{V}_{m},\mathcal{E}_{m})\}_{1\leq m\leq M}$ be a set of subgraphs that form a decomposition of $G=(\mathcal{V},\mathcal{E})$ (i.e., $\displaystyle\cup_{m=1}^{M}\mathcal{V}_{m}=\mathcal{V}$, $\displaystyle\cup_{m=1}^{M}\mathcal{E}_{m}=\mathcal{E}$). On each of these subgraphs, we define a local MRF with corresponding (unary and pairwise) potentials $\varphi^{m}=\left\{\{\varphi^{m}_{p}\}_{p\in\mathcal{V}_{m}},\{\sansmath\varphi^{m}_{e}\}_{e\in\mathcal{E}_{m}}\right\}$, whose cost function $f^{m}(x;\varphi^{m})$ is thus given by

Moreover, the sum (over $m$) of the potential functions $\varphi^{m}$ is ensured to give back the potentials $\varphi$ of the original MRF on $G$, i.e.,For instance, to ensure (48) we can simply set: $(\forall m\in\{1,\ldots,M\})$ $\varphi^{m}_{p}=\frac{\varphi_{p}}{|\{|m^{\prime}\mid p\in\mathcal{V}_{m^{\prime}}\}|}$ and $\sansmath\varphi^{m}_{e}=\frac{\sansmath\varphi_{e}}{|\{m^{\prime}\mid e\in\mathcal{E}_{m^{\prime}}\}|}$.

This guarantees that $f=\sum_{m=1}^{M}f^{m}$, thus allowing us to re-express problem (43) as follows

An assumption that often holds in practice is that minimizing separately each of the $f^{m}$ (over $x$) is easy, but minimizing their sum is hard. Therefore, to leverage this fact, we introduce, for every $m\in\{1,\ldots,M\}$, an auxiliary copy $x^{m}\in C_{G_{m}}$ for the variables of the local MRF defined on $G_{m}$, which are thus constrained to coincide with the corresponding variables in vector $x$, i.e., it holds $x^{m}=x_{|G_{m}}$, where $x_{|G_{m}}$ is used to denote the subvector of $x$ containing only those variables associated with vertices and edges of subgraph $G_{m}$. In this way, Problem (49) can be transformed into

By considering the dual of (IV-C), using a technique similar to the one described in framebox on page II-D, and noticing that

we finally end up with the following problem:

where, for every $m\in\{1,\ldots,M\}$, the dual variable $v^{m}$ consists of $\left\{\{v^{m}_{p}\}_{p\in\mathcal{V}_{m}},\{\mathsf{v}^{m}_{e}\}_{e\in\mathcal{E}_{m}}\right\}$ similarly to $\varphi^{m}$, and function $h^{m}$ is related to the following optimization of a slave MRF on $G_{m}$:

MASTER-SLAVE COMMUNICATION During dual decomposition a communication between a master process and the slaves (local subproblems) can be thought of as taking place, which can also be interpreted as a resource allocation/pricing stage. Resource allocation: At each iteration, the master assigns new MRF potentials (i.e., resources) $(\varphi^{m})_{1\leq m\leq M}$ to the slaves based on the current local solutions $(\widehat{x}^{m})_{1\leq m\leq M}$. Pricing: The slaves respond by adjusting their local solutions $(\widehat{x}^{m})_{1\leq m\leq M}$ (i.e., the prices) so as to maximize their welfares based on the newly assigned resources $(\widehat{x}^{m})_{1\leq m\leq M}$.

DECOMPOSITIONS AND RELAXATIONS Different decompositions can lead to different relaxations and/or can affect the speed of convergence. We show below, for instance, 3 possible decompositions for an MRF assumed to be defined on a $5\times 5$ image grid. Decompositions $\{G^{1}_{m}\},\{G^{2}_{m}\},\{G^{3}_{m}\}$ consist respectively of one suproblem per row and column, one subproblem per edge, and one subproblem per $2\times 2$ subgrid of the original $5\times 5$ grid. Both $\{G^{1}_{m}\}$ and $\{G^{2}_{m}\}$ (due to using solely subgraphs that are trees) lead to the same LP relaxation of (43), whereas $\{G^{3}_{m}\}$ leads to a relaxation that is tighter (due to containing loopy subgraphs). On the other hand, decomposition $\{G^{1}_{m}\}$ leads to faster convergence compared with $\{G^{2}_{m}\}$ due to using larger subgraphs that allow a faster propagation of information during message-passing.

Interestingly, if we choose to use a decomposition consisting only of subgraphs that are trees, then the resulting relaxation can be shown to actually coincide with the standard LP-relaxation of linear integer program (43) (generated by replacing the integrality constraints with non-negativity constraints on the variables). This also means that when this LP-relaxation is tight, an optimal MRF solution is computed. This, for instance, leads to the result that dual decomposition approaches can estimate a globally optimal solution for binary submodular MRFs (although it should be noted that much faster graph-cut based techniques exist for submodular problems of this type - see framebox on page IV-C). Furthermore, when using subgraphs that are trees, a minimizer to each slave problem can be computed efficiently by applying the Belief Propagation algorithm , which is a message-passing method. Therefore, in this case, Algorithm 10 essentially reduces to a continuous exchange of messages between the nodes of graph $G$. Such an algorithm relates to or generalizes various other message-passing approaches . In general, besides tree-structured subgraphs, other types of decompositions or subproblems can be used as well (such as binary planar problems, or problems on loopy subgraphs with small tree-width, for which MRF optimization can still be solved efficiently), which can lead to even tighter relaxations (see framebox on page IV-C) .

V Applications

Although the presented primal-dual algorithms can be applied virtually to any area where optimization problems have to be solved, we now mention a few common applications of these techniques.

For a long time, convex optimization approaches have been successfully used for solving inverse problems such as signal restoration, signal reconstruction, or interpolation of missing data. Most of the time, these problems are ill-posed and, in order to recover the signal of interest in a satisfactory manner, some prior information needs to be introduced. To do this, an objective function can be minimized which includes a data fidelity term modelling knowledge about the noise statistics and possibly involves a linear observation matrix (e.g. a convolutive blur), and a regularization (or penalization) term which corresponds to the additional prior information. This formulation can also often be justified statistically as the determination of a Maximum A Posteriori (MAP) estimate. In early developed methods, in particular in Tikhonov regularization, a quadratic penalty function is employed. Alternatively, hard constraints can be imposed on the solution (for example, bounds on the signal values), leading to signal feasibility problems. Nowadays, a hybrid regularization may be prefered so as to combine various kinds of regularity measures, possibly computed for different representations of the signal (Fourier, wavelets,…), some of them like total variation and its nonlocal extensions being taylored for preserving discontinuities such as image edges. In this context, constraint sets can be translated into penalization terms being equal to the indicator functions of these sets (see (2)). Altogether, these lead to global cost functions which can be quite involved, often with many variables, for which the splitting techniques described in Section III-F are very useful. An extensive literature exists on the use of ADMM methods for solving inverse problems (e.g., see ). With the advent of more recent primal-dual algorithms, many works have been mainly focused on image recovery applications . Two illustrations are now provided.

In , a generalization of the total variation is defined for an arbitrary graph in order to address a variety of inverse problems. For denoising applications, the optimization problem to be solved is of the form (19) where

Note that convex primal-dual proximal optimization algorithms have been applied to other fields than image recovery, in particular to machine learning , system identification , audio processing , optimal transport , empirical mode decomposition , seimics , database management , and data streaming over networks .

V-B Computer vision and image analysis

The great majority of problems in computer vision involve image observation data that are of very high dimensionality, inherently ambiguous, noisy, incomplete, and often only provide a partial view of the desired space. Hence, any successful model that aims to explain such data usually requires a reasonable regularization, a robust data measure, and a compact structure between the variables of interest to efficiently characterize their relationships. Probabilistic graphical models, and in particular discrete Markov Random Fields, have led to a suitable methodology for solving such visual perception problems . This type of models offer great representational power, and are able to take into account dependencies in the data, encode prior knowledge, and model (soft or hard) contextual constraints in a very efficient and modular manner. Furthermore, they offer the important ability to make use of very powerful data likelihood terms consisting of arbitrary nonconvex and non-continuous functions that are often crucial for accurately representing the problem at hand. As a result, MAP-inference for these models leads to discrete optimization problems that are (in most cases) highly nonconvex (NP-hard) and also of very large scale . These discrete problems take the form (42), where typically the unary terms $\varphi_{p}(\cdot)$ encode the data likelihood and the higher-order terms $\sansmath\varphi_{e}(\cdot)$ encode problem specific priors.

Primal-dual approaches can offer important computational advantages when dealing with such problems. One such characteristic example is the FastPD algorithm , which currently provides one of the most efficient methods for solving generic MRF optimization problems of this type, also guaranteeing at the same time the convergence to solutions that are approximately optimal. The theoretical derivation of this method relies on the use of the primal-dual schema described in Section IV, which results, in this case, into a very fast graph-cut based inference scheme that generalizes previous state-of-the-art approaches such as the $\alpha$-expansion algorithm (see Fig. 10). More generally, due to the very wide applicability of MRF models to computer vision or image analysis problems, primal-dual approaches can and have been applied to a broad class of both low-level and high-level problems from these domains, including image segmentation , stereo matching and 3D multi-view reconstruction , graph-matching , 3D surface tracking , optical flow estimation , scene understanding , image deblurring , panoramic image stitching , category-level segmentation , and motion tracking . In the following we mention very briefly just a few examples.

A primal-dual based optimization framework has been recently proposed in for the problem of deformable registration/fusion, which forms one of the most central and challenging tasks in medical image analysis. This problem consists of recovering a nonlinear dense deformation field that aligns two signals that have in general an unknown relationship both in the spatial and intensity domain. In this framework, towards dimensionality reduction on the variables, the dense registration field is first expressed using a set of control points (registration grid) and an interpolation strategy. Then, the registration cost is expressed using a discrete sum over image costs projected on the control points, and a smoothness term that penalizes local deviations on the deformation field according to a neighborhood system on the grid. One advantage of the resulting optimization framework is that it is able to encode even very complex similarity measures (such as normalized mutual information and Kullback-Leibler divergence) and therefore can be used even when seeking transformations between different modalities (inter-deformable registration). Furthermore, it admits a broad range of regularization terms, and can also be applied to both 2D-2D and 3D-3D registration, as an arbitrary underlying graph structure can be readily employed (see Fig. 11 for a result on 3D inter-subject brain registration).

where $f(u(s),s)$ is a data term favoring different depth values by measuring the absolute intensity differences of respective patches projected in the two input images, and the second term is a TV regularizer that promotes spatially smooth depth fields. The above problem is nonconvex (due to the use of the data term $f$), but it turns out that there exists an equivalent convex formulation obtained by lifting the original problem to a higher-dimensional space, in which $u$ is represented in terms of its level sets

In the above formulation, $\Sigma=\Omega\times\Gamma$, $\phi\colon\Sigma\rightarrow\{0,1\}$ is a binary function such that $\phi(s,\upsilon)$ equals $1$ if $u(s)>\upsilon$ and otherwise, and the feasible set is defined as $D=\left\{\phi\colon\Sigma\rightarrow\{0,1\}\mid(\forall s\in\Omega)\,\phi(s,\upsilon_{\rm min})=1,\phi(s,\upsilon_{\rm max})=0\right\}$. A convex relaxation of the latter problem is obtained by using D^{\prime}=\big{\{}\phi\colon\Sigma\rightarrow\mid(\forall s\in\Omega)\phi(s,\upsilon_{\rm min})=1,\phi(s,\upsilon_{\rm max})=0\big{\}} instead of $D$. A discretized form of the resulting optimization problem can be solved with the algorithms described in Section III-C. Fig. 12 shows a sample result of this approach.

Recently, primal-dual approaches have also been developed for discrete optimization problems that involve higher-order terms . They have been applied successfully to various tasks, like, for instance, in stereo matching . In this case, apart from a data term that measures similarity between corresponding pixels in two images, a discontinuity-preserving smoothness prior of the form $\sansmath\varphi(s_{1},s_{2},s_{3})=\min(|s_{1}-2s_{2}+s_{3}|,\kappa)$ with $\kappa\in\,\left]0,+\infty\right[$ has been employed as a regularizer that penalizes depth surfaces of high curvature. Indicative stereo matching results from an algorithm based on the dual decomposition principle described in Section IV-C are shown in Fig. 13.

It should be also mentioned that an advantage of all primal-dual algorithms (which is especially important for NP-hard problems) is that they also provide (for free) per-instance approximation bounds, specifying how far the cost of an estimated solution can be from the unknown optimal cost. This directly follows from the fact that these methods are computing both primal and dual solutions, which (in the case of a minimization task) provide respectively upper and lower limits to the true optimum. These approximation bounds are continuously updated throughout an algorithm execution, and thus can be directly used for assessing the performance of a primal-dual method with respect to any particular problem instance (and without essentially any extra computational cost). Moreover, often in practice, these sequences converge to a common value, which means that the corresponding estimated solutions are almost optimal (see, e.g., the plots in Fig. 13).

VI Conclusion

In this paper, we have reviewed a number of primal-dual optimization methods which can be employed for solving signal and image processing problems. The links existing between convex approaches and discrete ones were little explored in the literature and one of the contributions of this paper is to put them in a unifying perspective. Although the presented algorithms have been proved to be quite effective in numerous problems, there remains much room for extending their scope to other application fields, and also for improving them so as to accelerate their convergence. In particular, the parameter choices in these methods may have a strong influence on the convergence speed and it would be thus interesting to design automatic procedures for setting these parameters. Various techniques can also be devised for designing faster variants of these methods (preconditioning, activation of blocks of variables, combination with stochastic strategies, distributed implementations…). Another issue to pay attention to is the robustness to numerical errors although it can be mentioned that most of the existing proximal algorithms are tolerant to summable errors. Concerning discrete optimization methods, we have shown that the key to success lies in tight relaxations of combinatorial NP hard problems. Extending these methods to more challenging problems, e.g. those involving higher-order Markov fields or extremely large label sets, appears to be of main interest in this area. More generally, developing primal-dual strategies that further bridge the gap between continuous and discrete approaches, as well as for solving other kinds of nonconvex optimization problems such as those encountered in phase reconstruction or blind deconvolution opens the way to appealing investigations. So, the ground is yours now to play with duality!

References