Graph cut optimization

A pseudo-Boolean function ${\displaystyle f:\{0,1\}^{n}\to \mathbb {R} }$ is said to be representable if there exists a graph ${\displaystyle G=(V,E)}$ with non-negative weights and with source and sink nodes ${\displaystyle s}$ and ${\displaystyle t}$ respectively, and there exists a set of nodes ${\displaystyle V_{0}=\{v_{1},\dots ,v_{n}\}\subset V-\{s,t\}}$ such that, for each tuple of values ${\displaystyle (x_{1},\dots ,x_{n})\in \{0,1\}^{n}}$ assigned to the variables, ${\displaystyle f(x_{1},\dots ,x_{n})}$ equals (up to a constant) the value of the flow determined by a minimum cut ${\displaystyle C=(S,T)}$ of the graph ${\displaystyle G}$ such that ${\displaystyle v_{i}\in S}$ if ${\displaystyle x_{i}=0}$ and ${\displaystyle v_{i}\in T}$ if ${\displaystyle x_{i}=1}$.[8]

It is possible to classify pseudo-Boolean functions according to their order, determined by the maximum number of variables contributing to each single term. All first order functions, where each term depends upon at most one variable, are always representable. Quadratic functions

${\displaystyle f(\mathbf {x} )=w_{0}+\sum _{i}w_{i}(x_{i})+\sum _{i<j}w_{ij}(x_{i},x_{j}).}$

are representable if and only if they are submodular, i.e. for each quadratic term ${\displaystyle w_{ij}}$ the following condition is satisfied

${\displaystyle w_{ij}(0,0)+w_{ij}(1,1)\leq w_{ij}(0,1)+w_{ij}(1,0).}$

Cubic functions

${\displaystyle f(\mathbf {x} )=w_{0}+\sum _{i}w_{i}(x_{i})+\sum _{i<j}w_{ij}(x_{i},x_{j})+\sum _{i<j<k}w_{ijk}(x_{i},x_{j},x_{k})}$

are representable if and only if they are regular, i.e. all possible binary projections to two variables, obtained by fixing the value of the remaining variable, are submodular. For higher-order functions, regularity is a necessary condition for representability.[8]

Graph construction for a representable function is simplified by the fact that the sum of two representable functions ${\displaystyle f'}$ and ${\displaystyle f''}$ is representable, and its graph ${\displaystyle G=(V'\cup V'',E'\cup E'')}$ is the union of the graphs ${\displaystyle G'=(V',E')}$ and ${\displaystyle G''=(V'',E'')}$ representing the two functions. Such theorem allows to build separate graphs representing each term and combine them to obtain a graph representing the entire function.[8]

The graph representing a quadratic function of ${\displaystyle n}$ variables contains ${\displaystyle n+2}$ vertices, two of them representing the source and sink and the others representing the variables. When representing higher-order functions, the graph contains auxiliary nodes that allow to model higher-order interactions.

Binary terms

Example of a graph representing a quadratic term ${\displaystyle w_{ij}(x_{i},x_{j})}$ in case ${\displaystyle w_{ij}(1,0)-w_{ij}(0,0)>0}$ and ${\displaystyle w_{ij}(1,1)-w_{ij}(1,0)<0}$.

A quadratic (or binary) term ${\displaystyle w_{ij}}$ can be represented by a graph containing two non-terminal nodes ${\displaystyle v_{i}}$ and ${\displaystyle v_{j}}$. The term can be rewritten as

${\displaystyle w_{ij}(x_{i},x_{j})=w_{ij}(0,0)+k_{i}x_{i}+k_{j}x_{j}+k_{ij}\left((1-x_{i})x_{j}+x_{i}(1-x_{j})\right)}$

with

${\displaystyle {\begin{aligned}k_{i}&={\frac {1}{2}}(w_{ij}(1,0)-w_{ij}(0,0))\\k_{j}&={\frac {1}{2}}(w_{ij}(1,1)-w_{ij}(1,0))\\k_{ij}&={\frac {1}{2}}(w_{ij}(0,1)+w_{ij}(1,0)-w_{ij}(0,0)-w_{ij}(1,1)).\end{aligned}}}$

In this expression, the first term is constant and it is not represented by any edge, the two following terms depend on one variable and are represented by one edge, as shown in the previous section for unary terms, while the third term is represented by an edge ${\displaystyle v_{i}\rightarrow v_{j}}$ with weight ${\displaystyle w_{ij}(0,1)+w_{ij}(1,0)-w_{ij}(0,0)-w_{ij}(1,1)}$ (submodularity guarantees that the weight is non-negative).[8]

Ternary terms

A cubic (or ternary) term ${\displaystyle w_{ijk}}$ can be represented by a graph with four non-terminal nodes, three of them (${\displaystyle v_{i}}$, ${\displaystyle v_{j}}$ and ${\displaystyle v_{k}}$) associated to the three variables plus one fourth auxiliary node ${\displaystyle v_{ijk}}$.[note 1] A generic ternary term can be rewritten as the sum of a constant, three unary terms, three binary terms, and a ternary term in simplified form. There may be two different cases, according to the sign of ${\displaystyle p=w_{ijk}(0,0,0)+w_{ijk}(0,1,1)+w_{ijk}(1,0,1)+w_{ijk}(1,1,0)}$. If ${\displaystyle p>0}$ then

${\displaystyle w_{ijk}(x_{i},x_{j},x_{k})=w_{ijk}(0,0,0)+p_{1}(x_{i}-1)+p_{2}(x_{j}-1)+p_{3}(x_{k}-1)+p_{23}(x_{j}-1)x_{k}+p_{31}x_{i}(x_{k}-1)+p_{12}(x_{i}-1)x_{j}-px_{i}x_{j}x_{k}}$

Example of a graph representing the ternary term ${\displaystyle px_{i}x_{j}x_{k}}$ when ${\displaystyle p>0}$ (left) and when ${\displaystyle p<0}$ (right).

with

${\displaystyle {\begin{aligned}p_{1}&=w_{ijk}(1,0,1)-w_{ijk}(0,0,1)\\p_{2}&=w_{ijk}(1,1,0)-w_{ijk}(1,0,1)\\p_{3}&=w_{ijk}(0,1,1)-w_{ijk}(0,1,0)\\p_{23}&=w_{ijk}(0,0,1)+w_{ijk}(0,1,0)-w_{ijk}(0,0,0)-w_{ijk}(0,1,1)\\p_{31}&=w_{ijk}(0,0,1)+w_{ijk}(1,0,0)-w_{ijk}(0,0,0)-w_{ijk}(1,0,1)\\p_{12}&=w_{ijk}(0,1,0)+w_{ijk}(1,0,0)-w_{ijk}(0,0,0)-w_{ijk}(1,1,0).\end{aligned}}}$

If ${\displaystyle p<0}$ the construction is similarly, but the variables will have opposite value. If the function is regular, then all its projections of two variables will be submodular, implying that ${\displaystyle p_{23}}$, ${\displaystyle p_{31}}$ and ${\displaystyle p_{12}}$ are positive and then all terms in the new representation are submodular.

In this decomposition, the constant, unary and binary terms can be represented as shown in the previous sections. If ${\displaystyle p>0}$ the ternary term can be represented with a graph with four edges ${\displaystyle v_{i}\rightarrow v_{ijk}}$, ${\displaystyle v_{j}\rightarrow v_{ijk}}$, ${\displaystyle v_{k}\rightarrow v_{ijk}}$, ${\displaystyle v_{ijk}\rightarrow t}$, all with weight ${\displaystyle p}$, while if ${\displaystyle p<0}$ the term can be represented by four edges ${\displaystyle v_{ijk}\rightarrow v_{i}}$, ${\displaystyle v_{ijk}\rightarrow v_{j}}$, ${\displaystyle v_{ijk}\rightarrow v_{k}}$, ${\displaystyle s\rightarrow v_{ijk}}$ with weight ${\displaystyle -p}$.[8]

After building a graph representing a pseudo-Boolean function, it is possible to compute a minimum cut using one among the various algorithms developed for flow networks, such as Ford–Fulkerson, Edmonds–Karp, and Boykov–Kolmogorov algorithm. The result is a partition of the graph in two connected components ${\displaystyle S}$ and ${\displaystyle T}$ such that ${\displaystyle s\in S}$ and ${\displaystyle t\in T}$, and the function attains its global minimum when ${\displaystyle x_{i}=0}$ for each ${\displaystyle i}$ such that the corresponding node ${\displaystyle v_{i}\in S}$, and ${\displaystyle x_{i}=1}$ for each ${\displaystyle i}$ such that the corresponding node ${\displaystyle v_{i}\in T}$.

Max-flow algorithms such as Boykov–Kolmogorov's are very efficient in practice for sequential computation, but they are difficult to parallelise, making them not suitable for distributed computing applications and preventing them from exploiting the potential of modern CPUs. Parallel max-flow algorithms were developed, such as push-relabel[9] and jump-flood,[1] that can also take advantage of hardware acceleration in GPGPU implementations.[10][1][11]

The previous construction allows global optimization of pseudo-Boolean functions only, but it can be extended to quadratic functions of discrete variables with a finite number of values, in the form

${\displaystyle f(\mathbf {x} )=\sum _{i\in V}D(x_{i})+\sum _{(i,j)\in E}S(x_{i},x_{j})}$

where ${\displaystyle E\subseteq V\times V}$ and ${\displaystyle x_{i}\in \Lambda =\{1,\dots ,k\}}$. The function ${\displaystyle D(x_{i})}$ represents the unary contribution of each variable (often referred as data term), while the function ${\displaystyle S(x_{i},x_{j})}$ represents binary interactions between variables (smoothness term). In the general case, optimization of such functions is a NP-hard problem, and stochastic optimization methods such as simulated annealing are sensitive to local minima and in practice they can generate arbitrarily sub-optimal results.[note 2] With graph cuts it is possible to construct move-making algorithms that allow to reach in polynomial time a local minima with strong optimality properties for a wide family of quadratic functions of practical interest (when the binary interaction ${\displaystyle S(x_{i},x_{j})}$ is a metric or a semimetric), such that the value of the function in the solution lies within a constant and known factor from the global optimum.[12]

Given a function ${\displaystyle f:\Lambda ^{n}\to \mathbb {R} }$ with ${\displaystyle \Lambda =\{1,\dots ,k\}}$, and a certain assignment of values ${\displaystyle \mathbf {x} =(x_{1},\dots ,x_{n})\in \Lambda ^{n}}$ to the variables, it is possible to associate each assignment ${\displaystyle \mathbf {x} }$ to a partition ${\displaystyle P=\{P_{l}|l\in \Lambda \}}$ of the set of variables, such that, ${\displaystyle P_{l}=\{x_{i}|x_{i}=l\in \Lambda \}}$. Give two distinct assignments ${\displaystyle P}$ and ${\displaystyle P'}$ and a value ${\displaystyle \alpha \in \Lambda }$, a move that transforms ${\displaystyle P}$ into ${\displaystyle P'}$ is said to be an ${\displaystyle \alpha }$-expansion if ${\displaystyle P_{\alpha }\subset P'_{\alpha }}$ and ${\displaystyle P'_{l}\subset P_{l}\;\forall l\in \Lambda -\{\alpha \}}$. Given a couple of values ${\displaystyle \alpha }$ and ${\displaystyle \beta }$, a move is said to be an ${\displaystyle \alpha \beta }$-swap if ${\displaystyle P_{l}=P'_{l}\;\forall l\in \Lambda -\{\alpha ,\beta \}}$. Intuitively, an ${\displaystyle \alpha }$-expansion move from ${\displaystyle \mathbf {x} }$ assigns the value of ${\displaystyle \alpha }$ to some variables that have a different value in ${\displaystyle \mathbf {x} }$, while an ${\displaystyle \alpha \beta }$-swap move assigns ${\displaystyle \alpha }$ to some variables that have value ${\displaystyle \beta }$ in ${\displaystyle \mathbf {x} }$ and vice versa.

For each iteration, the ${\displaystyle \alpha }$-expansion algorithm computes, for each possible value ${\displaystyle \alpha }$, the minimum of the function among all assignments ${\displaystyle \mathrm {A} (\mathbf {x} )}$ that can be reached with a single ${\displaystyle \alpha }$-expansion move from the current temporary solution ${\displaystyle \mathbf {x} }$, and takes it as the new temporary solution.

${\displaystyle \mathbf {x} :={\text{arbitrary value in }}\Lambda ^{n}}$
${\displaystyle {\text{exit}}:=0}$
while ${\displaystyle {\text{exit}}\neq 1}$:
    ${\displaystyle {\text{exit}}=1}$
    foreach ${\displaystyle \alpha \in \Lambda }$:
        ${\displaystyle \mathbf {\hat {x}} :=\arg \min _{\mathbf {y} \in \mathrm {A} (\mathbf {x} )}f(\mathbf {y} )}$
        if ${\displaystyle f(\mathbf {\hat {x}} )<f(\mathbf {x} )}$:
            ${\displaystyle \mathbf {x} =\mathbf {\hat {x}} }$
            ${\displaystyle {\text{exit}}:=0}$

The ${\displaystyle \alpha \beta }$-swap algorithm is similar, but it searches for the minimum among all assignments ${\displaystyle \mathrm {A} \mathrm {B} (\mathbf {x} )}$ reachable with a single ${\displaystyle \alpha \beta }$-swap move from ${\displaystyle \mathbf {x} }$.

${\displaystyle \mathbf {x} :={\text{arbitrary value in }}\Lambda ^{n}}$
${\displaystyle {\text{exit}}:=0}$
while ${\displaystyle {\text{exit}}\neq 1}$:
    ${\displaystyle {\text{exit}}=1}$
    foreach ${\displaystyle (\alpha ,\beta )\in \Lambda ^{2}}$:
        ${\displaystyle \mathbf {\hat {x}} :=\arg \min _{\mathbf {y} \in \mathrm {A} \mathrm {B} (\mathbf {x} )}f(\mathbf {y} )}$
        if ${\displaystyle f(\mathbf {\hat {x}} )<f(\mathbf {x} )}$:
            ${\displaystyle \mathbf {x} =\mathbf {\hat {x}} }$
            ${\displaystyle {\text{exit}}:=0}$

In both cases, the optimization problem in the innermost loop can be solved exactly and efficiently with a graph cut. Both algorithms terminate certainly in a finite number of iterations of the outer loop, and in practice such number is small, with most of the improvement happening at the first iteration. The algorithms can generate different solutions depending on the initial guess, but in practice they are robust with respect to initialisation, and starting with a point where all variables are assigned to the same random value is usually sufficient to produce good quality results.[12]

The solution generated by such algorithms is not necessarily a global optimum, but it has strong guarantees of optimality. If ${\displaystyle S(x_{i},x_{j})}$ is a metric and ${\displaystyle \mathbf {x} }$ is a solution generated by the ${\displaystyle \alpha }$-expansion algorithm, or if ${\displaystyle S(x_{i},x_{j})}$ is a semimetric and ${\displaystyle \mathbf {x} }$ is a solution generated by the ${\displaystyle \alpha \beta }$-swap algorithm, then ${\displaystyle f(\mathbf {x} )}$ lies within a known and constant factor from the global minimum ${\displaystyle f(\mathbf {x} ^{*})}$:[12]

${\displaystyle f(\mathbf {x} )\leq 2{\frac {\max _{\alpha \neq \beta \in \Lambda }S(\alpha ,\beta )}{\min _{\alpha \neq \beta \in \Lambda }S(\alpha ,\beta )}}f(\mathbf {x} ^{*}).}$

Generally speaking, the problem of optimizing a non-submodular pseudo-Boolean function is NP-hard and cannot be solved in polynomial time with a simple graph cut. The simplest approach is to approximate the function with a similar but submodular one, for instance truncating all non-submodular terms or replacing them with similar submodular expressions. Such approach is generally sub-optimal, and it produces acceptable results only if the number of non-submodular terms is relatively small.[13]

In case of quadratic non-submodular functions, it is possible to compute in polynomial time a partial solution using algorithms such as QPBO.[13] Higher-order functions can be reduced in polynomial time to a quadratic form that can be optimised with QPBO.[14]

Quadratic functions are extensively studied and were characterised in detail, but more general results were derived also for higher-order functions. While quadratic functions can indeed model many problems of practical interest, they are limited by the fact they can represent only binary interactions between variables. The possibility to capture higher-order interactions allows to better capture the nature of the problem and it can provide higher quality results that could be difficult to achieve with quadratic models. For instance in computer vision applications, where each variable represents a pixel or voxel of the image, higher-order interactions can be used to model texture information, that would be difficult to capture using only quadratic functions.[15]

Sufficient conditions analogous to submodularity were developed to characterise higher-order pseudo-Boolean functions that can be optimised in polynomial time,[16] and there exists algorithms analogous to ${\displaystyle \alpha }$-expansion and ${\displaystyle \alpha \beta }$-swap for some families of higher-order functions.[15] The problem is NP-hard in the general case, and approximate methods were developed for fast optimization of functions that do not satisfy such conditions.[16][17]

• Boykov, Yuri; Veksler, Olga; Zabih, Ramin (2001). "Fast approximate energy minimization via graph cuts". IEEE Transactions on Pattern Analysis and Machine Intelligence. 23 (11): 1222–1239. CiteSeerX 10.1.1.439.2071. doi:10.1109/34.969114.
• Freedman, Daniel; Drineas, Petros (2005). Energy minimization via graph cuts: Settling what is possible (PDF). IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Vol. 2. pp. 939–946.
• Goldberg, Andrew V; Tarjan, Robert E (1988). "A new approach to the maximum–flow problem" (PDF). Journal of the ACM. 35 (4): 921–940. doi:10.1145/48014.61051. S2CID 52152408.
• Ishikawa, Hiroshi (2014). Higher–Order Clique Reduction Without Auxiliary Variables (PDF). IEEE Conference on Computer Vision and Pattern Recognition. IEEE. pp. 1362–1369.
• Hong, Li; Chen, George (2004). Segment-based stereo matching using graph cuts (PDF). Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Vol. 1. pp. 74–81.
• Kohli, Pushmeet; Kumar, M. Pawan; Torr, Philip H.S. (2009). "P³ & Beyond: Move Making Algorithms for Solving Higher Order Functions" (PDF). IEEE Transactions on Pattern Analysis and Machine Intelligence. 31 (9): 1645–1656. doi:10.1109/tpami.2008.217. PMID 19574624. S2CID 91470.
• Kim, Junhwan; Kolmogorov, Vladimir; Zabih, Ramin (2003). Visual correspondence using energy minimization and mutual information (PDF). Proceedings of the Ninth IEEE International Conference on Computer Vision. pp. 1033–1040.
• Kohli, Pushmeet; Ladicky, Lubor; Torr, PHS (2008). Graph cuts for minimizing robust higher order potentials (PDF) (Technical report). Oxford Brookes University. pp. 1–9.
• Kolmogorov, Vladimir; Rother, Carsten (2007). "Minimizing Nonsubmodular Functions: A Review". IEEE Transactions on Pattern Analysis and Machine Intelligence. 29 (7): 1274–1279. doi:10.1109/tpami.2007.1031. PMID 17496384. S2CID 15319364.
• Kolmogorov, Vladimir; Zabin, Ramin (2004). "What energy functions can be minimized via graph cuts?" (PDF). IEEE Transactions on Pattern Analysis and Machine Intelligence. 26 (2): 1645–1656. doi:10.1109/TPAMI.2004.1262177. hdl:1813/5842. PMID 15376891.
• Lombaert, Herve; Cheriet, Farida (2012). Simultaneous image de-noising and registration using graph cuts: Application to corrupted medical images (PDF). 11th International Conference on Information Science, Signal Processing and their Applications. pp. 264–268.
• Peng, Yi; Chen, Li; Ou-Yang, Fang-Xin; Chen, Wei; Yong, Jun-Hai (2015). "JF-Cut: a parallel graph cut approach for large-scale image and video". IEEE Transactions on Image Processing. 24 (2): 655–666. Bibcode:2015ITIP...24..655P. doi:10.1109/TIP.2014.2378060. PMID 25494510. S2CID 1665580.
• Rother, Carsten; Kolmogorov, Vladimir; Blake, Andrew (2004). Grabcut: Interactive foreground extraction using iterated graph cuts (PDF). ACM transactions on graphics. Vol. 23. pp. 309–314.
• So, Ronald WK; Tang, Tommy WH; Chung, Albert CS (2011). "Non-rigid image registration of brain magnetic resonance images using graph-cuts". Pattern Recognition. 44 (10–11): 2450–2467. doi:10.1016/j.patcog.2011.04.008.
• Stich, Timo (2009). Graph Cuts with CUDA (PDF). GPU Technology Conference.
• Tang, Tommy WH; Chung, Albert CS (2007). Non-rigid image registration using graph-cuts (PDF). International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 916–924. doi:10.1007/978-3-540-75757-3_111.
• Vineet, Vibhav; Narayanan, PJ (2008). CUDA cuts: Fast graph cuts on the GPU (PDF). IEEE Computer Society Conference on Computer Vision and Pattern Recognition Workshops. pp. 1–8.

• Implementation (C++) of several graph cut algorithms by Vladimir Kolmogorov.
• GCO, graph cut optimization library by Olga Veksler and Andrew Delong.