Reflection methods for user-friendly submodular optimization

Introduction

Submodular functions underlie the goals of numerous problems in machine learning, computer vision and signal processing . Several problems in these areas can be phrased as submodular optimization tasks: notable examples include graph cut-based image segmentation , sensor placement , or document summarization . A longer list of examples may be found in .

The theoretical complexity of submodular optimization is well-understood: unconstrained minimization of submodular set functions is polynomial-time while submodular maximization is NP-hard. Algorithmically, however, the picture is different. Generic submodular maximization admits efficient algorithms that can attain approximate optima with global guarantees; these algorithms are typically based on local search techniques . In contrast, although polynomial-time solvable, submodular function minimization (SFM) which seeks to solve

poses substantial algorithmic difficulties. This is partly due to the fact that one is commonly interested in an exact solution (or an arbitrarily close approximation thereof), and “polynomial-time” is not necessarily equivalent to “practically fast”.

Submodular minimization algorithms may be obtained from two main perspectives: combinatorial and continuous. Combinatorial algorithms for SFM typically use close connections to matroid and maximum flow methods; the currently theoretically fastest combinatorial algorithm for SFM scales as $O(n^{6}+n^{5}\tau)$, where $\tau$ is the time to evaluate the function oracle (for an overview of other algorithms, see e.g., ). These combinatorial algorithms are typically nontrivial to implement.

Continuous methods offer an alternative by instead minimizing a convex extension. This idea exploits the fundamental connection between a submodular function $F$ and its Lovász extension $f$ , which is continuous and convex. The SFM problem (1) is then equivalent to

The Lovász extension $f$ is nonsmooth, so we might have to resort to subgradient methods. While a fundamental result of Edmonds demonstrates that a subgradient of $f$ can be computed in $O(n\log n)$ time, subgradient methods can be sensitive to choices of the step size, and can be slow. They theoretically converge at a rate of $O(1/\sqrt{t})$ (after $t$ iterations). The “smoothing technique” of does not in general apply here because computing a smoothed gradient is equivalent to solving the submodular minimization problem. We discuss this issue further in Section 2.

An alternative to minimizing the Lovász extension directly on $^{n}$ is to consider a slightly modified convex problem. Specifically, the exact solution of the discrete problem $\min_{S\subseteq V}F(S)$ and of its nonsmooth convex relaxation $\min_{x\in^{n}}f(x)$ may be found as a level set $S_{0}=\{k\mid x^{\ast}_{k}\geqslant 0\}$ of the unique point $x^{*}$ that minimizes the strongly convex function :

We will refer to the minimization of (3) as the proximal problem due to its close similarity to proximity operators used in convex optimization . When $F$ is a cut function, (3) becomes a total variation problem (see, e.g., and references therein) that also occurs in other regularization problems . Two noteworthy points about (3) are: (i) addition of the strongly convex component $\tfrac{1}{2}\|x\|^{2}$; (ii) the ensuing removal of the box-constraints $x\in^{n}$. These changes allow us to consider a convex dual which is amenable to smooth optimization techniques.

Typical approaches to generic SFM include Frank-Wolfe methods that have cheap iterations and $O(1/t)$ convergence, but can be quite slow in practice (Section 5); or the minimum-norm-point/Fujishige-Wolfe algorithm that has expensive iterations but finite convergence. Other recent methods are approximate . In contrast to several iterative methods based on convex relaxations, we seek to obtain exact discrete solutions.

To the best of our knowledge, all generic algorithms that use only submodularity are several orders of magnitude slower than specialized algorithms when they exist (e.g., for graph cuts). However, the submodular function is not always generic and given via a black-box, but has known structure. Following , we make the assumption that $F(S)=\sum_{i=1}^{r}F_{i}(S)$ is a sum of sufficiently “simple” functions (see Sec. 3). This structure allows the use of (parallelizable) dual decomposition techniques for the problem in Eq. (2), with or without Nesterov’s smoothing technique, or with direct smoothing techniques. But existing approaches typically have two drawbacks: (a) they use smoothing or step-size parameters whose selection may be critical and quite tedious; and (b) they still exhibit slow convergence (see Section 5).

These drawbacks arise from working with formulation (2). Our main insight is that, despite seemingly counter-intuitive, the proximal problem (3) offers a much more user-friendly tool for solving (1) than its natural convex counterpart (2), both in implementation and running time. We approach Problem (3) via its dual. This allows decomposition techniques which combine well with orthogonal projection and reflection methods that (a) exhibit faster convergence, (b) are easily parallelizable, (c) require no extra hyperparameters, and (d) are extremely easy to implement.

The main three algorithms that we consider are: (i) dual block-coordinate descent (equivalently, primal-dual proximal-Dykstra), which was already shown to be extremely efficient for total variation problems that are special cases of Problem (3); (ii) Douglas-Rachford splitting using the careful variant of , which for our formulation (Section 4.2) requires no hyper-parameters; and (iii) accelerated projected gradient . We will see these alternative algorithms can offer speedups beyond known efficiencies. Our observations have two implications: first, from the viewpoint of solving Problem (3), they offers speedups for often occurring denoising and reconstruction problems that employ total variation. Second, our experiments suggest that projection and reflection methods can work very well for solving the combinatorial problem (1).

In summary, we make the following contributions:

Review of relevant results from submodular analysis

The relevant concepts we review here are the Lovász extension, base polytopes of submodular functions, and relationships between proximal and discrete problems. For more details, see .

The power set $2^{V}$ may be naturally identified with the vertices of the hypercube, i.e., $\{0,1\}^{n}$. The Lovász extension $f$ of any set function is defined by linear interpolation, so that for any $S\subset V$, $F(S)=f(1_{S})$. It may be computed in closed form once the components of $x$ are sorted: if $x_{\sigma(1)}\geqslant\cdots\geqslant x_{\sigma(n)}$, then $f(x)=\sum_{k=1}^{n}x_{\sigma(k)}\big[F(\{\sigma(1),\dots,\sigma(k)\})-F(\{\sigma(1),\dots,\sigma(k-1)\})\big]$ . For the graph cut function, $f$ is the total variation.

In this paper, we are going to use two important results: (a) if the set function $F$ is submodular, then its Lovász extension $f$ is convex, and (b) minimizing the set function $F$ is equivalent to minimizing $f(x)$ with respect to $x\in^{n}$. Given $x\in^{n}$, all of its level sets may be considered and the function may be evaluated (at most $n$ times) to obtain a set $S$. Moreover, for a submodular function, the Lovász extension happens to be the support function of the base polytope $B(F)$ defined as

that is $f(x)=\max_{y\in B(F)}y^{\top}x$ . A maximizer of $y^{\top}x$ (and hence the value of $f(x)$), may be computed by the “greedy algorithm”, which first sorts the components of $w$ in decreasing order $x_{\sigma(1)}\geqslant\cdots\geqslant x_{\sigma(n)}$, and then compute $y_{\sigma(k)}=F(\{\sigma(1),\dots,\sigma(k)\})-F(\{\sigma(1),\dots,\sigma(k-1)\})$. In other words, a linear function can be maximized over $B(F)$ in time $O(n\log n+n\tau)$ (note that the term $n\tau$ may be improved in many special cases). This is crucial for exploiting convex duality.

Dual of discrete problem.

We may derive a dual problem to the discrete problem in Eq. (1) and the convex nonsmooth problem in Eq. (2), as follows:

where $(y)_{-}=\min\{y,0\}$ applied elementwise. This allows to obtain dual certificates of optimality from any $y\in B(F)$ and $x\in^{n}$.

Proximal problem.

The dual problem of Problem (3) reads as follows:

where primal and dual variables are linked as $x=-y$. Observe that this dual problem is equivalent to finding the orthogonal projection of $0$ onto $B(F)$.

Divide-and-conquer strategies for the proximal problems.

Decomposition of submodular functions

The key to the algorithms presented here is to be able to minimize $\tfrac{1}{2}\|x-z\|_{2}^{2}+f_{j}(x)$, or equivalently, to orthogonally project $z$ onto $B(F_{j})$: $\min\tfrac{1}{2}\|y-z\|_{2}^{2}$ subject to $y\in B(F_{j})$.

We next sketch some examples of functions $F$ and their decompositions into simple functions $F_{j}$. As shown at the end of Section 2, projecting onto $B(F_{j})$ is easy as soon as the corresponding submodular minimization problems are easy. Here we outline some cases for which specialized fast algorithms are known.

where $g_{j}(\lambda_{j})=\min_{S\subset V}F_{j}(S)-\lambda_{j}(S)$ is a nonsmooth concave function.

The dual is the maximization of a nonsmooth concave function over a convex set, onto which it is easy to project: the projection of a vector $y$ has $j$-th block equal to $y_{j}-\frac{1}{r}\sum_{k=1}^{r}y_{k}$. Moreover, in our setup, functions $g_{j}$ and their subgradients may be computed efficiently through SFM.

We consider several existing alternatives for the minimization of $f(x)$ on $x\in^{n}$, most of which use Lemma 1. Computing subgradients for any $f_{j}$ means calling the greedy algorithm, which runs in time $O(n\log n)$. All of the following algorithms require the tuning of an appropriate step size.

Dual smoothing with entropy also admits coordinate descent methods that exploit the decomposition, but we do not compare to those here.

2 Dual decomposition methods for proximal problems

We may also consider Eq. (3) and first derive a dual problem using the same technique as in Section 3.1. Lemma 2 (proved in Appendix A) formally presents our dual formulation as a best approximation problem. The primal variable can be recovered as $x=-\sum_{j}y_{j}$.

The dual of Eq. (3) may be written as the best approximation problem

We can actually eliminate the $\lambda_{j}$ and obtain the simpler looking dual problem

Such a dual was also used in . In Section 5, we will see the effect of solving one of these duals or the other. For the simpler dual (8) the case $r=2$ is of special interest; it reads

We write Problem (9) in this suggestive form to highlight its key geometric structure: it is, like (7), a best approximation problem: i.e., the problem of finding the closest point between the polytopes $B(F_{1})$ and $-B(F_{2})$. Notice, however, that (7) is very different from (9)—the former operates in a product space while the latter does not, a difference that can have impact in practice (see Section 5). We are now ready to present algorithms that exploit our dual formulations.

Algorithms

We describe a few competing methods for solving our smooth dual formulations. We describe the details for the special 2-block case (9); the same arguments apply to the block dual from Lemma 2.

Perhaps the simplest approach to solving (9) (viewed as a minimization problem) is to use a block coordinate descent (BCD) procedure, which in this case performs the alternating projections:

The iterations for solving (8) are analogous. This BCD method (applied to (9)) is equivalent to applying the so-called proximal-Dykstra method to the primal problem. This may be seen by comparing the iterates. Notice that the BCD iteration (10) is nothing but alternating projections onto the convex polyhedra $B(F_{1})$ and $B(F_{2})$. There exists a large body of literature studying method of alternating projections—we refer the interested reader to the monograph for further details.

However, despite its attractive simplicity, it is known that BCD (in its alternating projections form), can converge arbitrarily slowly depending on the relative orientation of the convex sets onto which one projects. Thus, we turn to a potentially more effective method.

2 Douglas-Rachford splitting

The Douglas-Rachford (DR) splitting method includes algorithms like ADMM as a special case . It avoids the slowdowns alluded to above by replacing alternating projections with alternating “reflections”. Formally, DR applies to convex problems of the form

subject to the qualification $\relint(\dom\phi_{1})\cap\relint(\dom\phi_{2})\neq\varnothing$. To solve (11), DR starts with some $z_{0}$, and performs the three-step iteration (for $k\geq 0$):

where $\gamma_{k}\in$ is a sequence of scalars that satisfy $\sum\nolimits_{k}\gamma_{k}(2-\gamma_{k})=\infty$. The sequence $\{x_{k}\}$ produced by iteration (12) can be shown to converge to a solution of (11) [4; Thm. 25.6].

and setting $\gamma_{k}=1$, the DR iteration (12) may be written in a more symmetric form as

Applying DR to the duals (7) or (9), requires first putting them in the form (11), either by introducing extra variables or by going back to the primal, which is unnecessary. This is where the special structure of our dual problem proves crucial, a recognition that is subtle yet remarkably important.

Instead of applying DR to (9), consider the closely related problem

where $\delta_{1}$, $\delta_{2}^{-}$ are indicator functions for $B(F_{1})$ and $-B(F_{2})$, respectively. Applying DR directly to (14) does not work because usually $\relint(\dom\delta_{1})\cap\relint(\dom\delta_{2})=\varnothing$. Indeed, applying DR to (14) generates iterates that diverge to infinity [5; Thm. 3.13(ii)]. Fortunately, even though the DR iterates for (14) may diverge, Bauschke et al. 2004 show how to extract convergent sequences from these iterates, which actually solve the corresponding best approximation problem; for us this is nothing but the dual (9) that we wanted to solve in the first place. Theorem 3, which is a simplified version of [5; Thm. 3.13], formalizes the above discussion.

Let $\mathcal{A}$ and $\mathcal{B}$ be nonempty polyhedral convex sets. Let $\Pi_{\mathcal{A}}$ ($\Pi_{\mathcal{B}}$) denote orthogonal projection onto $\mathcal{A}$ ($\mathcal{B}$), and let $R_{\mathcal{A}}:=2\Pi_{\mathcal{A}}-\id$ (similarly $R_{\mathcal{B}}$) be the corresponding reflection operator. Let $\{z_{k}\}$ be the sequence generated by the DR method (13) applied to (14). If $\mathcal{A}\cap\mathcal{B}\neq\varnothing$, then $\{z_{k}\}_{k\geq 0}$ converges weakly to a fixed-point of the operator $T:=\tfrac{1}{2}[R_{\mathcal{A}}R_{\mathcal{B}}+\id]$; otherwise $\|{z_{k}}\|_{2}\to\infty$. The sequences $\{x_{k}\}$ and $\{\Pi_{\mathcal{A}}\Pi_{\mathcal{B}}z_{k}\}$ are bounded; the weak cluster points of either of the two sequences

are solutions best approximation problem $\min_{a,b}\|a-b\|$ such that $a\in\mathcal{A}$ and $b\in\mathcal{B}$.

The key consequence of Theorem 3 is that we can apply DR with impunity to (14), and extract from its iterates the optimal solution to problem (9) (from which recovering the primal is trivial). The most important feature of solving the dual (9) in this way is that absolutely no stepsize tuning is required, making the method very practical and user friendly (see also Appendix D).

Experiments

We empirically compare the proposed projection methods Code and data corresponding to this paper are available at https://sites.google.com/site/mloptstat/drsubmod to the (smoothed) subgradient methods discussed in Section 3.1. For solving the proximal problem, we apply block coordinate descent (BCD) and Douglas-Rachford (DR) to Problem (8) if applicable, and also to (7) (BCD-para, DR-para). In addition, we use acceleration to solve (8) or (9) . The main iteration cost of all methods except for the primal subgradient method is the orthogonal projection onto polytopes $B(F_{j})$, and therefore the number of iterations is a suitable criterion for comparisons. The primal subgradient method uses the greedy algorithm in each iteration, which runs in $O(n\log n)$. However, as we will see, its convergence is so slow to counteract any benefit that may arise from not using projections. We do not include Frank-Wolfe methods here, since FW is equivalent to a subgradient descent on the primal and converges correspondingly slowly.

As benchmark problems, we use (i) graph cut problems for segmentation, or MAP inference in a 4-neighborhood grid-structured MRF, and (ii) concave functions similar to those used in , but together with graph cut functions. The segmentation problems (i) are set up in a fairly standard way on a 4-neighbor grid graph, with unary potentials derived from Gaussian Mixture Models of color features. The weight of graph edge $(i,j)$ is a function of $\exp(-\|y_{i}-y_{j}\|^{2})$, where $y_{i}$ is the RGB color vector of pixel $i$. The functions in (i) decompose as sums over vertical and horizontal paths. All horizontal paths are independent and can be solved together in parallel, and similarly all vertical paths. The functions in (ii) are constructed by extracting regions $R_{j}$ via superpixels and, for each $R_{j}$, defining the function $F_{j}(S)=|S||R_{j}\setminus S|$. We use 200 and 500 regions. The problems have size $640\times 427$. Hence, for (i) we have $r=640+427$ (but solve it as $r=2$) and for (ii) $r=640+427+500$ (solved as $r=3$).

For algorithms working with formulation (7), we compute an improved smooth duality gap of a current primary solution $x=-\sum_{j}y_{j}$ as follows: find $y^{\prime}\in\argmax_{y\in B(F)}x^{\top}y$ (then $f(x)=x^{\top}y^{\prime}$) and find an improved $x^{\prime}$ by minimizing $\min_{z}z^{\top}y^{\prime}+\tfrac{1}{2}\|z\|^{2}$ subject to the constraint that $z$ has the same ordering as $x$ . The constraint ensures that $(x^{\prime})^{\top}y^{\prime}=f(x^{\prime})$. This is an isotonic regression problem and can be solved in time $O(n)$ using the “pool adjacent violators” algorithm . The gap is then $f(x^{\prime})+\tfrac{1}{2}\|x^{\prime}\|^{2}-(-\tfrac{1}{2}\|y^{\prime}\|^{2})$.

For computing the discrete gap, we find the best level set $S_{i}$ of $x$ and, using $y^{\prime}=-x$, compute $\min_{i}F(S_{i})-y^{\prime}_{-}(V)$.

Two functions ($r=2$). Figure 2 shows the duality gaps for the discrete and smooth (where applicable) problems for two instances of segmentation problems. The algorithms working with the proximal problems are much faster than the ones directly solving the nonsmooth problem. In particular DR converges extremely fast, faster even than BCD which is known to be a state-of-the-art algorithms for this problem . This, in itself, is a new insight for solving TV. We also see that the discrete gap shrinks faster than the smooth gap, i.e., the optimal discrete solution does not require to solve the smooth problem to extremely high accuracy. Figure 1 illustrates example results for different gaps.

More functions ($r>2$). Figure 3 shows example results for four problems of sums of concave and cut functions. Here, we can only run DR-para. Overall, BCD, DR-para and the accelerated gradient method perform very well.

If we aim for parallel methods, then again DR outperforms BCD. Figure 4 (right) shows the speedup gained from parallel processing for $r=2$. Using 8 cores, we obtain a 5-fold speed-up.

Running time compared to graph cuts

Table 1 shows the running times of our DR method (implemented in Matlab/C++) and the Maxflow code of (using the wrapper ) for the four graph cut (segmentation) instances above on a MacBook Air with a 2 GHz Intel Core i7. The running times are averages over 5 repetitions. DR was run for 10, 10, 21, and 20 iterations, respectively.

DR is by a factor of 2-9 slower than the specialized code. Given that, as opposed to the combinatorial algorithm, DR solves the full regularization path, is parallelizable, generic and straightforwardly extends to a variety of functions, this is remarkable.

In summary, our experiments suggest that projection methods can be extremely useful for solving the combinatorial submodular minimization problem. Of the tested methods, DR, cyclic BCD and accelerated gradient perform very well. For parallelism, applying DR on (9) converges much faster than BCD on the same problem.

Conclusion

We have presented a novel approach to submodular function minimization based on the equivalence with a best approximation problem. The use of reflection methods avoids any hyperparameters and reduce the number of iterations significantly, suggesting the suitability of reflection methods for combinatorial problems. Given the natural parallelization abilities of our approach, it would be interesting to perform detailed empirical comparisons with existing parallel implementations of graph cuts (e.g., ). Moreover, a generalization beyond submodular functions of the relationships between combinatorial optimization problems and convex problems would enable the application of our framework to other common situations such as multiple labels (see, e.g., ).

This research was in part funded by the Office of Naval Research under contract/grant number N00014-11-1-0688, by NSF CISE Expeditions award CCF-1139158, by DARPA XData Award FA8750-12-2-0331, and the European Research Council (SIERRA project), as well as gifts from Amazon Web Services, Google, SAP, Blue Goji, Cisco, Clearstory Data, Cloudera, Ericsson, Facebook, General Electric, Hortonworks, Intel, Microsoft, NetApp, Oracle, Samsung, Splunk, VMware and Yahoo!. We would like to thank Martin Jaggi, Simon Lacoste-Julien and Mark Schmidt for discussions.

References

Appendix A Derivations of Dual Problems

To derive the non-smooth dual problem, we follow and use Lagrangian duality:

where $g_{j}(\lambda_{j})=\min_{A\subset V}F_{j}(A)-\lambda_{j}(A)$ is a nonsmooth concave function, which may be computed efficiently through submodular function minimization. ∎

A.2 Proof of Lemma 2

The proof follows a similar saddle-point approach.

Writing (16) as a minimization problem and ignoring constants completes the proof. ∎

Appendix B Divide-and-conquer algorithm for parametric submodular minimization

Here, we extend the approach of Tarjan et al. 2006 for parametric max-flow to all submodular functions and all monotone strictly convex functions beyond the square functions used in the main paper. More precisely, we consider a submodular function $F$ defined on $V=\{1,\dots,n\}$ and $n$ differentiable strictly convex functions $h_{i}$ such that their Fenchel-conjugates $h_{i}^{\ast}$ have full domain, for $i\in\{1,\dots,n\}$. The functions $h^{\ast}_{i}$ are then differentiable. We consider the following problem:

$-y_{i}=h_{i}^{\prime}(x_{i})\Leftrightarrow x_{i}=(h_{i}^{\ast})^{\prime}(-y_{i})$.

Algorithm 1 is a divide-and-conquer algorithm. In each recursive call, it takes an interval $[\lambda_{\min},\lambda_{\max}]$ in which all components of the optimal solution lie and either (a) shortens the search interval for any break point, (b) finds the optimal (constant) value of $x$ on a range of elements, or (c) recursively splits the problem into a set $S$ and $V\setminus S$ with corresponding ranges for the values of $x^{*}$ and finds the optimal values of $x$ on the two subsets.

B.2 Review of related results

The goal of this appendix is to show Proposition 4 below. We first start by reviewing existing results regarding separable problems on the base polyhedron (see for details).

It is known that if $y\in B(F)$, then $y_{k}\in\big[F(V)-F(V\backslash\{k\}),F(\{k\})\big]$; thus, the optimal solution $x$ is such that $x_{k}\in\big[(h_{k}^{\ast})^{\prime}(-F(\{k\})),(h_{k}^{\ast})^{\prime}(F(V\backslash\{k\})-F(V))\big].$ We therefore set the initial search range to

Let $\alpha<\beta$ and $S^{\alpha}$ be any minimizer of $F(S)+h^{\prime}(\alpha)(S)$ and $S^{\beta}$ any minimizer of $F(S)+h^{\prime}(\beta)(S)$. Then $S^{\beta}\subseteq S^{\alpha}$.

The coordinates $x_{j}^{*}$ ($j\in V$) of the unique optimal solution $x^{*}$ of Problem 17 are

where $S^{\lambda}$ is any minimizer of $F(S)+h^{\prime}(\lambda)(S)$.

Propositions 1 and 2 imply that the level sets of $x^{*}$ form a chain $\emptyset=S_{0}\subset S_{1}\subset\ldots\subset S_{k}=V$ of maximal minimizers for the critical values of $\lambda$ (which are the entries of $x^{*}$). (Each $S_{i}=S^{\lambda}$ for some $\lambda=x^{*}_{j}$.)

Then $x^{*}_{j}=y_{j}$ for $j\in T$ and $x^{*}_{j}=z_{j}$ for $j\in V\setminus T$.

The algorithm uses Proposition 3 recursively.

Let $\lambda$ be the value in $x^{*}$ defining $S_{i}=S^{\lambda}$. It is easy to see that the restriction $F_{T}$ and the contraction $F^{T}$ are both submodular. Hence, Propositions 1 and 2 hold for them.

Since the restriction on $T$ is equivalent to the original function for any $S\subseteq T$, $F(S)+h(\lambda)(S)=F_{T}(S)+h_{T}(\lambda)(S)$ for any $S\subseteq T$. With this, Propositions 1 and 2 imply that for any $\alpha>\lambda$, $F(S)+h(\lambda)(S)=F_{T}(S)+h_{T}(\lambda)(S)$ and therefore $x^{*}_{j}=y_{j}$ for $j\in T$.

Similarly, for any $S\in V\setminus T$, it holds that $F(S\cup T)+h(\lambda)(S\cup T)=F^{T}(S)+(h^{\prime}(\lambda))^{T}(S)+F(T)+h^{\prime}(\lambda)(T)$. Due to the monotonicity property of the optimizing sets, $S^{\alpha}\supseteq T$ for all $\alpha<\lambda$, and therefore the maximal minimizer $U^{\alpha}$ of $F^{T}(S)+(h^{\prime}(\alpha))^{T}(S)$ satisfies $U^{\alpha}\cup T=S^{\alpha}$ (the terms $F(T)+h^{\prime}(\lambda)(T)$ are constant with respect to $U$). Hence Poposition 2 implies that $x^{*}_{j}=z_{j}$ for $j\in V\setminus T$. ∎

These propositions imply that there is a set of at most $n$ values of $\lambda=\alpha$ that define the level sets $S^{\alpha}$ of the optimal solution $x^{*}$. If we know these break point values, then we know $x^{*}$. Algorithm 1 interleaves an unbalanced split strategy that may split the seach interval in an unbalanced way but converges in $O(n)$ recursive calls, and a balanced split strategy that always halves the search intervals but is not finitely convergent.

B.3 Proof of convergence

We now prove the convergence rate for Algorithm 1.

The minimum of $f(x)+\sum_{i=1}^{n}h_{i}(x_{i})$ may be obtained up to coordinate-wise accuracy $\epsilon$ within

The proof relies on Propositions 1, 2 and 3.

The choice of $\lambda$ in the unbalanced splitting strategy corresponds to solving a simplified version of the dual problem. Indeed, by convex duality, the following two problems are dual to each other:

Problem (21) replaces the constraint that $y\in B(F)$ by $y(V)=F(V)$, dropping the constraint that $y(S)\leq F(S)$ for all $S\subseteq V$. Testing whether $y$ satisfies all constraints of (19), i.e., $y\in B(F)$ is equivalent to testing whether $F(S)-y(S)\geq 0$. We do this implicitly by our choice of $\lambda$: Convex duality implies that the the optimal solutions of Problems (21) and (22) satisfy $y_{i}=-h^{\prime}_{i}(\lambda)$. This holds in particular for the chosen (unique optimal) $\lambda$ in the algorithm.

Let $T$ be a minimizer of $F(S)+h^{\prime}(\lambda)(S)=F(S)-y(S)$. If $T=\emptyset$ or $T=V$, then $y\in B(F)$ and an optimal solution for the full dual problem (19). Hence, $y$ and $x=\lambda\mathbf{1}_{V}=(h^{*})^{\prime}(-y)$ form a primal/dual optimal pair for (19).

If $\emptyset\subset T\subset V$ and $F(T)-y(T)<0$, then $y\notin B(F)$, and we perform a split with the same argumentation as above. This splitting strategy is exactly that of and splits at most $n$ times. Hence, this strategy yields the global optimum (to machine precision) in the time of $O(n)$ times solving a submodular minimization on $V$. If $n$ is large, this may be computationally expensive.

If we only do balanced splits, we end up approaching the break points more and more closely (but typically never exactly). Unbalanced splits always find an exact break point, but with potentially little progress in reducing the intervals. Algorithm 1 thus interleaves both strategies where we store intervals of allowed values for subsets of components of $A$. At step $d$ there are at most $\min\{n,2^{d}\}$ different intervals (as there cannot be more intervals than elements of $V$). To split these intervals, submodular function minimization problems have to be solved on each of these intervals, with total complexity less than a single submodular function optimization problem on the full set. At each iteration, intervals corresponding to a singleton may be trivially completely solved, and components which are already found are discarded. Hence, at each recursive level, the total computation time is bounded above by $\tau(V)$.

Finally, we adress the precision for the special case that $h_{i}(x_{i})=\tfrac{1}{2}x_{i}^{2}$ for all $i\in V$. If the interval lengths are smaller than the smallest gap between any two break points (components of $x^{*}$), then each interval contains at most one break point and the algorithm converges after at most two unbalanced splits. Hence, we here consider $\epsilon$ to be one half times the smallest gap between any two break points. Let $\emptyset=S_{0}\subset S_{1}\subset\ldots\subset S_{k}=V$ be the chain of level sets of $x^{*}$. By the optimality conditions discussed above for unbalanced splits, any constant part $T=S_{i}\setminus S_{i-1}$ of $x^{*}$ takes value $\lambda\mathbf{1}=-y_{j}\mathbf{1}$ ($j\in T$), where $y(T)=F_{S_{i-1}}(T)$, and hence

Therefore, the (absolute) difference between any two such values is loosely lower bounded by

Note that in the case of flows, the algorithm is not exactly equivalent to the flow algorithm of , which updates flows directly.

Appendix C BCD and proximal Dykstra

We consider the best approximation problem

Let us show the details for only the two block case. The general case follows similarly.

Clearly, this problem contains the two-block best approximation problem as a special case (by setting $f$ and $h$ to be suitable indicator functions). Now introduce two variables $z,w$ that equal $x$; then the corresponding Lagrangian is

From this Lagrangian, a brief calculation yields the dual optimization problem

We solve this dual problem via BCD, which has the updates

Thus, $0\in\nu_{k+1}+\mu_{k}-y+\partial f^{*}(\nu_{k+1})$ and $0\in\nu_{k+1}+\mu_{k+1}-y+\partial h^{*}(\mu_{k+1})$. The first optimality condition may be rewritten as

Similarly, we second condition yields $\mu_{k+1}=y-\nu_{k+1}-\prox_{h}(y-\nu_{k+1})$. Now use Lagrangian stationarity

to rewrite BCD using primal and dual variables to obtain the so-called proximal-Dykstra method:

We discussed the more general problem (25) because it contains the smoothed primal as a special case, namely with $y=0$ in (25), $f=f_{1}$, and $h=f_{2}$, we obtain

for which BCD yields the proximal-Dykstra method that was previously used in for two-dimensional TV optimization.

Appendix D Recipe: Submodular minimization via reflections

To be precise, we summarize here how to solve Problem (17) via reflections. As we showed above, the dual is of the form

The vector $y$ consists of $r$ parts $y_{j}\in B(F_{j})$. We first solve the dual by starting with any $z^{(0)}\in\mathcal{H}^{r}$, and iterate

Upon convergence to a point $z^{*}$, we extract the components