Graph cut optimization

Graph cut optimization is a combinatorial optimization method applicable to a family of functions of discrete variables, named after the concept of cut in the theory of flow networks. Thanks to the max-flow min-cut theorem, determining the minimum cut over a graph representing a flow network is equivalent to computing the maximum flow over the network. Given a pseudo-Boolean function formula_1, if it is possible to construct a flow network with positive weights such that
then it is possible to find the global optimum of formula_1 in polynomial time by computing a minimum cut of the graph. The mapping between cuts and variable assignments is done by representing each variable with one node in the graph and, given a cut, each variable will have a value of 0 if the corresponding node belongs to the component connected to the source, or 1 if it belong to the component connected to the sink.
Not all pseudo-Boolean functions can be represented by a flow network, and in the general case the global optimization problem is NP-hard. There exist sufficient conditions to characterise families of functions that can be optimised through graph cuts, such as submodular quadratic functions. Graph cut optimization can be extended to functions of discrete variables with a finite number of values, that can be approached with iterative algorithms with strong optimality properties, computing one graph cut at each iteration.
Graph cut optimization is an important tool for inference over graphical models such as Markov random fields or conditional random fields, and it has applications in computer vision problems such as image segmentation, denoising, registration and stereo matching.
Representability.
A pseudo-Boolean function formula_8 is said to be "representable" if there exists a graph formula_9 with non-negative weights and with source and sink nodes formula_10 and formula_11 respectively, and there exists a set of nodes formula_12 such that, for each tuple of values formula_13 assigned to the variables, formula_14 equals (up to a constant) the value of the flow determined by a minimum cut formula_15 of the graph formula_16 such that formula_17 if formula_18 and formula_19 if formula_20.
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
are representable if and only if they are submodular, i.e. for each quadratic term formula_22 the following condition is satisfied
Cubic functions
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.
Graph construction.
Graph construction for a representable function is simplified by the fact that the sum of two representable functions formula_25 and formula_26 is representable, and its graph formula_27 is the union of the graphs formula_28 and formula_29 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.
The graph representing a quadratic function of formula_30 variables contains formula_31 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.
Unary terms.
A unary term formula_32 depends only on one variable formula_33 and can be represented by a graph with one non-terminal node formula_34 and one edge formula_35 with weight formula_36 if formula_37, or formula_38 with weight formula_39 if formula_40.
Binary terms.
A quadratic (or binary) term formula_22 can be represented by a graph containing two non-terminal nodes formula_34 and formula_43. The term can be rewritten as
with
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 formula_46 with weight formula_47 (submodularity guarantees that the weight is non-negative).
Ternary terms.
A cubic (or ternary) term formula_48 can be represented by a graph with four non-terminal nodes, three of them (formula_34, formula_43 and formula_51) associated to the three variables plus one fourth auxiliary node formula_52. 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 formula_53. If formula_54 then
with
If formula_57 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 formula_58, formula_59 and formula_60 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 formula_54 the ternary term can be represented with a graph with four edges formula_62, formula_63, formula_64, formula_65, all with weight formula_66, while if formula_57 the term can be represented by four edges formula_68, formula_69, formula_70, formula_71 with weight formula_72.
Minimum cut.
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 formula_73 and formula_74 such that formula_75 and formula_76, and the function attains its global minimum when formula_18 for each formula_78 such that the corresponding node formula_79, and formula_20 for each formula_78 such that the corresponding node formula_19.
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 and jump-flood, that can also take advantage of hardware acceleration in GPGPU implementations.
Functions of discrete variables with more than two values.
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
where formula_84 and formula_85. The function formula_86 represents the unary contribution of each variable (often referred as "data term"), while the function formula_87 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. 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 formula_87 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.
Given a function formula_89 with formula_90, and a certain assignment of values formula_91 to the variables, it is possible to associate each assignment formula_3 to a partition formula_93 of the set of variables, such that, formula_94. Give two distinct assignments formula_95 and formula_96 and a value formula_97, a move that transforms formula_95 into formula_96 is said to be an formula_100-expansion if formula_101 and formula_102. Given a couple of values formula_100 and formula_104, a move is said to be an formula_105-swap if formula_106. Intuitively, an formula_100-expansion move from formula_3 assigns the value of formula_100 to some variables that have a different value in formula_3, while an formula_105-swap move assigns formula_100 to some variables that have value formula_104 in formula_3 and vice versa.
For each iteration, the formula_100-expansion algorithm computes, for each possible value formula_100, the minimum of the function among all assignments formula_117 that can be reached with a single formula_100-expansion move from the current temporary solution formula_3, and takes it as the new temporary solution.
 formula_120
 formula_121
 while formula_122:
 formula_123
 foreach formula_97:
 formula_125
 if formula_126:
 formula_127
 formula_121
The formula_105-swap algorithm is similar, but it searches for the minimum among all assignments formula_130 reachable with a single formula_105-swap move from formula_3.
 formula_120
 formula_121
 while formula_122:
 formula_123
 foreach formula_137:
 formula_138
 if formula_126:
 formula_127
 formula_121
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.
The solution generated by such algorithms is not necessarily a global optimum, but it has strong guarantees of optimality. If formula_87 is a metric and formula_3 is a solution generated by the formula_100-expansion algorithm, or if formula_87 is a semimetric and formula_3 is a solution generated by the formula_105-swap algorithm, then formula_6 lies within a known and constant factor from the global minimum formula_149:
Non-submodular functions.
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.
In case of quadratic non-submodular functions, it is possible to compute in polynomial time a partial solution using algorithms such as QPBO. Higher-order functions can be reduced in polynomial time to a quadratic form that can be optimised with QPBO.
Higher-order functions.
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.
Sufficient conditions analogous to submodularity were developed to characterise higher-order pseudo-Boolean functions that can be optimised in polynomial time, and there exists algorithms analogous to formula_100-expansion and formula_105-swap for some families of higher-order functions. 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.