An Algorithmic Theory of Dependent Regularizers, Part 1: Submodular Structure

We present an exploration of the rich theoretical connections between several classes of regularized models, network flows, and recent results in submodular function theory. This work unifies key aspects of these problems under a common theory, leading to novel methods for working with several important models of interest in statistics, machine learning and computer vision. In Part 1, we review the concepts of network flows and submodular function optimization theory foundational to our results. We then examine the connections between network flows and the minimum-norm algorithm from submodular optimization, extending and improving several current results. This leads to a concise representation of the structure of a large class of pairwise regularized models important in machine learning, statistics and computer vision. In Part 2, we describe the full regularization path of a class of penalized regression problems with dependent variables that includes the graph-guided LASSO and total variation constrained models. This description also motivates a practical algorithm. This allows us to efficiently find the regularization path of the discretized version of TV penalized models. Ultimately, our new algorithms scale up to high-dimensional problems with millions of variables.


Introduction
High-dimensional data is a central focus of modern research in statistics and machine learning. Recent technological advances in a variety of scientific fields for gathering and generating data, matched by rapidly increasing computing power for analysis, has attracted significant research into the statistical questions surrounding structured data. Numerous models and computational techniques have emerged recently to work with these data types.
Our primary contribution is a theoretical framework enabling the efficient optimization of models that work with high-dimensional models in which the predictors are believed to be dependent, correlated, or sparse. In particular, we propose a new approach for estimation of certain types of dependent predictors, in which the model incorporates a prior belief that many of the predictors take similar or identical values. This has been a hot topic of research in recent years, with applications in genomics, image analysis, graphical models, and several other areas. However, efficient estimation involving structured and dependent predictors has proven to be quite challenging. Our main contribution, which includes a number of related theoretical results, is a theoretical framework and algorithmic approach that unlocks a large and particularly thorny class of these models.
After laying out the context for our research in the next section, we outline our main results as well as the needed background in section 3. Our contributions build on recent results in optimization theory and combinatorial optimization, which we describe in detail. In 4, we present a clear theoretical connection between submodular optimization theory, network flows, and the proximal operator in regularized models such as the graph-guided LASSO. Then, in 5, we extend the underlying submodular optimization theory to allow a weighted version of the underlying submodular problem. In the context of regularized regression, this allows the use of additional unary terms in the regularizer.

Overview of Statistical Models and Regularization
To begin, consider a simple regression model, in which we have N observations y = (y 1 , y 2 , ..., y N ) and n predictors u = (u 1 , u 2 , ..., u n ), with where ε 1 , ε 2 , ..., ε N are independent random noise vectors with E ε i = 0; typically, these are assumed to be i.i.d. Gaussian. In many high dimensional contexts, the estimation ofû may be problematic using classical methods. For example, n may be much larger than N , making the problem ill-posed as many possible values of u map to the same response. (n and N are used here instead of the more common p and n to be consistent with the optimization literature we connect this problem to.) Additionally, many predictors of interest may have a negligible or nearly identical affect on y; eliminating or grouping these predictors is then desirable. To handle this, a number of sophisticated approaches have been proposed to incorporate variable selection or aggregation into the statistical estimation problem. One common and well-studied approach is to add a penalty term, or regularizer, to the loglikelihood that enforces some prior belief about the structure of u [Bickel et al., 2006, Hastie et al., 2009. In this setup, estimating the predictor u involves finding the minimizer of a log-likelihood term or loss function L plus a regularization term Φ: where λ controls the strength of the regularization λΦ(u).
Many well studied and frequently used models fall into this context. For example, if L is the log-likelihood from a multivariate Gaussian distribution, one of the oldest regularization techniques is the L 2 -norm, which gives us the classic technique of ridge regression [Hoerl and Kennard, 1970], whereû ridge = argmin u∈R n y − Au 2 2 + λ u 2 2 . (3) In this example, the regularization term is typically used as an effective way of dealing with an ill-conditioned inverse problem in the least squares context or as a simple way of preventing u from over-fitting the model [Bishop, 2006]. In the Bayesian context, (2) often corresponds directly to finding the maximum a posteriori estimate in the classic likelihood-prior formulation, i.e. p(u | y) ∝ p(y | u)p(u ; λ, θ) ∝ e −L(u,y) × e −λΦ θ (u) .
where λ and θ are hyperparameters controlling the behavior of the prior distribution. In this formulation, the prior captures the belief encoded by the regularization term. For example, the ridge regression problem of (3) corresponds to using a standard multivariate Gaussian likelihood but assumes a 0-mean spherical Gaussian prior distribution over u.
More recently, the Least Absolute Shrinkage and Selection Operator (LASSO), is used frequently to promote sparsity in the resulting estimator [Tibshirani, 1996]. This approach uses the L 1 -norm as the regularization term, i.e. u lasso = argmin u∈R n y − Au 2 2 + λ u 1 The LASSO problem has been analyzed in detail as a method for simultaneous estimation and variable selection, as the L 1 penalty tends to set many of the predictors to 0. The consistency and theoretical properties of this model, for both estimation and variable selection, are well studied [Meinshausen and Yu, 2009, Bunea et al., 2007, Van De Geer, 2008, van de Geer, 2007, Bickel et al., 2006. Many variants of this problem have also been proposed and analyzed. These include the elastic net, in which Φ(u) = α u 1 + (1 − α) u 2 2 , with α ∈ [0, 1]; this generalizes both ridge regression and LASSO [Zou and Hastie, 2005]. Data-adaptive versions of these estimators are analyzed in [Zou, 2006, Zou and Hastie, 2005, Zou and Zhang, 2009]. In addition, there are many efficient algorithms to quickly solve these problems, which is one of the reasons they are used frequently in practice [Friedman et al., 2008a, Efron et al., 2004, Beck and Teboulle, 2009].

Graph Structured Dependencies
While there are other possible regularization strategies on individual parameters, of most interest to us is the incorporation of correlation structures and dependencies among the predictors into the prior distribution or regularization term. In many cases, this translates into penalizing the differences of predictors. This has recently been attracting significant algorithmic and theoretical interest as a way to effectively handle dependency structures in the data.
In our particular context, we are interested in the MAP estimate of models with Markov Random Field prior. The resulting log-linear models consist of a collection of unary and pairwise terms that capture the dependency structure of the problem. In particular, the priors we are interested in enforce similarity between neighboring variables, where "neighbor" is defined according to the graph structure of the MRF prior. The type of dependency we are looking at is determined by the The simplest form of model dependency in the parameters is captured by the Fused LASSO problem , which penalizes pairwise differences between terms in an ordered problem to be u * fused = argmin u∈R n = Au − y 2 2 + λ w 1 u 1 + w 2 where w 1 , w 2 ≥ 0 control the balance between the L 1 -norm controlling variable sparsity and the sum of ordered differences that effectively penalizes changepoints. This problem has gained some recent attention in the statistics community in the context of non-parametric regression [Dümbgen and Kovac, 2009, Cho and Fryzlewicz, 2011, Davies and Meise, 2008, group sparsity [Tibshirani et al., 2005, Bleakley andVert, 2011], and change-point detection [Bleakley and Vert, 2011]. From the optimization perspective, several algorithms to find the solution to this problem have been proposed; we refer the reader to Liu et al. [2010], Ye and Xie [2011], Bach et al. [2012] or Friedman et al. [2007] for discussions of this particular problem. The generalized problem, in which the pairwise terms are not required to be ordered, is of significant practical interest in both the statistics and machine learning communities. Here, the pairwise interactions are controlled with an arbitrary graph of weighted difference penalties: where w 1 ≥ 0 controls the sparsity of the individual predictors, and w ij ≥ 0 penalizes differences, typically between predictors known to be correlated. This problem has become known as the graph guided LASSO Yu [2013]. This model -often with w 1 = 0 -arises in the context of "roughness" penalized regression in which differences between neighboring regions are penalized [Belkin et al., 2004]. Besag, Green, Higdon, and Mengersen [1995] examined pairwise interactions on Markov random fields, which describe a prior distribution on the edges of an undirected graph; the proposed models were then solved with MCMC methods. Similar problems also arise in the estimation of sparse covariance matrices as outlined by Friedman, Hastie, and Tibshirani [2008b]; there, structure learning is the goal and the pairwise interaction terms are not formed explicitly.
Similarly, Kovac and Smith [2011] analyzes this model as an approach to nonparametric regression on a graph, although with the simpler case of A = I. He proposes an efficient active-region based procedure to solve the resulting optimization problem. In the genomics community, a version of (8) has been proposed as the Graph-Guided Fused LASSO problem by Kim, Sohn, and Xing [2009]. There, the pairwise interactions were estimated from the correlation structure in the quantitative traits. Their model is essentially identical to this one, though with w 1 = 0. There, the authors used a general quadratic solver to tackle the problem and showed good statistical performance in detecting genetic markers. The use of recent results in optimization, in particular proximal operator methods [Nesterov, 2007], have been proposed as a way of optimizing this model Chen, Kim, Lin, Carbonell, and Xing [2010], Chen, Lin, Kim, Carbonell, and Xing [2012], Bach, Jenatton, Mairal, and Obozinski [2012]. The resulting algorithms, discovered independently, are similar to the ones we propose, although we improve upon them in several ways.
Some recent approaches focus on approximate methods with provable bounds on the solution. Nesterov proposed a smooth approximation approach which can be applied to the graph-guided LASSO in Nesterov [2007]. More recently, Yu [2013] proposed a non-smooth approximation method that decomposes the summation in (8).
A very similar model to (8) is proposed in Sharma, Bondell, and Zhang [2013]; there, however, the authors include an additional regularization term penalizing |u i + u j |. In addition, they give several adaptive schemes for choosing the weights in a way that depends on the correlation structures of the design matrix. This work is noteworthy in that the authors give a detailed analysis of the asymptotic convergence and estimation rate of the model when the regularization weights are chosen adaptively. Their optimization method, however, involves a simple quadratic programming setup, which can severely limit the size of the problems for which their estimator can be used.
In the context of regularized regression with pairwise dependencies, our contribution is a theoretical treatment and several algorithms for a class of models that generalizes (8). Namely, we allow the |u i | penalty in the L 1 norm to be replaced with an arbitrary convex piecewise linear function ξ i (u i ). In particular, we examine the estimator which generalizes many of the above models. While our primary result is a representation of the theoretical structures underlying this optimization problem, our work motivates efficient and novel algorithms for working with these types of structure. The theory we develop leads to proofs of the correctness of our algorithms. In addition, we are able to present an algorithmic description that gives the entire regularization path for one version of this problem.

Total Variation Models
Total variation models are effectively an extension of the above pairwise regularization schemes to the estimation of continuous functions. In this context, we wish to estimate a function u in which the total variation of u -the integral of the norm of the gradient -is controlled. This model is used heavily in developing statistical models of images [Chambolle et al., 2010] and the estimation of density functions [Bardsley and Luttman, 2009], but arises in other contexts as well.
When dealing with total variation models, we work theoretically with general L 2 -measurable functions with bounded variation. Our response becomes a function f , and our predictor becomes a function u. However, in practice, we discretize the problem by working with the functions f and u only on a discrete lattice of points, as is common in these problems.
Let Ω ⊂ R d be a compact region, and let u be an L 1 -integrable continuously differentiable function defined on Ω. In this context, the total variation of u is given by More general definitions exist when u is not continuously differentiable Ambrosio and Di Marino [2012], Giusti [1984]; for simplicity, we assume this condition. Let F(Ω) represent the space of L 1integrable continuously differentiable functions of bounded variation. Here, bounded variation can be taken as the condition TV(u) < +∞.
Again, much more general definitions exist, but this suffices for our purposes. Formally, then, the total variation problem seeks to find an estimation function u * that to a  function f : R d ⊃ Ω → R under a constraint or penalty on the total norm of the gradient of u.
This model -often called the Rudin-Osher-Fatimi (ROF) model after the authors -was originally proposed in Rudin et al. [1992] as a way of removing noise from images corrupted by Gaussian white noise; however, it has gained significant attention in a variety of other areas.
In the image analysis community, it has been used often as a model for noise removal, with the solution of (10) being the estimation of the noise-free image. In addition, the value of the total variation term in the optimal solution is often used as a method of edge detection, as the locations of the non-zero gradients tend to track regions of significant change in the image.
Beyond image analysis, a number of recent models have used total variation regularization as a general way to enforce similar regions of the solution to have a common value. In this way, it is a generalization of the behavior of the Fused LASSO in one dimension; it tends to set regions of the image to constant values. As such, it is used in MRI reconstruction, electron tomography [Goris et al., 2012, Gao et al., 2010, MRI modeling [Michel et al., 2011, Keeling et al., 2012, CT reconstruction [Tian et al., 2011], and general ill-posed problems [Bardsley and Luttman, 2009].
The total variation regularizer tends to have a smoothing effect around regions of transition, while also setting similar regions to a constant value. As such, it has proven to be quite effective for removing noise from images and finding the boundaries of sufficiently distinct regions. The locations where the norm of the gradient TV(u) is nonzero correspond to boundaries of the observed process in which the change is the greatest. A rigorous treatment of the theoretical aspects surrounding using this as a regularization term can be found in Bellettini, Caselles, and Novaga [2002], Ring [2000] and Chan and Esedoglu [2005]. A complete and treatment of this topic in practice can be found in a number of surveys, including Darbon and Sigelle [2006], Allard [2007Allard [ , 2008Allard [ , 2009; in particular, see Caselles, Chambolle, and Novaga [2011] and Chambolle, Caselles, Cremers, Novaga, and Pock [2010]. As our purpose in this work is to present a treatment of the underlying theoretical structures, we refer the reader to one of the above references for use practical image analysis. Still, several points are of statistical interest, which we discuss now.
Several generalizations of the traditional Rudin-Osher-Fatimi model have been proposed for different models of the noise and assumptions about the boundaries of the images. Most of these involve changes to the loss or log-likelihood term but preserve the total variation term as the regularizer. In general, we can treat u * as the estimation of where L is a smooth convex loss function. One option is to use L 1 loss for L, namely as a way of making the noise more robust to outliers. This model is explored in Bect et al. [2004], Chan and Esedoglu [2005] and the survey papers above. These perform better under some types of noise and have also become popular [Goldfarb and Yin, 2009]. Similarly, many models for image denoising -often when the image comes from noisy sensing processes in physics, medical imaging, or astronomy -involve using a Poisson likelihood for the observations. Here, we wish to recover a density function where we observe Poisson counts in a number of cells with rate proportional to the true density of the underlying process. Total variation regularization can be employed to deal with low count rates. Here, total variation regularization is used to recover the major structures in the data. In this case, the optimal loss function is given by Research on this, including optimization techniques, can be found in Bardsley and Luttman [2009], Bardsley [2008], Sawatzky et al. [2009]. Also of interest in the statistics community, Wunderli [2013] replaces the squared-norm fit penalty in (13) with the quantile regression penalty from Koenker and Bassett Jr [1978], where loss function is replaced with for which the author proves the existence, uniqueness, and stability of the solution.
One of the variations of total variation minimization is the Mumford-Shah model. In this model, used primarily for segmenting images, regions of discontinuity are handled explicitly Mumford and Shah [1989]. In this model, the energy which is minimized is where Γ is a collection of curves and Γ is the total length of these curves. Thus the function is allowed to be discontinuous on Γ; these are then taken to be the segmentations of the image. This model has also attracted a lot of attention in the vision community as a way of finding boundaries in the image [Chan and Vese, 2000, Pock et al., 2009, El Zehiry et al., 2007. The Bayesian model of this is as a type of mixture model in which the mean function and noise variance differs between regions [Brox and Cremers, 2007]. Unlike general TV regularization, however, it is typically focused exclusively on finding curves for image segmentation, rather than on estimating the original image. While it shares a number of close parallels to our problem, particularly in the graph cut optimizations used [El Zehiry et al., 2007], our method does not appear to generalize to this problem. The optimization of the total variation problem has garnered a lot of attention as well. An enormous number of approaches have been proposed for this problem. These largely fall into the category of graph based algorithms, active contour methods, or techniques based around partial differential equations. These are discussed in more depth in part 2 of this paper, where we mention them in introducing our approach to the problem.

Background and Outline of Main Results
The primary purpose of our work is to present a theory of the underlying structure connecting the optimization of the above models over dependent variables and recent results in combinatorial optimization, particularly submodular function minimization. Our theory connects and makes explicit a number of connections between known results in these fields, many of which we extend in practically relevant ways. The insights gained from our theory motivate a family of novel algorithms for working with these models. Furthermore, we are able to give a complete description of the structure of the regularization path; this result is also completely novel.
One primary practical contribution is the development of methods to efficiently optimize functions of the following form, which we denote as ( B ) with a ∈ R n , λ is a non-negative regularization parameter, w ij ≥ 0 controls the regularization of dependent variables, and ξ i (u i ) is a convex piecewise linear function, possibly the common L 1 penalty ξ i (u i ) = |u i |.
We develop a novel algorithm for ( B ) based on network flows and submodular optimization theory, and completely solve the regularization path as well over λ > 0. The crux of the idea is to construct a collection of graph-based binary partitioning problems indexed by a continuous parameter. With our construction, we show that the points at which each node exactly gives the solution to ( B ). We develop algorithms for this that scale easily to millions of variables, and show that the above construction also unlocks the nature of the regularization path.
These functions arise in two active areas of interest to the statistics and machine learning communities. The first is in the optimization of penalized regression functions where we wish to find where y is a collection of observed response and L(u, y) is a smooth, convex log-likelihood term. Throughout our work, we denote this problem as ( L ). In the general case, we handle this problem through the use of the proximal gradient methods, a variant of sub-gradient methods, which we described below. The inner routine involves solving ( B ).
Finally, ( B ) also extends to the continuous variational estimation problem of Total Variation minimization from section 2.2. where we wish to minimize where F(Ω) and TV(u) are defined in (10) -(13). We show in the second part of our work that this problem can be solved using the methods developed for ( B ). Thus we not only present an efficient algorithm for ( TV ), we also present the first algorithm for finding the full regularization path of this problem.

Outline
This paper is laid out as followed. The rest of this section presents a survey of the background to the current problem, describing the current state of the research in the relevant areas and laying out the reasons it is of interest to the statistics and machine learning communities. Section 4 lays out the basic framework connecting network flows, submodular function minimization, and the optimization of ( B ). While some of our results were discovered independently, our theory pulls them together explicitly in a common framework, and provides more direct proofs of their correctness. Additionally, our approach motivates a new algorithm to exactly solve the resulting optimization problem.
Section 5 presents entirely novel results. We extend the basic theory of section 4 using a novel extension lemma to allow for a general size-biasing measure on the different components of the input. The theoretical results are proved for general submodular functions, extending the state-of-the-art in this area. In the context of our problem, we use these results to extend the algorithmic results of section 4. In particular, this opens a way to handle the ξ i (u i ) term in ( B ) above, significantly extending the results of  related to section 4.

Optimization of Regularized Models
Regularized models have gained significant attention in both the statistics and machine learning communities. Intuitively, they provide a way of imposing structure on the problem solution. This facilitates the accurate and tractable estimation of high-dimensional variables or other features, such as change points, that are difficult to do in a standard regression context. Furthermore, these methods have been well studied both from the theoretical perspective -correct estimation is guaranteed under a number of reasonable assumptions -and the optimization perspective -efficient algorithms exist to solve them. Not surprisingly, these approaches have gained significant popularity in both the statistics and machine learning communities. We will describe several relevant examples below.
In the general context, we consider finding the minimum of where L(u, y) is a smooth convex Lipshitz loss function, Φ(u) is a regularization term, and λ ≥ 0 controls the amount of the regularization. In statistical models, L is typically given by the loglikelihood of the model. The regularization term Φ(u) is required to be convex but is not required to be smooth; in practice it is is often a regularization penalty that enforces some sort of sparseness or structure. A common regularization term is simply the L 1 norm on u, i.e.
which gives the standard LASSO problem proposed in [Tibshirani, 1996] and described in section 2. A number of algorithms have been proposed for the optimization of this particular problem Friedman et al. [2010], Efron et al. [2004]. Most recently, a method using proximal operators has been proposed called the Fast Iterative Soft-Threshold Algorithm [Beck and Teboulle, 2009]; this achieves theoretically optimal performance -identical to the O 1 k 2 lower bound of the convergence rate by iteration of optimization of general smooth, convex functions [Nesterov, 2005]. The method is based on the proximal operators we discuss next, and is described in detail, along with theoretical guarantees, in Beck and Teboulle [2009]. This is the method that immediately fits with our theoretical results.

Proximal Operators
Proximal gradient methods Nesterov [2007], Beck and Teboulle [2009] are a very general approach to optimizing (19). At a high level, they can be understood as a natural extension of gradient and sub-gradient based methods designed to deal efficiently with the non-smooth component of the optimization. The motivation behind proximal methods is the observation that the objective γ(u) is the composition of two convex functions, with the non-smooth part isolated in one of the terms. In this case, the proximal problem takes the form where we abbreviate f (u) = L(u, y). L is the Lipshitz constant of f , which we assume is differentiable. We rewrite this as The minimizing solution of (22) forms the update rule. Intuitively, at each iteration, proximal methods linearize the objective function around the current estimate of the solution,û, then update this with the solution of the proximal problem. As long as (21) can be optimized efficiently, the proximal operator methods give excellent performance. While this technique has led to a number of algorithmic improvements for simpler regularizers, it has also opened the door to the practical use of much more complicated regularization structures. In particular, there has been significant interest in more advanced regularization structures, particularly ones involving dependencies between the variables. These types of regularization structures are the primary focus of our work.

Related Work
Finally, of most relevance to our work here, is a series of recent papers by Francis Bach and others, namely Bach, Jenatton, Mairal, and Obozinski [2012], Mairal, Jenatton, Obozinski, and Bach [2010], , that independently discover some of our results. In these papers, combined with the more general work of , Bach [2010b] and Bach [2010a], ?, 2011], several of the results we present here were independently proposed, although from a much different starting point. In particular, he proposes using parametric network flows to solve the proximal operator for Φ(u) = i,j w ij |u i − u j | and that this corresponds to calculating the minimum norm vector of the associated submodular function. These results are discussed in more detail in section 4. While our results were discovered independently from these, we extend them in a number of important ways. First, we propose a new algorithm to exactly solve the linear parametric flow problem; this is a core algorithm in these problems. Second, through the use of the weighting scheme proposed in section 5, we develop a way to include general convex linear functions directly in the optimization procedure. This expands the types of regularization terms that can be handled as part of the proximal operator scheme.

Combinatorial Optimization and Network Flows
The final piece of background work we wish to present forms the last pillar on which our theory is built. Core to the theory and the algorithms is the simple problem of finding the minimum cost partitioning of a set of nodes V in which relationships between these nodes are defined by a graph connecting them. Our results are built on a fundamental connection between a simple combinatorial problem -finding a minimum cost cut on a graph -and the convex optimization problems described in section 3.2.
A fundamental problem in combinatorial optimization is finding the minimum cost cut on a directed graph. This problem is defined by (V, E, c, s, t), where V is a set of nodes and E is a set of directed edges connecting two nodes in V. We consistently use n to denote the number of nodes and, without loss of generality, assume V = {1, 2, ..., n}, i.e. the nodes are refered to using the first n counting numbers. Similarly, the edges are denoted using pairs of vertices; Here, (i, j) denotes an edge in the graph going from i to j. For an undirected graph, we simply assume that (i, j) ∈ E if and only if (j, i) ∈ E, and that c ij = c ji for all i, j ∈ V. c associates a cost with each of the edges in E. I.e. maps from a pair of edges to a non-negative cost, i.e. c ij ∈ R + for (i, j) ∈ E.
The letters s and t here denote specific nodes not in the node set V. For reasons that will become clear in a moment, s is called the source node and t is called the sink node. We treat s and t as valid nodes in the graph and define the cost associated with an edge from s to i as c si ; analogously, c it denotes the cost of an edge from i to the sink t. These nodes are treated specially in the optimization, however; the minimum cut problem is the problem of finding a cut in the graph separating s and t such that the total cost of the edges cut is minimal. Formally, we wish to find a set S * ⊆ V satisfying where Argmin with a capital A returns the set of minimizers, as there may be multiple partitions S * achieving the minimum cost. The above lays out the basic definitions needed for our work with graph structures. This problem is noteworthy, however, as it can be solved easily by finding the maximal flow on the graph -one of the fundamental dualities in combinatorial optimization is the fact that the cut edges defining the minimum cost partition are the saturated edges in the maximum flow from s to t in the equivalent network flow problem. It is also one of the simplest practical examples of a submodular function, another critical component of our theory. We now describe these concepts.

Network Flows
The network flow problem is the problem of pushing as much "flow" as possible through the graph from s to t, where the capacity of each edge is given by the cost mapping c above. The significance of this problem is the fact that it map directly to finding the minimum cost partition in a graph [Dantzig and Fulkerson, 1955, Cormen et al., 2001, Kolmogorov and Zabih, 2004. The saturating edges of a maximizing network flow -those edges limiting any more flow from being pushed through the graph -defines the optimal cut in the minimum cut partitioning. This result is one of the most practical results of combinatorial optimization, as many problems map to the partitioning problem above, and one and simple algorithms exist for solving network flow problems.
In the network flow problem, we wish to construct a mapping z that represents "flow" from s to t. A mapping z is a valid flow if, for all nodes in V, the flow going into the each node is the same as the flow leaving that node. Specifically, where for convenience, we assume that z ii = 0 and z ij = c ij = 0 if (i, j) / ∈ E. Furthermore, we say that an edge (i, j) is saturated if z ij = c ij , i.e. the flow on that edge cannot be increased.
Since we frequently talk about flows in a more informal sense, we use the term "valid flow" to reference this formal definition.
It is easy to show that the amount of flow leaving s is the same as the amount of flow leaving t, this total flow is what we wish to maximize in the maximum flow problem, i.e. we wish to find a maximal flow z * such that The canonical min-cut-max-flow theorem [Dantzig andFulkerson, 1955, Cormen et al., 2001] states that the set of edges saturated by all possible maximal flows defines a minimum cut in the sense of (23) above. The immediate practical consequence of this duality is that we are able to find the minimum cut partitioning quickly as fast algorithms exist for finding a maximal flow z * . We outline one of these below, but first we generalize the idea of a valid flow in two practically relevant ways.

Preflows and Pseudoflows
Two extensions of the idea of a flow z on a graph relaxes the equality in the definition of a flow. Relaxing the inequality is key to several algorithms, and forms some interesting connections to the theory we outline. The equality constraint is effectively replaced with an excess(·) function that gives the excess flow at each node; if the total amount of flow into and out of a node i is equal, then excess(i) = 0. Formally, With this definition, we define the excess(·) function as It is easy to see that if z is a preflow, excess(i) ≥ 0 for all nodes i. Maintaining a valid preflow is one of the key invariants in the common push-relabel algorithm discussed below.
A Pseudoflow is the weakest definition of a flow; it relaxes (29) completely, allowing there to be both excesses and deficits on the nodes of the graph. However, it does add in an additional constraint, namely that all edges connected to the source and sink are saturated. Formally, Note that now the excess(i) function can be either positive or negative. A pseudoflow has some interesting geometric and algorithmic properties, It was described in Hochbaum [1998] and an efficient algorithm for solving network flows based on pseudoflows was presented in Hochbaum [2008]. We mention it here to preview our results in section 4: the pseudoflow matches up directly with the base polytope, one of the fundamental structures in submodular function optimization.

Network Flow Algorithms
Finding the minimum cut or the maximum flow on a network is one of the oldest combinatorial optimization problems, and, as such, it is also one of the most well studied. It would be impossible to detail the numerous algorithms to solve the network flow problem here - Schrijver [2003] lists over 25 survey papers on different network flow algorithms. Many of these are tailored for different types of graphs (sparse vs. dense) and other variations of the problem.
Perhaps the most well-known algorithm for Network flow analysis is the push-relabel method. Variants of it achieve the best known performance guarantees for a number of problems of interest. It is simple to implement; essentially, each node has a specific height associated with it that represents (informally) the number of edges in the shortest unsaturated path to the sink. Each node can have an excess amount of flow from the source. At each iteration, a node with excess either pushes flow along an unsaturated edge to a neighbor with a lesser height, or increments its height. Termination occurs when the source node is incremented to |V| + 1, where |V| is the number of nodes -at this point, there is no possible path to the sink and the maximum cut can be found. For more information on this, along with other information on this algorithm, see Cormen et al. [2001] or Schrijver [2003]. This algorithm runs at the core of our network flow routines. This algorithm is guaranteed to run in polynomial time; in the general case, it runs in O |V| 2 |E| time, where V is the set of vertices and E is the set of edges. However, it can be improved to O |V| |E| log( V 2 / |E|) using dynamic tree structures [Schrijver, 2003].

The Parametric Flow Problem
One variation of the regular network flow problem is the parametric flow problem [Gallo et al., 1989]. The general parametric flow problem is a simple modification of the maximum flow problem given above, except that now c si is replaced with a monotonic non-decreasing function c si (β), Gallo et al. [1989] showed that this problem could be solved for a fixed sequence β 1 < β 2 < · · · < β m in time proportional to the time of a single run of the network flow problem by exploiting a nestedness property of the solutions. Namely, as β increase, the minimum cut of the optimal solution moves closer on the graph to the source. More formally, if β 1 < β 2 , then where S * (β) is an optimal cut at β in the sense of (23). If c si (β) is a linear function with positive slope, we call this problem the the linear parametric flow problem. In sections 4 and 5, we develop an exact algorithm for this problem that does not require a predetermined sequence of β to solve. This method is one of the key routines in the algorithms for general optimization.

Connections to Statistical Problems
One of the recent applications of graph partitioning is in finding the minimum energy state of a binary Markov random field. In particular, if we can model the pairwise interactions over binary variables x 1 , ..., x n on the graph using unary and pairwise potential functions E 1 i (x i ) and E 2 ij (x i , x j ), then the minimum energy solution can be found using the maximum cut on a specially formulated graph, provided that E 2 ij (0, 0) + E 2 ij (1, 1) ≤ E 2 ij (0, 1) + E 2 ij (1, 0) for all i, j. This application of network flow solvers has had substantial impact in the computer vision community, where pairwise potential functions satisfying this condition are quite useful.
We discuss the details of this application in section 4, where we start with the theory behind these connections as first building block of our other results. Ultimately, we show that network flow algorithms can be used to solve not only this problem, but the much more general optimization problems described in section 3.2 as well. Through these connections, we hope to bring the power of network flow solvers into more common use in the statistical community.

Submodular Functions
The problem of finding a minimal partition in a graph is a special case of a much larger class of combinatorial optimization problems called submodular function minimization. Like the problem of finding a minimal-cost partitioning of a graph given in (23), this optimization problem involves finding the minimizing set over a submodular function f : 2 V → R, where 2 V denotes the collection of all subsets of a ground set V. Given this, f is submodular if, for all sets S, T ⊆ V, This condition is seen as the discrete analogue of convexity. It is sufficient to guarantee that a minimizing set S * can be found in polynomial time complexity. An alternative definition of submodularity involves the idea of diminishing returns [Fujishige, 2005]. Specifically f is submodular if and only if for all S ⊆ T ⊂ V, and for all i ∈ V\T , This property is be best illustrated with a simple example. Let A 1 , A 2 , ..., A n ⊆ A be n subsets of a larger set A, and define so f coverage (S) measures the coverage of the set ∪ i∈S A i . In this context, it is easy to see that f coverage (S) satisfies (37). f coverage (S) is an example of a monotone submodular function as adding new elements to S is only going to increase the value of f coverage (S). However, many practical examples do not fall into this category. In particular, the graph cut problem given in equation (23) above is non-monotone submodular: We examine this particular example in more detail in section 4, where we prove f gc is indeed submodular. Submodular function optimization has gained significant attention lately in the machine learning community, as many other practical problems involving the sets can be phrased as submodular optimization problems. It has been used for numerous applications in computer vision, language modeling Bilmes, 2010, 2012], clustering [Narasimhan, Jojic, andBilmes, 2005, Narasimhan andBilmes, 2007], computer vision Bilmes, 2011, 2010], and many other domains. This field is quite active, both in terms of theory and algorithms, and we contribute some novel results to both areas in section 5.
The primary focus of our work has been developing a further connection between these graph problems and submodular optimization theory. Our approach, however, is the reverse of much of the previous work. Ultimately, we attempted to map the network flow problems back to the submodular optimization problems. Surprisingly, this actually opened the door to several new theoretical results for continuous optimization, and, in particular, to efficient solutions of the optimization problems given in section 2.

Some Formalities
The theory surrounding submodular optimization, and combinatorial optimization in general, is quite deep. Many of the results underlying submodular optimization require a fairly substantial tour of the theory of the underlying structures; good coverage of these results is found in [Fujishige, 2005] and [Schrijver, 2003]. Most of these results are not immediately relevant to our work, so we leave them to the interested reader. However, several additional results are needed for some of the proofs we use later.
As with the discussion of problem, we assume that V = {1, 2, 3, ..., n}; our notation intentionally matches that of the minimum cut problem defined above and is consistent throughout our work. We refer to V as the ground set. Formally, f can be restricted to map from a collection of subsets of 2 V , which we refer to consistently as D, with D ⊆ 2 V . In the context of submodular functions, D must be closed under union and intersection and include V as an element [Fujishige, 2005]. In general, and except for some of our proofs in section 5, D can be thought of as 2 V .

Submodular Function Optimization
The first polynomial time algorithm for submodular function optimization was described in [Grötschel et al., 1993]; it used the ellipsoid method from linear programming [Chvtal, 1983]. While sufficient to prove that the algorithm can be solved in polynomial time, it was impractical to use on any real problems. The first strongly polynomial time algorithms -polynomial in a sense that does not depend on the values in the function -were proposed independently in Schrijver [2000] and Iwata, Fleischer, and Fujishige [2001]. The algorithm proposed by Schrijver runs in O n 8 + γn 7 , where γ refers to the complexity of evaluating the function. The latter algorithm is O γn 7 log n , which may be better or worse depending on γ. The weakly polynomial version of this algorithm runs in O γn 5 log M , where M is the difference between maximum and minimum function values. Research in this area, however, is ongoing -a strongly polynomial algorithm that runs in O n 6 + n 5 γ has been proposed by Orlin [2009].
In practice, the minimum norm algorithm -also called the Fujishige-Wolfe Algorithm -is generally much faster, although it does not have a theoretical upper bound on the running time [Fujishige, 2005]. We discuss this algorithm in detail, as it forms the basis of our work. However, there are clear cases that occur in practice where this algorithm does not seem to improve upon the more complicated deterministic ones -in Jegelka, Lin, and Bilmes [2011], a running time of O n 7 was reported. Thus the quest for practical algorithms for this problem continues; it is an active area of research.
Additionally, several other methods for practical optimization have been proposed for general submodular optimization or for special cases that are common in practice. Stobbe and Krause [2010] proposed an efficient method for submodular functions that can be represented as a decomposable sum of smaller submodular functions given by g i (|U i ∩ S|), where g i is convex. In this particular case, the function can be mapped to a form that permits the use of nice numerical optimization techniques; however, many submodular functions cannot be minimized using this technique, and it can also be slow . Along with analyzing the deficiencies of existing methods, Jegelka, Lin, and Bilmes [2011] propose a powerful approach that relies on approximating the submodular functions with a sequence of graphs that permit efficient optimization. This method is quite efficient on a number of practical problems, likely because the graph naturally approximates the underlying structure present in many real-world problems.
Most recently, in Iyer, Jegelka, and Bilmes [2013], another practical optimization method is proposed; at its core is a framework for both submodular minimization and maximization based on a notion of discrete sub-gradients and super-gradients of the function. While not having a theoretical upper bound itself, it performs efficiently in practice. This method is also noteworthy in that it can constrain the solution space in which other exact solvers operate, providing substantial speedups.

Geometrical Structures and the Minimum Norm Algorithm
A number of geometrical structures underpin the theory of submodular optimization. In particular, an associated polymatroid is defined as the set of points in |V|-dimensional Euclidean space with sums of sets of the dimensions constrained by the function value of the associated set.
For notational convenience, for a vector x ∈ R n and set S ⊆ V, define In this way, x(S) forms a type of unnormalized set measure. Now the polymatroid associated with f is defined as P (f ) is fundamental to the theory behind submodular function optimization. In our theory, we work primarily with the base of the polymatroid P (f ), denoted by B(f ). This is the extreme (|V| − 1)-dimensional face of P (f ); it is defined as In the case where the domain D = 2 V , B(f ) is a linear, convex compact set. Thus it is often refered to as the Base Polytope of f as it is compact.
The base B(f ) is particularly important for our work, as one of the key results from submodular function optimization is the minimum norm algorithm, which states effectively that sets S * minimizing f are given by the sign of the point in B(f ) closest to the origin. This surprising result, while simple to state, takes a fair amount of deep theory to prove for general f ; we state it here: and let S * 1 = {i : y * i < 0} and S * 2 = {i : y * i ≤ 0}. Then both S * 1 and S * 2 minimize f (S) over all subsets S ⊆ V. Furthermore, for all S † ⊆ V such that f (S † ) = min S∈V f (S), S * 1 ⊆ S † ⊆ S * 2 . A central aspect of our theory relies on this result. In general, it as one way of working the geometric structure of the problem. It turns out that in the case of graph partitioning, B(f ) takes on a particularly nice form. This allows us to very quickly solve the minimum norm problem here. We show the resulting minimum norm vector also gives us the full solution path over weighting by the cardinality of the minimizing set. While this basic result was independently discovered in , our approach opens several doors for theoretical and algorithmic improvements, which we outline in the subsequent sections.

The Combinatorial Structure of Dependent Problems
The foundational aspect of our work is the result established in this section, namely an exploration of network flow minimization problems in terms of their geometric representation on the base polytope of a corresponding submodular problem. This representation is not new; it is explored in some depth as an illustrative example in [Fujishige, 2005] and connected to the minimum norm problems in . Our contribution, however minor, is based on a simple transformation of the original problem that yields a particularly intuitive geometric form. The fruit of this transformation, however, is a collection of novel results of theoretic and algorithmic interest; in particular, we are able to exactly find the optimal u * over the problem Recall that this problem was discussed in depth in section 3.2; here, a ∈ R n is given, and λ ∈ R + and the weights w ij ∈ R + control the regularization.
In this section, we lay the theoretical foundation of this work, denoting connections to other related or previously known results. The contribution at the end is a strongly polynomial algorithm for the solution of a particular parametric flow problem; this algorithm follows naturally from this representation. In the next section, we extend the theory to naturally allow a convex piecewise-linear function ξ i (u i ) to be included in (44) to match ( B ).
Before presenting the results for the real-valued optimization of (44), we must first present a number of theoretical results using the optimization over sets V = {1, 2, ..., n}. We begin with the simplest version of this optimization -finding the minimizing partition of a graph -and extend this result to general submodular functions later.

Basic Equivalences
We proceed by defining and establishing equivalences between three basic forms of the minimum cut problem on a graph. Recall from section 3.3 that the minimum cut problem is the problem of finding a set S * ⊆ V satisfying As we mentioned earlier, several other important problems can be reduced to this form; in particular, finding the lowest energy state of a binary Markov random field, when the pairwise potential functions satisfy the submodularity condition, is equivalent to this problem. Our task now is to make this explicit. We here show equivalences between three versions of the problem. The first is P E , which gives the standard energy minimization formulation, i.e. finding the MAP estimate of which is equivalent to finding P Q gives the formulation as a quadratic binary minimization problem; here, the problem is to find x ∈ {0, 1} n that minimizes x T Qx for an n × n matrix Q. This version forms a convenient form that simplifies much of the notation in the later proofs. Finally, we show it is equivalent to the classic network flow formulation, denoted P N . This states the original energy minimization problem as the minimum st-cut on a specially formulated graph structure. The equivalence of these representations is well known [Kolmogorov and Zabih, 2004] and widely used, particularly in computer vision applications. We here present a different parametrization of the problem which anticipates the rest of our results.
The other unique aspect of our problem formulation is the use of a size biasing term; specifically, we add a term β |S| to the optimization problem, where β ∈ R can be positive or negative. This term acts similarly to a regularization term in how it influences the optimization, but to think of it this way would lead to confusion as the true purpose of this formulation is revealed later in this section -ultimately, we show equivalence between the values of β at which the set membership of a node flips and the optimal values of the continuous optimization problem of ( B ).
Theorem 4.1. Let G = V, E) be an undirected graph, where V is the set of vertices (assume V = {1, ..., n}) and E ⊆ {(i, j) : i, j ∈ V} is the set of edges. Without loss of generality, assume that i < j ∀ (i, j) ∈ E. Define S * E (β), S * Q (β), and S * N (β) as the sets of optimizing solutions to the following three problems, respectively: Energy Minimization: P E (β). Given an energy function and let S * Quadratic Binary Formulation: P Q (β). Given E 1 and E 2 as in P E (β), define the n × n matrix Q = [q ij ] as: Suppose q ij ≤ 0 for i = j, and let and let S * Minimum Cut Formulation: P N (β). Let G = (V , E ) be an augmented undirected graph with V = V ∪ {s, t}, where s and t represent source and sink vertices, respectively, and Define capacities c ij , (i, j) ∈ E on the edges as: where Then the set of minimum cut solutions S * N (β) is given by Then P E (β), P Q (β), and P N (β) are equivalent in the sense that any minimizer of one problem is also a minimizer of the others; specifically, Proof. A reformulation of known results [Kolmogorov and Zabih, 2004], but given in section A for convenience.
The primary consequence of this theorem is that solving the energy minimization problem can be done efficiently and exactly due to several types of excellent network flow solvers that make solving problems with millions of nodes routine [Boykov and Kolmogorov, 2004, Cormen et al., 2001, Schrijver, 2003]. Because of this, numerous applications for graphcuts have emerged in recent years for computer vision and machine learning. Our purpose, in part, is to expand the types of problems that can be handled with network flow solvers, and problems of interest in statistics in particular.
In our work, we alternate frequently between the above representations. For construction, the energy minimization problem is nicely behaved. In the theory, we typically find it easiest to work with the quadratic binary problem formulation due to the algebraic simplicity of working in that form. Again, however, each of these is equivalent; when it does not matter which form we use, we refer to the problem an solution set as S * (β) and P(β) respectively.

Connections to Arbitrary Network Flow Problems
While theorem 4.1 lays out the equivalence between P E (β) and P Q (β) and a specific form of network flow problem P N (β), for completeness we show that any network flow problem can be translated into the form P N (β) and thus P E (β) and P Q (β). The two distinguishing aspects of P N (β) are the facts that each node is connected to either the source or the sink, and that all the edges are symmetric, i.e. c ij = c ji for all i = j. It is thus sufficient to show that an arbitrary flow problem can be translated to this form. For conciseness, assume that β = 0; the results can be adapted for other β easily.
Theorem 4.2. Any minimum cut problem on an arbitrary, possibly directed graph can be formulated as a quadratic binary problem of the form P Q (β = 0) as follows: 1. For all edges (i, j) such that c ij > c ji , add a path s → j → i → t with capacity c ij − c ji . That edge is now undirected in the sense that both directions have the same capacity, and the edges in the minimum cut are unchanged, as these paths will simply be eliminated by flow along that path. 2. Set q ij = c ij for i < j and q ij = 0 for i > j.
3. Given q ij , set q ii = (c si − c it ) − 2 i :i <i q i ,i + j : i<j q ij . Using these steps, any minimum cut problem can be translated to P Q in the sense that the set of minimizing solutions is identical.
Proof. Network flows are additive in the sense that increasing or decreasing the capacity of each edge in any path from s to t by a constant amount does not change the set of minimizing solutions of the resulting problem, even if new edges are added [Cormen et al., 2001, Schrijver, 2003]. Thus step (1) is valid. The rest follows from simple algebra.

Network Flows and Submodular Optimization
Recall from section 3.4 that the minimum cut problem is a subclass of the more general class of submodular optimization problems. In the context of P Q (β), it is easy to state a direct proof of this fact, additionally showing that here the submodularity of the pairwise terms is also necessary for general submodularity.
Theorem 4.3. The minimization problem P Q (β) can be expressed as minimization of a function Then f β is submodular if and only if q ij ≤ 0 ∀ i, j ∈ S, i < j.
Proof. One immediate proof of the if part follows from the fact that P Q (β) can be expressed as the sum of submodular pairwise potential terms, and the sum of pairwise submodular functions is also submodular [Fujishige, 2005]. The direct proof, including both directions, is a simple reformulation of known results see [Kolmogorov and Zabih, 2004], but given in section A on page 38 for convenience.

Geometric Structures
We are now ready to present the theory that explicitly describes the geometry of P Q (β), which extends to both P E (β) and P N (β), in the context of submodular function optimization. This theory is not new; several abstract aspects of it have been thoroughly explored by Schrijver [2003] and Fujishige [2005]. In our case, however, the exact form of the problem presented in theorem 4.1 was carefully chosen to yield nice properties when this connection is made explicit. From these, a number of desirable properties follow immediately. The rest of this section is arranged as follows. First, we show that the base polytope B(f β ) is a reduction of all pseudoflows (see definition 3.3) on the form of the cut problem P N (β) from theorem 4.1. B(f β ), described in section 3.4, has special characteristics for our purposes, as the minimum norm algorithm described in section 3.4 provides a convenient theoretical tool to examine the structure of P Q (β). Our key result is to show that the minimum norm vector -the L 2 -projection of the origin onto B(f β ) -depends on β only through a simple, constant shift. Thus this vector allows us to immediately compute the minimum cut directly for any β, and we are thus able to compute the minimum cut solution as well for any β as well.
Theorem 4.4 (Structure of B(f β )). Let f β be defined in theorem 4.3 (58), and let Let r i (α), i = 1, ..., n, be defined as follows: and denote r(α) = (r 1 (α), ..., r n (α)) ∈ R n . Then the base of the polymatroid associated with f β , B(f β ), is given by and the full polymatroid polytope is given by Proof. Proceeds with straightforward albeit tedious algebra. Proved in section A on page 37.
In light of this, the minimum norm problem on the graph structure is as follows. The immediate corollary to the min-norm theorem is that α * (β) yields the optimum cut S * (β) of the corresponding network flow problem, P N (β): Theorem 4.5 (Minimum Norm Formulation of P N (β)). Let A and r(α) be given by equations (59) and (60), respectively. Then the min-norm problem P N (β) associated with P Q (β) is defined by Then any optimal cut S * (β) solving P G (β), as given in theorem 4.1, satisfies: Proof. Follows immediately from theorems 4.1 and 4.3 characterizing the cut problem as a submodular function optimization problem, theorem 4.4 describing the structure of this problem, and theorem 3.4 to characterize the solution.
The optimal α * in the above formulation has some surprising consequences that motivate the rest of our results. In particular, we show that the optimal α * in equation (63) is independent of β; this is the key observation that allows us to find the entire regularization path over β. Formally, this result is given in theorem 4.9. However, we need other results first.

Connection to Flows
The representation in terms of α is significant partly as the values of α effectively form a pseudoflow in the sense of Hochbaum [2008] (see section 3.3.2). Recall that a pseudoflow extends the concept of a flow by allowing all the nodes to have both excesses and deficits. In addition, a pseudoflow assumes that all edges from the source to nodes in the graph are saturated, possibly creating excesses at these nodes, and all edges connected to the sink are similarly saturated, possibly creating deficits at these nodes.
Theorem 4.6. Consider the problem P N (β). For any α ij ∈ R, with (i, j) ∈ E, let α ij represent the flow on each edge c ij , with α ij > 0 indicating flow from i to j, and α ij < 0 indicating flow from j to i. Then α ∈ A defines a pseudoflow on the graph structure indexed by A. Furthermore, r i (α) − β = excess(i) is the (possibly negative) excess at node i in the sense of (30).
Proof. The pseudoflow condition that all edges from the source node and to the sink node are saturated is immediately implied by the fact that r i (0) − βw i = (c si − c it ) − β. The flow conditions, then follow from the edge capacity being c ij = c ji = −q ij and −|q ij | ≤ α ij ≤ |q ij |. Proof. Follows immediately from theorem 4.6.

Structure of the Complete Solution
The theorem above has a number of important consequences detailed in the next few sections. The most immediate consequence comes when we consider the structure of the optimal solution of the minimum norm algorithm specialized to the network flow problem P N (β); this effectively allows us to derive a way of solving P N (β) for all β.
Lemma 4.8 (Optimal solutions to P N (β)). Then α * is an optimal solution to P N (β) if and only if for all i, j, i < j, the following condition holds: In particular, the optimum value of α * in this case is independent of β.
Proof. First, A is convex as each dimension α ij is bounded independently. Thus the objective of P N (β) is minimizing a convex function over a convex domain. Therefore, it suffices to prove that equation (65) can be satisfied if and only if α * is a local minimum of r(α) − β 2 2 . As r(α) − β 2 2 is differentiable w.r.t. α, this is equivalent to showing that either the gradient is 0 or α * is on the boundary of A and all coordinate-wise derivatives point outside the domain A. First, define g ij (α) as the gradient of r(α) − β 2 2 w.r.t. α ij : Thus, and the following condition holds for all α: Matching these conditions to those in (65) shows that α * as given defines a local, and thus global, optimum of P N (β). In particular, note that this criterion is independent of β, completing the proof.

Invariance to β
The invariance of the optimal α * to β allows us to characterize the solution space of optimal cuts as the level sets of r(α * ). The core result, as well as our algorithm, is based on this intuition. I. α * is the optimal solution to (63) if and only if for all β ∈ R, all optimal cuts S * β ∈ S * (β) for P(β) satisfy II. Furthermore, for all β, U 1 (β) is the unique smallest minimizer of P(β) and U 2 (β) is the unique largest minimizer.
Proof. Part (I) follows as a direct consequence of theorem 4.5 and the invariance of the minimum norm problem to changes in β as given in lemma 4.8. Part (II) is then an immediate consequence of the minimum norm algorithm.
This theorem is intuitively important as the key values that permit a connection to the continuous problems are values of β at which the membership of the different nodes change. This theorem tells us that these points are given by the values of the minimum norm vector, here given as r(α * ).

Monotonicity
One immediate corollary of theorem 4.9 is a monotonicity property on the optimal sets, used in the continuous optimization theory we present later: Corollary 4.10. Let β 1 < β 2 . Then for all S * 1 ∈ S * (β 1 ) and S * 2 ∈ S * (β 2 ), Proof. Follows immediately from the equivalence of the optimizing sets to the level sets of the minimum norm vector given in theorem 4.9.

Beyond Network Flows
Theorem 4.9 above was discovered independently for the full case of general submodular functions by Nagano et al. [2011]. There, the authors similarly showed that the level sets of the minimum norm algorithm give the solutions for the f (S) − β |S| problem. While the approach those authors take is different and more involved, we give a shorter, alternative proof. We use a simple argument following from the fact that B(f ) constrains the minimum norm vector y to a constant total sum. The result is that the constant offset in the minimum norm objective drops out of the optimization. More formally: Theorem 4.11 (Invariance of General Submodular Functions to β).
I. Let f be a general submodular function. Then y * is the optimal solution to the minimum norm problem if and only if ∀ β ∈ R, the set of optimizing solutions II. Furthermore, for all β, U 1 (β) is the unique minimal solution to (63) and U 2 (β) is the unique maximal solution in S * (β).
Proof. Consider the submodular function f β (S) = f (S) − β |S|. It is easy to show that Denote by y * (β) the minimum norm solution for f β . Then the submodular problem for y * (β) is given by where steps (78)-(79) follow by definition of B(f ), causing the terms dependent on β to drop out by way as constants under the optimization. From this, we have that the optimal minimum norm vector for f β is just the minimum norm vector for f shifted by β, immediately implying part (I). Similarly, part (II) follows immediately from the minimum norm theorem.
This result is used in several other sections as well, and has a number of practical implications for size-constrained optimizations and related problems. For a full treatment of related implications, see [Nagano et al., 2011].

Exact Algorithm for Constant Parametric Flows
Using the above theory, we now wish to present a viable algorithm to calculate the reduction, and hence all cuts, for each node. The idea is simple and follows immediately from the similarity of the minimum cut problem P N (β) to the structure of B(f β ) as detailed in theorem 4.6. As each level set of r(α * ) defines an optimum cut in the graph, we can adjust all of the unary potentials by r µ = mean i∈S , chosen to bisect the reduction values, and solve the resulting cut problem. By the max-flow min-cut theorem, all edges crossing the cut are saturated. Specifically, As these α ij are optimal in the final solution α * , they can be fixed by permanently adding their values to the corresponding r i (α) and removing them from consideration in the optimization. This then bisects the nodes, allowing us to treat these two subsets separately when solving for the rest Algorithm 1: AlphaReduction Input: Submodular Q. Output: α * , the minimizer in A of r(α) 2 2 . // Begin by calling BisectReductions below on the full set V to get α. return BisectReductions (T = V, α = 0, Q) // α [T ] denotes α restricted to edges with both nodes in T .
if r µ = r min then return α [T ] // Done; We are on a single level set.
// Fix the flow on edges in the cut by adjusting the source/sink capacities of each node, then removing those edges.
of the bisections. The validity of this bisection can also be seen by the optimality of the minimum norm solution as described by theorem 4.9.
Algorithm 1 can be summarized as follows. We first start by considering the entire set of nodes, setting the working set S = V. At each step, we recursively partition the working set S using a minimum cut as follows: 1. If r(α) is constant in S, then return. We're done. 2. Otherwise, set up a network flow problem to bisect the nodes and find a minimum cut. Once a minimum cut is found, set all the edges in the cut to their saturated values. 3. Repeat on the two resulting subsets of nodes. Pseudocode for this algorithm is presented in Algorithm 1.
Theorem 4.12 (Correctness of Algorithm 1.). After the termination of Algorithm 1, all values of α are set such that r(α) 2 is minimized over α ∈ A.
Proof. Let i, j, i = j, be any pair of nodes such that q ij = 0, and let α † be the solution returned by Algorithm 1.
Next, suppose that r i (α) = r j (α), and first suppose that r i (α) > r j (α). Then, by the termination condition of the recursion in BisectReductions, nodes i and j must have been separated by a valid minimum cut for some r µ . However, as all edges crossing a minimum cut are saturated by the max-flow-min-cut theorem, α ij = |q ij |. Thus condition (65) in Lemma 4.8 is satisfied. Similarly, if r i (α) < r j (α), then α ij = −|q ij |, indicating a flow of |q ij | from j to i; again, this satisfies (65).
As the above holds for any pairs of nodes i, j, the optimality criteria of Lemma 4.8 is satisfied globally, proving the correctness of the algorithm.
This algorithm is quite efficient in practice, and it forms an core routine of the total variation minimization algorithm given in the second part of our work, where we present full experiments and some comparisons with existing approaches.

Extensions to Real-valued Variables
One of the intriguing consequences of the above theory, and one that opens new doors to efficiently optimizing several other classes of functions, comes as the result of being able to map other correlated data to this framework. In general, interactions between terms can be very difficult to work with in practice. However, the above theory allows us to exactly find the optimizer of a large class of general functions. These functions may not necessarily be smooth.
Our approach connects closely to several recent results discovered independently by Mairal  and Bach [Bach, 2010a], which connect some of these problems to an older result by Hochbaum [Hochbaum and Hong, 1995]. The last of these papers effectively establishes an equivalence between a class of quadratic objective functions and some types of network flow algorithms, although the equivalence to parametric flows isn't really explored. However, this result was used by Bach in Bach [2010a] to note that the minimum norm problem of the network flow problem can be solved using classical methods for solving parametric flow problems [Gallo et al., 1989], and the implications of this for structured sparse recovery are explored in . The end result, discovered independently, parallels the theorem we present below, albeit with a different algorithm.
In contrast, while less general, the algorithm we presented in 1 gives the exact change points immediately, and the theoretical framework presented surrounding this problem is more thoroughly explored here. However, the primary practical improvement we provide comes in the next chapter when we incorporate the use of piecewise-linear convex penalty terms as well.
Theorem 4.13. Suppose γ : R n → R + can be expressed as where u ∈ R n is the variable we wish to optimize over and a ∈ R n , λ > 0, and w ij ≥ 0 are given. Without loss of generality, assume that i < j. Then the minimizer can be found exactly using Algorithm 1 with Then u * = r(α * ).
Proof. See appendix A, page 41.

Unary Regularizers and Non-Uniform Size Measures
In many contexts, it is helpful to use regularization terms or weighting terms on the solution to control the final behavior of the result. In the previous section, we discussed the simplest case in the submodular function context, namely weighting the problem against the cardinality of the solution set. We proved this for the case of graph-based submodular problems and extended that argument to the full submodular function context. In this section, we extend this result to the case where the size biasing term β |S| term is replaced with a weighted size biasing term βw(S). Here, is the weighted size biasing term on the function. We then wish to find for all β ∈ R and all positive weight measures w.
The results of the rest of these papers are entirely novel. We prove that the optimal solution to an alternate construction of the minimum norm theorem yields the entire solution path for all β. Like the last section, we find a vector z * on the base B(f ) such that every solution is given by a zero-crossing of [z * − βw]. We then show that this allows us to include more detailed structures in the continuous optimization problem; in particular, we are able to incorporate the piecewise-linear convex penalty term ξ i (u i ) in ( B ) directly into our optimization.
Our result is based around a very simple technique to augment the original graph such that auxiliary variables "attract" parts of the regularization influence of β and transfer it to associated variables in the original problem. We then show that it is possible to translate this result into arbitrary weights by taking several well-controlled limits. The end result is an algorithm for solving (89) for arbitrary positive weights.
As this result is novel and holds for general submodular functions, we prove it for the general case first, then extend it to the special case of the linear parametric flow problem. Analogously to algorithm 1, the algorithm we develop here solves this problem exactly. The fact that the solution is exact also allows us to solve ( B ).

Encoding Weights by Augmentation
The primary tool used for introducing weights into the optimization of the level sets of the function is to augment the original problem with additional variables. When β = 0, these variables do not contribute to the solution values of the base set of nodes. In the network flow interpretation, they have no connection to the source or sink -but they are subject to the influence by β in the resulting solutions.
With the proper construction, it is possible to guarantee that these augmented nodes always have the same reduction value as the nodes they are augmenting; this allows us to construct a graph such that these values then translate back into weights on the β terms. The primary tool used is the following lemma, which forms the basis of the rest of our results.
where Argmin returns the set of minimizing sets, f w (T ) is submodular and given by and K i is a block of indices of length w i − 1, given by II. Define Then D w is a distributed lattice and (90) can be replaced by III. Furthermore, Proof. See appendix A, page 42.
The above lemma is noteworthy as it provides a way to theoretically augment the original problem in a way that alters the original problem such that the relative influence of the β scaling can be altered. In particular, in the augmented problem, the size of the evaluation set |T | includes these augmented nodes -since they are included deterministically based on the values in the unaugmented set V, the unaugmented node is effectively counted multiple times. This allows us to weight the nodes separately.
In our context, when dealing with graph structures, this corresponds to adding a collection of single nodes with no connections other than an effectively infinite capacity edge connecting each to one of the base nodes. As this edge ties the nodes together in any cut solution, the influence of the global weighting parameter β on this auxiliary node is simply transferred to the attached node. The next two theorems extend this result to the minimum norm vector y * , and an approximation lemma extends this to general weight vectors.

Optimization Structure and The Weighted Minimum Norm Problem
The central result of this section is a weighted version of the minimum norm problem. Under this construction, the solution to the original minimum norm problem is the same as problem a β |S| weighting term, but with additional nodes. However, the level sets of the resulting vector yield the optimal minimizing sets f (S) − βw(S) for all values of the parameter β.
Definition 6.1 (Weighted Minimum Norm Problem). For a submodular function f defined on D ⊆ 2 V , and positive weights w ∈ R n , w > 0, the weighted minimum norm problem is given by and we call the solution vector z * the Weighted Minimum Norm Vector. Furthermore, if the weights w are restricted to be positive integers, then this is called the Integer Weighted Minimum Norm Problem.
We use the solution to this problem in an analogous way to the use of the minimum norm vector of section 4. The theorems in this section show that z * gives the entire solution path to f (S)−βw(S) over β.
Theorem 6.2 (Structure). Let f be a submodular function on D, and let w, V w , f w ,K i , and D w be as defined in lemma 5.1 (in particular, the elements of w are integers). Let The polymatroid associated with f w is given by and the associated base polymatroid is given by II. Let y * = argmin y∈B(fw) Then III. Furthermore, let z * be the solution to the integer weighted minimum norm problem, i.e. Then Furthermore a vector z * is the optimal solution to (101) if and only y * , as given by (102), is the optimal solution to (99).
Proof. See appendix A, page 44.
The important concept behind this theorem and its corollaries is that it demonstrates a direct connection between the minimum norm problem on the augmented problem f w and the original problem f . This connection allows us to build the theory of the weighted problem directly upon the original theory, essentially using those results.

On the Use of General Positive Weights
The goal of this section is to extend the above results on integer w to all positive real numbers. This allows us to do a number of interesting things, particularly in the case of network flow algorithms. It also extends the state of the known theory on general submodular function minimization outside of the cases we are interested in. We here state the form of theorem 6.2 for general w, then discuss some of the implications for the case of network flows and the continuous optimization problems introduced earlier. In particular, this theorem allows us a way to include the piecewise linear ξ i (u i ) term in ( B ).
Theorem 6.3. Let f be a submodular function defined on D ⊂ 2 V , and let w ∈ R |V| be strictly positive, finite weights. Then I. Let z * be the optimal solution to the weighted minimum norm problem, i.e.
and, for all β ∈ R, let and let Then z * is the optimal solution to (103) if and only if, for all β ∈ R and all S * ∈ S * (β, w), II. Furthermore, for all β, U 1 (β) is the unique minimal solution to (106) and U 2 (β) is the unique maximal solution in S * (β, w).
Proof. This result is somewhat involved and quite technical. We present a proof of it, along with supporting lemmas, in section A.1 on page 45.
The above problem allows us to generalize the previous results of size-penalized submodular optimization to general weighted penalties. This result may have several significant practical implications; several of these we explore later in the context of the network flow analysis results.
An interesting corollary to the above theorem is that the original formulation of the minimum norm problem is still valid when the norm being optimized over is reweighted. It may be that this would open up an way to remove some of the numerical difficulties often encountered with the minimum norm problem . More formally, Corollary 6.4 (Validity of Weighted Minimum Norm.). Let z * be the optimal value of the weighted minimum norm problem, with w ∈ R |V| , w > 0. Let U 1 = {i : z * i < 0} and U 2 = {i : z * i ≤ 0}, Then both U 1 and U 2 minimize f (U ) over all subsets U ∈ D. Furthermore, for all U † ∈ D such that f (U † ) = min U ∈V , U 1 ⊆ U † ⊆ U 2 . In other words, the weighted minimum norm vector z * may be substituted for the original minimum norm vector.

Handling the Case of w i = 0
One of the challenging aspects here is that we might be interested in the case of w i = 0. In theory, this can be easily handled by simply allowing w i to be so small that its effect on the problem is negligibly different from w i = 0; in other words, we can see it as the limit w i 0. In practice, this leads to numerical issues. Thus we propose here a numerically stable method to work with w i = 0 by investigating the limiting behavior.
Theorem 6.5. Let w ≥ 0 and define Q = {i ∈ V : w i = 0}. Then let and let where z[Q] denotes the vector of elements of z in Q. Then, for we have that z * ε → z * as ε 0 (111) Proof. Consider the form of the optimization problem in (110). For ε sufficiently small, we have that As ε −1 ∞, the optimum value of z is constrained to be on the simplex in which i∈Q z 2 i is minimal. This value is given by C Q above, and this is the constraint that is enforced explicitly in (109). Since everything is continuous, it is valid to take the limit as ε 0. Thus the theorem is proved.
It is outside the current realm of our investigation how to implement this in the inner workings of the general minimum norm algorithm; however, we will revisit this issue later when proving the correctness of the network flow version of the weighted reduction algorithm.

Network Flow Solutions to and Unary Regularizers
We now turn our attention to the specific case of network flows. The linear parametric flow problem is similar to the flow problem described earlier, except that now we allow the capacity functionsanalogous to the unary energy terms -to be a non-decreasing linear function of the weighting term β. Previously, we treated this global weighting term as having equal influence on all nodes. This section considers the case where the influence of β has a different weight on each node. Specifically, we replace β with βw i in theorem 4.1. In this case, we are still able to compute the entire path directly. This result, while interesting in its own right, also sets the stage for our later results for total variation minimization.
To be specific, we extend the problems in theorem 4.1 as follows: Energy Minimization: P E (β, w). Let E 1 and E 2 be defined as in theorem 4.1, and let and let Quadratic Binary Formulation: P Q (β, w). Given Q defined as in theorem 4.1, let and let Minimum Cut Formulation: P C (β, w). Let the graph structure be defined as in theorem 4.1, and define capacities c ij , (i, j) ∈ E on the edges as: Then the minimum cut solution S * C (β) is given by It is a simple matter to show that the above are equivalent: Theorem 7.1. P * E (β, w), P * Q (β, w), and P * C (β, w) from the above description are equivalent in the sense that any minimizer of one problem is also a minimizer of the others, i.e.
Proof. Replace β with βw i or βw as appropriate in the proof of theorem 4.1.
However, it is a much more complicated endeavor to show that this formulation can be solved exactly in a similar manner to the previous result. In the end, we prove the following result: analogously to before, we find an optimal pseudoflow α * such that the zero crossings of r(α * ) − βw, with r(α) defined as in 4.4, give the level sets of the augmented problem. The purpose of the current section is to define these relationships explicitly.
The above result works as well for general network flow solutions as well. In this case, we have the following immediate corollary to theorem 6.3: Corollary 7.2. Let w be a collection of positive weights. Then α * is the optimal pseudoflow solution to the weighted minimum norm problem if and only if ∀ β ∈ R, the optimal cut S * (β, w) for P G satisfies Proof.
Replace B(f ) with its representation for network flow problems as given in theorem 4.4.
An alternate version of this, which is handy for the proofs, is the following: Corollary 7.3. α † is optimal for a problem if and only if, for all β ∈ R, the following condition holds for each ordered pair i, j ∈ V, i < j: Proof. Follows immediately from noting that criteria (122) in corollary 7.2 specifies that if there exists a β that separates two scaled reductions, then there is valid cut in the associated flow problem that separates these nodes. This, however, is equivalent to the condition on the α terms given by (123) or (124).

Anchoring Nodes with Differing or Infinite Weights
While the original interpretation of the weights is to vary the influence of each node in the regularization path, a natural extension in the network flow context is to use these weights to fix nodes at particular reduction values -informally, to make them less responsive to the influence of the external flow during the optimization. This is done by scaling both the base reduction value q ii and the β term by w i instead of just β. In this context, then, replace r(α) with r w (α), defined as and the node membership changes sign at As a result, it is possible to anchor nodes at a particular value by setting the corresponding weights to a high value and adjusting the corresponding source-to-sink weights to match. In particular, taking the limit as w i → ∞ in (127) causes the membership to be fixed at the internal reduction value. In this case, it effectively fixes these nodes at a pre-determined reduction value, even though the node still influences other nodes in the flow. In particular, this fact is used below to handle the case where we wish to look at the L 1 norm of a reduction value, as we explain later.

Algorithm
In the unweighted algorithm, one of the key routines was shifting the reduction level of all the nodes by the average on that region. This allowed us to effectively bisect the problem at that value of β; this in turn allowed us to set some of the edges as saturated, fixing them at their extreme values.
When the resulting solution on a region is that the flow solution sets all reductions to the average value -zero in the scaled version -that region is connected in the final solution, and we can set it to its common value. Now that we are putting weights on the nodes, the idea of an "average" is redefined. Now, it is the value of the β µ such that, if all the reductions values were in the same contiguous block, r i (α) − β µ w i would equal 0. This means that when β is increased or decreased, all of these nodes would flip at the same time.
Formally, this β value on a set T is given by: This behavior is consistent with the augmentation scheme discussed earlier -a node with weight 2 would count the same as 2 nodes in the "average." Now, when solving the resulting flow equation, it is important that the reduction values are scaled away from the mean such that shifts in the flow affect the reduction value on the correct scale. Since all we care about are the zero-crossings of the nodes, we then rescale all the starting reduction values away from the mean by w −1 i . In other words, we end up setting the unary potential value ρ i used in calculating the cut to Note that i∈T w i ρ i = 0, as we would expect, and that when used to split the regions, β µ will always separate the starting reductions ρ i into positive and negative components, ensuring a non-trivial cut solution when this is incorporated into the algorithm. As before, the resulting cut dictates the edges that will be saturated; this is how the solution in the rest of the problem is tracked. In the case of w i = 0, these nodes simply do not enter into the average; In rescaling them away from the mean, we set them to extreme positive or negative numerical values. This approach is consistent with theorem 6.5; it will force the squared reduction value in these nodes to be minimized.
Thus we can define a weighted version of the bisection routine, WeightedBisectionCut, to work with weighted nodes: Note that what is important once we are done with the algorithm is not the resulting values of r µ (T ), but that the edges spanning the cut are saturated. These are fixed per the same logic as in section 4.5.
Proof. This proof is nearly identical to the proof of algorithm 1; the difference is that we instead work with the condition of optimality given in corollary 7.3.
Note first that for w i = 0, the method used to set the values of ρ i on these nodes effectively sets them to their numerical extremes, ensuring that all operations seek to pull them as far as possible towards zero. Thus we are consistent with the behavior of theorem 6.5.
Let i, j, i < j, be any pair of nodes such that q ij = 0, and let α † be the solution returned by Algorithm 3.
First, suppose there existed a β such that Then trivially, the optimality criterion of Lemma 7.3 is satisfied.
Algorithm 2: WeightedBisectionCut Input: Submodular Q, weights w, subset of nodes T . Output: Bisecting cut S ⊆ T or ∅ if all nodes in same partition. ρ ← 0 if i∈T w i = 0 then // In this case, the optimal cut that divides the reductions into the positive or negative components is all that matters; w i = 0 means this result is true for all β.
// M here denotes the largest numerically stable number.
if ρ [T ] = 0 then return ∅ // These nodes are all at a common reduction level already.
, with capacities formed from Q by theorem 4.1. return S * T Next, suppose that there exists a β such that Then, according to the termination criterion of WeightedBisectionCut -in which it returns S * T = T or S * T = ∅ -the edge between nodes i and j will be saturated in such a way that the reduction value of r i (α † ) is maximized and the reduction value of r j (α † ) is minimized, i.e. α ij = −|q ij |. Thus condition (123) in corollary 7.3 is satisfied. Similarly, if the order in 131 is reversed, then condition (124) is satisfied.
As the above holds for any pairs of nodes i, j, the optimality criteria of corollary 7.3 are satisfied globally, proving the correctness of the algorithm.

General Optimization with Unary Regularizers
This section improves upon the result of theorem 4.13 given at the end of section 4: using the weighting method above, it is possible to anchor weights at specific values in the optimization. This indicates that the reduction levels of these nodes do not vary in the optimization even though they are treated the same way as the other nodes. This provides a simple way of incorporating an L 1 fit penalty into our optimization. Formally, the following corollary to theorem 4.13 lays this out: Theorem 8.1. Consider ( B ): suppose γ : R n → R + can be expressed as Algorithm 3: FindWeightedReductions Input: Submodular Q, vector of non-negative weights w. Output: α * satisfying condition 121.
// Begin by calling WeightedBisectReductions below on the full set V to get α. return WeightedBisectReductions (T = V, α = 0, Q) // Note: α [T ] denotes α restricted to edges with both nodes in T .
WeightedBisectReductions (T , α, Q) We are done on this set. return α [T ] // Fix the flow on edges in the cut by adjusting the source/sink capacities of each node, then removing those edges.
with a, b ∈ R n , λ ≥ 0, and w ij ≥ 0, and where ξ i is a convex piecewise-linear function. Then the minimizer u * = argmin can be found exactly using Algorithm 3.
Proof. First, to accommodate the ξ i (u i ) functions, suppose, without loss of generality, that each function is defined by a sequence of m i inflection points and an associated line slope θ ij between each, so Then, since ξ i (x) is linear, it is easy to show that at any point x, ξ i (x) can be expressed as the integral over sums of step functions; if c − i and c + i are the values of these step functions, The exact values of c + ij and c − ij can be found easily by solving a simple linear system, and they can satisfy c + ij + c − ij = 0 by adjusting the constant term in front of the integral; this ensures these pairwise terms satisfy the submodularity condition. Thus ξ i (·), can be encoded on the graph by creating m i auxiliary nodes u i,j fixed at reduction values b i1 , ..., b i,mi using the method described in section 7.0.2. Setting the overall weight of the edge using ensures that the cost of a cut at any region between b ij and b i,j+1 incurs cost Thus, in the final integral, the reduction level u i incurs the penalty const + ξ(u i ). The proof then follows identically to theorem 4.13 with some minor modifications -namely, for the auxiliary nodes, we are actually working with: but this easily converts to 1 {un+i≤β} .
As the L 2 fit terms on the auxiliary nodes become constant in the limit, the rest of the proof of theorem 4.13 goes through for this form as well, giving us the desired result.
In this chapter,

Conclusion
In this work, we outlined the theoretical connections between network flows, results in submodular optimization, and regularized regression problems such as the graph-guided LASSO. We rigorously established the submodular structure underlying the optimization of this problem. This motivated several novel algorithms, and extended some of the theory surrounding the minimum norm algorithm. We extended the existing theory of size-constrained submodular optimization first proposed by [Nagano et al., 2011] to the weighted case. This theoretical tool has several important consequences.
In the case of network flows, it gives us the ability to make nodes more or less affected by the optimization process. This opens the door to the general optimization problem given in theorem 8.1.
In part 2, we extend these results to develop a full treatment of the entire regularization path over all λ. In the second part of this work, we will explore further implications of this theory, namely a technique to algorithmically recycle solutions for network flows through the use of the structure presented here.

A Proofs
In this section, we give a number of the proofs needed for the previous theorems; they are given here for readability.
Proof of Theorem 4.1. To show the equivalence between P E (β) and P Q (β), it is sufficient to verify for arbitrary x ∈ {0, 1} n . To simplify notation, assume that E ij (x i , x j ) = 0 for all (i, j) / ∈ E. Then Thus the values of the objective function differ by a constant amount for all values x, so the set of minimizers is identical. Now, to verify the equivalence of P Q (β) and P N (β), let x i be the associated indicator vector, given by: It then suffices to verify that for all x ∈ {0, 1} n and corresponding cuts {s} ⊆ S ⊆ V : Again, the values of the objective functions differ by a constant amount for all x, proving the theorem.
Proof of Theorem 4.3. From [Fujishige, 2005], we know that f β is submodular if and only if, for all For convenience, assume that the β is absorbed into the diagonal elements of q ii . Now The theorem immediately follows; this is non-negative if q ij ≤ 0 ∀ i < j, and T, U, v can be easily chosen to make it negative if ∃ i, j such that q ij > 0.
Proof of theorem 4.4. For ease of notation, assume that βw i is absorbed into q ii . For ease of algebra, we first prove a simpler result, then show it immediately implies the desired one. Define (Recall that q ij ≤ 0.) Then, we show that the associated polymatroid P (f β ) can be represented by To show that P (f β ) is the base of f , we must show that ∀ y ∈ P (f β ) and ∀ U ⊆ V, i∈U y i ≤ f (U ) and that ∀ U ⊆ V, there exists a point y ∈ P (f β ) such that i∈U y i = f (U ) [Schrijver, 2003].
To show equation (160), we first prove that ∀ y ∈ P (f β ), i∈U y It remains to show that ∀ U ⊆ V, ∃ y ∈ P (f β ) such that i∈U y i = f (U ). For this, set w * ij = q ij i ∈ U, i < j 0 otherwise (164) Clearly, w * ∈ W. Setting y * to be the extreme point of P (f β ) using w = w * , we have, for all U ⊆ V, Thus (160) defines the polymatroid polytope of f β . Similarly, we show that the base of P (f β ) can be given by To show this, note that, trivially, B(f β ) ⊂ P (f β ). It remains to show that for all y ∈ B(f β ), i y i = f (V) for all w ∈ W: Finally, the theorem follows by replacing w ij with 1 2 (α ij − q ij ) in the above results and bringing the βw i outside of q ii and r(α).
Proof of theorem 4.13. First, we can write γ(u) as Now, define 1 {u≤β} as the cut vector of u at β, given by Now, assume that u i ∈ [−M, M ] for i = 1, 2, ..., n, where M is sufficiently large. This is reasonable, as λ −1 u T u dominates the optimization asymptotically. This enables us to write the terms in the above expression as  (176) However, the first term in the given approach maps directly to the unary q ii − β terms in the graph problem earlier, with the 1 {u≤β} term indexing the cut from source to node or from node to sink. Thus, for a given β, 1 {ui≤β} is a binary variable with 1 {ui≤β} indicating that the edge from source to node i is cut, incurring the appropriate cost in the overall function. Similarly, I 1 {ui≤β} = 1 {uj ≤β} is 1 for the range of β in which the edge connecting nodes i and j is cut; i.e. where the indicator variables of set membership differ. The cost of cutting this edge is given by 1 2 λw ij , and it is incurred by all points in the range of β between min(u i , u j ) and max(u i , u j ). Integrating this gives us 1 2 λw ij |u i − u j |. This establishes (176). Now, it remains to show that finding the minimum energy cut solution of the above problem at all β is equivalent to finding u * . For this step, define a set of functions mapping R to {0, 1} as: However, by corollary 4.10, we know that for any β 1 < β 2 , with x * (β 1 ) ∈ Argmin x∈{0,1} n h(β 1 , x) x * (β 2 ) ∈ Argmin x∈{0,1} n h(β 2 , x), we have the following implications for all i: if x * i (β 1 ) = 0, then x * i (β 2 ) = 0 (185) if x * i (β 2 ) = 1, then x * i (β 1 ) = 1.
This, however, is exactly the constraint implied in the set S . Thus (184) holds with equality: However, this minimization is exactly the problem solved earlier; namely, by theorem 4.9, we find the optimal set for each β. Thus, the net result of this operation is that u * = r(α * ). Letting M → ∞ completes the theorem.
Proof of lemma 5.1. ((I)) First, we establish that f w is indeed a submodular function defined on a distributed lattice. As given, f w is defined on D = T ∈ 2 Vw : T ∩ V ∈ D .
where (199) for all S ∈ D, proving the first part of the theorem. ((II)) To show (II), note that the membership of each item in V w \V exactly matches the membership of an item in V; as D is a distributed lattice, D w is also a distributed lattice that is closed under intersection and union. It remains to show that (94) is equivalent to (90). This follows directly from noting that for M sufficiently large, D w can be written as As we have already argued that in any optimal solution T * of f w (T ), the M terms cancel out, so any minimizer of f w (T ) must be in D w , proving (II). ((III)) Finally, in this case, for all {j, k} such that ∃ S ∈ D w with {j, k} ⊆ S, g jk (S) = 0, proving (95) and completing the proof.
Proof of theorem 6.2. ((I)) Recall that the definition of the polymatroid associated with the submodular function f w is defined as Now, by equation (95) in lemma 5.1, for all T ∈ D w . Furthermore, for all j ∈ K i , i ∈ T if and only if j ∈ T . Thus the condition and the condition x(V w ) = f w (V w ) is equivalent to This proves part (I) ((II)-(III)) Now, to prove part (II), consider minimizing y The objective ( ) in (211) above can then be expressed as ( ) = minimize i∈V h i (y, y(κ i )).
This, however, has an analytic solution, as the L 2 -norm above is minimized when y is equal on the set κ i , i.e. Thus This proves equation (100), and allows us to rewrite (211) as which is exactly identical to the optimization in equation (101). Given that y is constant on each κ i , we have already proved equation (102) by setting t = z i in equation (214).

A.1 Proof of Correctness for the Weighted Minimum Norm Problem
Before proving theorem 6.3, we must establish a lemma that shows convergence of the mapping in (103) under various perturbations. The reason that this lemma is needed is that the crux of the proof of theorem 6.3 involves showing that all positive w ∈ R n can be seen as the limit of a sequence of problems with integer w ∈ Z n . The crux of this lemma is found in Wets [2003], where a convenient theorem shows that mappings of the type (103) are locally Lipshitz-continuous with respect to perturbations of the objective function. This continuity is sufficient to ensure that the limiting argument employed in the proof of theorem 6.3 is valid.