Runtime Guarantees for Regression Problems

Introduction

The problem of recovering a discrete, clear signal from noisy data is an important problem in signal processing. One general approach to this problem is to formulate an objective based on required properties of the answer, and then return its minimizer via optimization algorithms. The power of this method was first demonstrated in image denoising, where the total variation minimization approach by Rudin, Osher and Fatemi [ROF92] had much success. More recent works on sparse recovery led to the theory of compressed sensing [Can06], which includes approaches such as the least absolute shrinkage and selection operator (LASSO) objective due to Tibshirani [Tib96]. These objective functions have proven to be immensely powerful tools, applicable to problems in signal processing, statistics, and computer vision. In the most general form, given vector y and a matrix A, one seeks to minimize:

It can be shown to be equivalent to the following by introducing a Lagrangian multiplier, $\lambda$:

Many of the algorithms used to minimize the LASSO objective in practice are first order methods [NES07, BCG11], which updates a sequence of solutions using well-defined vectors related to the gradient. These methods are guaranteed to converge well when the matrix A is “well-structured”. The formal definition of this well-structuredness is closely related to the conditions required by the guarantees given in the compressed sensing literature [Can06] for the recovery of a sparse signal. As a result, these methods perform very well on problems where theoretical guarantees for solution quality are known. This good performance, combined with the simplicity of implementation, makes these algorithms the method of choice for most problems.

However, LASSO type approaches have also been successfully applied to larger classes of problems. This has in turn led to the use of these algorithms on a much wider variety of problem instances. An important case is image denoising, where works on LASSO-type objectives predates the compressed sensing literature [ROF92]. The matrices involved here are based on the connectivity of the underlying pixel structure, which is often a $\sqrt{n}\times\sqrt{n}$ square mesh. Even in a unweighted setting, these matrices tend to be ill-conditioned. In addition, the emergence of non-local formulations that can connect arbitrary pairs of vertices in the graph also highlights the need to handle problems that are traditionally considered ill-conditioned. We show in Appendix A that the broadest definition of LASSO problems include well-studied problems from algorithmic graph theory:

Both the $s$-$t$ shortest path and $s$-$t$ minimum cut problems in undirected graphs can be solved by minimizing a LASSO objective.

Although linear time algorithms for unweighted shortest path are known, finding efficient parallel algorithms for this has been a long-standing open problem. The current state of the art parallel algorithm for finding $1+\epsilon$ approximate solutions, due to Cohen [Coh00], is quite involved. Furthermore, as the reductions done in Lemma A.1 are readily parallelizable, an efficient algorithm for LASSO minimization would also lead to an efficient parallel shortest path algorithm. This suggests that algorithms for minimizing LASSO objectives, where each iteration involve simple, parallelizable operations, are also difficult. Finding a minimum $s$-$t$ cut with nearly-linear running time is also a long standing open question in algorithm design. In fact, there are known hard instances where many algorithms do exhibit their worst case behavior [JM93]. The difficulty of these problems and the non-linear nature of the objective are two of the main challenges in obtaining fast run time guarantees for grouped least squares minimization.

Previous run time guarantees for minimizing LASSO objectives rely on general convex optimization routines [BV04], which take at least $\Omega(n^{2})$ time. As the resolution of images are typically at least $256\times 256$, this running time is prohibitive. As a result, when processing image streams or videos in real time, gradient descent or filtering based approaches are typically used due to time constraints, often at the cost of solution quality. The continuing increase in problem instance size, due to higher resolution of streaming videos, or 3D medical images with billions of voxels, makes the study of faster algorithms an increasingly important question.

While the connection between LASSO and graph problems gives us reasons to believe that the difficulty of graph problems also exists in minimizing LASSO objectives, it also suggests that techniques from algorithmic graph theory can be brought to bear. To this end, we draw upon recent developments in algorithms for maximum flow [CKM+11] and minimum cost flow [DS08]. We show that relatively direct modifications of these algorithms allows us to solve a generalization of most LASSO objectives, which we term the grouped least squares problem. Our algorithm is similar to convex optimization algorithms in that each iteration of it solves a quadratic minimization problem, which is equivalent to solving a linear system. The speedup over previous algorithms come from the existence of much faster solvers for graph related linear systems [ST06], although our approaches are also applicable to situations involving other underlying quadratic minimization problems.

The organization of this paper is as follows: In Section 2 we provide a unified optimization problem that encompasses LASSO, fused LASSO, and grouped LASSO. We then discuss known applications of the grouped least squares minimization in Section 3 and existing approaches in Section 4. In Section 5 we give two algorithms that use solving quadratic minimization problems as underlying routines: An approximate algorithm based on the maximum flow algorithm of Christiano et al. [CKM+11], and an almost-exact algorithm that rely on interior point algorithms.

Background and Formulations

The formulation of our main problem is motivated by the total variation objective from image denoising. This objective has its origin in the seminal work by Mumford and Shah [MS89]. There are two conflicting goals in recovering a smooth image from a noisy one, namely that it must be close to the original image, while having very little noise. The Mumford-Shah function models the second constraint by imposing penalties for neighboring pixels that differ significantly. These terms decrease with the removal of local distortions, offsetting the higher cost of moving further away from the input. However, the minimization of this functional is computationally difficult and subsequent works focused on minimizing functions that are close to it.

The total variation objective is defined for a discrete, pixel representation of the image and measures noise using a smoothness term calculated from differences between neighboring pixels. This objective leads naturally to a graph $G=(V,E)$ corresponding to the image with pixels. The original (noisy) image is given as a vertex labeling s, while the goal of the optimization problem is to recover the ‘true’ image x, which is another set of vertex labels. The requirement of x being close to s is quantified by $||\textbf{x}-\textbf{s}||_{2}^{2}$, specifically summing over the squares of what we identify as noise while the smoothness term is a sum over absolute values of difference between adjacent pixels’ labels:

This objective can be viewed as an instance of the fused LASSO objective [TSR+05]. As the orientation of the underlying pixel grid is artificially imposed by the camera, this method can introduce rotational bias in its output. One way to correct this bias is to group the differences of each pixel with its 4 neighbors, giving terms of the form:

where $v$ and $w$ are the horizontal and vertical neighbor of $u$.

Our generalization of these objectives is based on the key observation that $\sqrt{(x_{u}-x_{v})^{2}+(x_{u}-x_{w})^{2}}$ and $|x_{u}-x_{v}|$ are both $L_{2}$ norms of a vector containing differences of values of adjacent pixels. Each such difference can be viewed as an edge in the underlying graph, and the grouping gives a natural partition of the edges into disjoint sets $S_{1}\dots S_{k}$:

It can be checked that when each $S_{i}$ contain a single edge, this formulation is identical to the objective in Equation 2.3 since $\sqrt{(x_{u}-x_{v})^{2}}=|x_{u}-x_{v}|$. Then all the terms can be written as quadratic positive semi-definite terms involving x and $\textbf{s}_{i}$, where we now allow for different fixed values in the groups. Specifically when given symmetric positive semidefinite (PSD) matrices $\textbf{L}_{0}\ldots\textbf{L}_{k}$, the objective can be rewritten as:

If we use $||\cdot||_{\textbf{L}_{i}}$ to denote the norm induced by the PSD matrix $\textbf{L}_{i}$, each of the terms can be written as $||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}$. To make the first term resemble the other terms in our objective, we will take a square root of it – as we prove in Appendix D, algorithms that give exact minimizers for this variant still captures the original version of the problem. This simplification allows us to define our main problem:

Input: $n\times n$ matrices $\textbf{L}_{1}\ldots\textbf{L}_{k}$ and fixed potentials $\textbf{s}_{1}\ldots\textbf{s}_{k}\in\Re^{n}$.

Note that this objective allows for the usual definition of LASSO involving terms of the form $|\textbf{x}_{u}|$ by having one group for each such variable with $\textbf{s}=\textbf{0}$. It is also related to group LASSO [YL06], which incorporates similar assumptions about closer dependencies among some of the terms. To our knowledge grouping has not been studied in conjunction with fused LASSO, although many problems such as the ones listed in Section 3 require this generalization.

Our algorithmic approach to the group least squares problem crucially depends on solving a related quadratic minimization problem. Specifically, we solve linear systems involving a weighted combination of the $\textbf{L}_{i}$ matrices. Let $\textbf{w}_{1}\ldots\textbf{w}_{k}\in\Re^{+}$ denote weights, where $\textbf{w}_{i}$ is the weight on the $i$th group. Then the quadratic minimization problem that we consider is:

We will use $\text{OPT2}(\textbf{w})$ to denote the minimum value that is attainable. This minimizer, x, can be obtained using the following Lemma:

$\mathcal{OBJ}2(\textbf{x},\textbf{w})$ is minimized for x such that

Applications

A variety of problems ranging from computer vision to statistics can be formulated as grouped least squares. We describe some of them below, starting with classical problems from image processing.

As mentioned in Section 1, one of the earliest applications of these objectives was in the context of image processing. More commonly known as total variation minimization in this setting [CS05], various variants of the objective have been proposed with the anisotropic objective the same as Equation 2.3 and the isotropic objective being the one shown in Equation 2.5.

2 Denoising with Multiple Colors

Most works on image denoising deals with images where each pixel is described using a single number corresponding to its intensity. A natural extension would be to colored images, where each pixel has a set of $c$ attributes (in the RGB case, $c=3$). One possible analogue of $|x_{i}-x_{j}|$ in this case would be $||x_{i}-x_{j}||_{2}$, and this modification can be incorporated by replacing a cluster involving a single edge with clusters over the $c$ edges between the corresponding pixels.

This type of approach can be viewed as an instance image reconstruction algorithms using Markov random fields. Instead of labeling each vertex with a single attribute, a set of $c$ attributes are used instead and the correlation between vertices is represented using arbitrary PSD matrices. It’s worth remarking that when such matrices have bounded condition number, it was shown in [KMP12] that the resulting least squares problem can still be solved efficiently by preconditioning with SDD matrices, yielding a similar overall running time.

3 Poisson Image Editing

The Poisson Image Editing method of Perez, Gangnet and Blake [PGB03] is a very popular method for image blending. This method aims to minimize the difference between the gradient of the image and a guidance field vector v. We show here that the grouped least square problem can be used for minimizing objectives from this framework. The objective function given in equation (6) of [PGB03]

This term can be rewritten as $((x_{p}-x_{q})-(v_{pq}-0))^{2}$. So if we let $\textbf{s}_{i}$ be the vector where $s_{i,p}=v_{pq}$ and $s_{i,q}=0$, and $\textbf{L}_{i}$ be the graph Laplacian for the edge connecting $p$ and $q$, then the term equals to $||\textbf{x}-\textbf{s}_{i}||^{2}_{\textbf{L}_{i}}$. The other terms on the boundary will have $x_{q}$ as a constant, leading to terms of the form $||x_{i,p}-s_{i,p}||^{2}_{2}$ where $s_{i,p}=x_{q}$. Therefore the discrete Poisson problem of minimizing the sum of these squares is an instance of the quadratic minimization problem as described in Section 2.1. Perez et al. in Section 2 of their paper observed that these linear systems are sparse, symmetric and positive definite. We make the additional observation here that the systems involved are also symmetric diagonally dominant. The use of the grouped least squares framework also allows the possibility of augmenting these objectives with additional $L_{1}$ or $L_{2}$ terms.

4 Clustering

Hocking et al. [HVBJ11] recently studied an approach for clustering points in $d$ dimensional space. Given a set of points $x_{1}\ldots x_{n}\in\Re^{d}$, one method that they proposed is the minimization of the following objective function:

Previous Algorithmic Results

Due to the importance of optimization problems motivated by LASSO there has been much work on efficient algorithms for them. We briefly describe some of the previous approaches for LASSO minimization below.

2 Graph Cuts

It’s worth mentioning that both of these algorithms requires extracting the minimum cut in order to construct the problems for subsequent iterations. As a result, it’s not clear whether recent advances on fast approximations of maximum flow and minimum $s$-$t$ cuts [CKM+11] can be used as a black box with these algorithms. Extending this approach to the non-linear isotropic objective also appears to be difficult.

3 Iterative Reweighted Least Squares

An approach similar to convex optimization methods, but has much better observed rates of convergence is the iterative reweighted least squares (IRLS) method. This method does a much more aggressive adjustment each iteration and to give good performances in practice [WR07].

4 First Order Methods

The method of choice in practice are first order methods such as [NES07, BCG11]. Theoretically these methods are known to converge rapidly when the objective function satisfies certain Lipschitz conditions. Many of the more recent works on first order methods focus on lowering the dependency of $\epsilon$ under these conditions. As discussed in Section 1 and Appendix A, this direction can be considered orthogonal to our guarantees as the grouped least squares problem is a significantly more general formulation.

Solving Grouped Least Squares Using Quadratic Minimizations

In this section, we show two algorithms for the grouped least squares problem based on direct adaptations of state-of-the-art algorithms for maximum flow [CKM+11] and minimum cost flow [DS08]. Our guarantees can be viewed as reductions to the quadratic minimization problems described in Section 2.1. As a result, they imply efficient algorithms for problems where fast algorithms are known for the corresponding least squares problems. The analyses of these algorithms are intricate, but mostly follow the presentations in [CKM+11, DS08, BV04]. They’re presented in Appendices B and C.

Our first algorithm is based on the electrical flow based maximum flow and minimum cut algorithm by Christiano et al. [CKM+11]. Recall that the minimum $s$-$t$ cut problem - equivalent to an $L_{1}$-minimization problem - is a special case of the grouped least squares problem where each edge belongs to its own group(i.e., $k=m$). As a result, it’s natural to extend the approach of [CKM+11] to the whole spectrum of values of $k$ by treating each group as an ’edge’.

One view of the cut algorithm from [CKM+11] is that it places a weight on each group, and minimizes a quadratic, or $L_{2}^{2}$ problem involving terms of the from $\frac{1}{\textbf{w}_{i}}||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}^{2}$. These weights are adjusted using the multiplicative weights update framework [AHK05, LW94] based on the energy of each group. The terms $||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}$ and $\frac{1}{\textbf{w}_{i}}||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}^{2}$ are equal when $\textbf{w}_{i}=||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}$. Therefore, one natural view of this routine is that it gradually adjusts the weights to become a scaled copy of $||\bar{\textbf{x}}-\textbf{s}_{i}||_{\textbf{L}_{i}}$. This leads to a simplification of the Christiano et al. algorithm, whose update requires a flow obtained from the dual of the quadratic minimization problem. Pseudocode of the algorithm is shown in Algorithm 1.

The main difficulty of analyzing this algorithm is that the analysis of minimum $s$-$t$ cut algorithm of [CKM+11] relies strongly on the existence of a solution where x is either or $1$. Our analysis extends this potential function into the fractional setting via. a function based on the Kulback-Liebler (KL) divergence [KL51]. To our knowledge the use of this potential with multiplicative weights was first introduced by Freund and Schapire [FS99], and is common in learning theory. This function can be viewed as measuring the KL-divergence between $\textbf{w}_{i}^{(t)}$ and $||\bar{\textbf{x}}-\textbf{s}_{i}||_{\textbf{L}_{i}}$ over all groups, where $\bar{\textbf{x}}$ an optimum solution to the grouped least squares problem. Formally, the KL divergence between these two vectors is:

Subtracting the constant term given by $\sum_{i}||\bar{\textbf{x}}-\textbf{s}_{i}||_{\textbf{L}_{i}^{(t)}}\log(||\bar{\textbf{x}}-\textbf{s}_{i}||_{\textbf{L}_{i}^{(t)}})$ and multiplying by $-1/\text{OPT}$ gives us our key potential function, $\nu^{(t)}$:

It’s worth noting that in the case of the cut algorithm, this function is identical to the potential function used in [CKM+11]. We show the convergence of our algorithm by proving that if the solution produced in some iteration is far from optimal, $\nu^{(t)}$ increases substantially. Upper bounding it with a term related to the sum of weights, $\mu^{(t)}$ allows us to prove convergence. The full proof is given in Appendix B.

To simplify the analysis, we assume that the guess that we’re trying to solve the decision problem on, OPT, all entries of s, and spectrum of $\textbf{L}_{i}$ are polynomially bounded in $n$. That is, there exist some constant $d$ such that $-n^{d}\leq\textbf{s}_{i,u}\leq n^{d}$ and $n^{-d}\textbf{I}\preceq\sum_{i}\textbf{L}_{i}\preceq n^{d}\textbf{I}$ where $\textbf{A}\preceq\textbf{B}$ means $\textbf{B}-\textbf{A}$ is PSD. Some of these assumptions can be relaxed via. analyses similar to Section 2 of [CKM+11].

On input of an instance of $\mathcal{OBJ}$ with edges partitioned into $k$ sets. If all parameters polynomially bounded between $n^{-d}$ and $n^{d}$, running ApproxGroupedLeastSquares with $\rho=2k^{1/3}\epsilon^{-2/3}$ returns a solution x with such that $\mathcal{OBJ}(\textbf{x})\leq\max\{(1+10\epsilon)\text{OPT},n^{-d}\}$ where OPT is the value of the optimum solution.

The additive $n^{-d}$ case is included to deal with the case where $\text{OPT}=0$, or is close to it. We believe it should be also possible to handle this case via. condition numbers restrictions on $\sum_{i}\textbf{L}_{i}$.

2 Almost-Exact Algorithm

There are various ways to solve the grouped least squares problem using interior point algorithms. We follow the log-barrier method, as presented in Boyd and Vandenberghe [BV04] here for simplicity. This formulation defines one extra variable $y_{i}$ for each group and enforces $y_{i}\geq||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}$ using the barrier function $\phi_{i}(\textbf{x},y_{i})=\log(y_{i}^{2}-||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}^{2})$. Minimizing $t\cdot(\sum_{i}y_{i})$ for gradually increasing values of $t$ gives the following sequence of functions to minimize:

Various interior point algorithms have been proposed, one commonality that they have is finding an update direction by solving a linear system. The iteration guarantees for recovering almost exact solution can be characterized as follows:

These systems are examined in detail in Appendix C. Since the $t\sum_{i}y_{i}$ term is linear, it can be omitted from the Hessian, leaving $\sum_{i}\nabla^{2}\phi(\textbf{x},y_{i})$. We then check that the barrier term $y_{i}$ creates a low rank perturbation to the $\textbf{L}_{i}$ term, which is the Hessian for $||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}^{2}$. By taking Schur complements and applying the Sherman-Morrison-Woodbury identity on inverses for low rank perturbations, we arrive at the following observation in Appendix C:

Suppose there is an algorithm for solving linear systems of the form $\sum_{i}\alpha_{i}\textbf{L}_{i}$ in $T(n,m)$ time where $m$. For any choice of $\textbf{x},\textbf{y}$, a linear system involving the Hessian of $\phi(\textbf{x},\textbf{y})$, $\nabla^{2}\phi(\textbf{x},\textbf{y})$ can be solved in $O(k^{\omega}+kT(n,m)+k^{2}n)$ time.

Experimental Results

We performed a series of experiments using the approximate algorithm described in Section 5.1. The SDD linear systems that arise in the quadratic minimization problems were solved using the combinatorial multigrid (CMG) solver [KM09, KMT09]. One side observation from our experiments is that for the sparse SDD linear systems that arise from image processing, the CMG solver yields good results both in accuracy and running time.

Total Variational Denoising is the concept of applying Total Variational Minimization as denoising process. This was pioneered by Rudin, Osher and Fatemi [ROF92] and is commonly known as the ROF image model [CS05]. Our algorithm from Section 5.1 yields a simple way to solve the ROF model and most of its variants. In Figure 1, we present a simple denoising experiment using the standard image processing data set, ‘Leena’. The main goal of the experiment is to show that our algorithm is competitive in terms of accuracy, while having running times comparable to first-order methods. On a $512\times 512$ grayscale image, we introduce Additive White Gaussian Noise (AWGN) at a measured Signal to Noise Ratio (SNR) of 2. AWGN is the most common noise model in photon capturing sensors from consumer cameras to space satellites systems. Error is measured as the intensity difference from the original, uncorrupted image summed over all pixels

Our experiments were conducted on a single core 64-bit Intel(R) Xeon(R) E5440 CPU @ 2.83GHz. The non-solver portion of the algorithm was implemented in Matlab(R). On images of size $256\times 256$, $512\times 512$ and $1024\times 1024$, the average running times are $2.31$, $9.70$ and $47.61$ seconds respectively.

It’s worth noting is that on average 45% of the total running time is from solving the SDD linear systems using the CMG solver. The rest is due to reweighting edges in Matlab, and should be handled as part of the CMG solver routine in more optimized versions. More importantly, in all of our experiments the weights are observed to converge in under 15 iterations, even for larger images of size up to $3000\times 3000$. This is much better than the guarantees given for either of our algorithms in Section 5.

2 Image Processing

As exemplified by the denoising with colors application in Section 3.2, the grouped least squares framework has great flexibility in expressing image processing tasks. We applied our denoising algorithm to Optical Coherence Tomography (OCT) images of the retina as a preprocessing step for segmentation. Here the key is to preserve the sharpness between the nerve fiber layers and this is achieve by using a $L_{1}$ regularization term.

Variations of this formulation allows one to model a large variety of established image preprocessing applications. For example, Gaussian blurred images can be obtained using a $L_{2}$ penalty term. This simulates camera zoom and is similar to the preprocessing step in the popular Scale Invariant Feature Transform (SIFT) algorithm. By mixing and matching penalty terms on the groups, we can preserve global features while favoring the removal of small artifacts that are often the result of sensor noise. Our approaches also extend naturally to multichannel images, (RGB or multi-spectral), with little modification to the underlying algorithm.

An example of Poisson Image Editing mentioned in Section 3.3 is shown in Figure 2. The specific application is seamless cloning as described in Section 3 of [PGB03], which aims to insert complex objects into another image. Given two images, they are blended by solving the discrete poisson equation based on a mix of their gradients and boundary values. We also added $L_{2}$ constraints on different parts of the image to give a smoother result.

Remarks

We believe that the ability of our algorithm to encompass many of the current image processing algorithms represents a major advantage in practice. It allows the use of a common data structure (the underlying graph) and subroutine (the SDD solver) for many different tasks in the image processing pipeline. Theoretically, this feature is also interesting as it represents a smooth interpolation between $L_{2}$ and $L_{1}$ problems.

The performances of our algorithms depend on $k$, which is the number of groups in the formulation given in Definition 2.1. , it gives them provable runtime guarantees. Two settings of $k$ are helpful for comparison to previous works. When $k=1$, the problem becomes the electrical flow problem, and the running time of both algorithms are similar to directly solving the linear system. This is also the case when there is a small (constant) number of groups. The other extremum is when each edge belongs to its own group, aka. $k=m$. Here our approximate algorithm is the same as the minimum $s$-$t$ cut algorithm given in [CKM+11], but our analysis for our almost exact algorithm gives a much worse running time. This is due to the interior point algorithm generating much more complicated linear systems, and actually occur when most groups contain a small number of edges. As a result, finding faster algorithms with better error guarantees for problems with intermediate to large values of $k$ is an interesting direction for future work. Also, preliminary experimental results such as the ones from Section 6 show that more aggressive reweightings of edges lead to much faster convergence than what we showed for our two algorithms. Therefore a suite of examples where the theoretical guarantees are tight would also give a better understanding of the interplay between multiplicative updates and quadratic minimization.

One other consequence of this dependency on $k$ is that although the problem with smaller number of groups is no longer captured by linear optimization, the minimum $s$-$t$ cut problem - that still falls within the framework of linear optimization is in some sense the hardest problem in this class. Therefore we believe that the grouped least squares problem is a natural interpolation between the $L_{1}$ and $L_{2}$ problems, and has potential as an intermediate subroutine in graph algorithms.

Acknowledgements

The authors would like to thank Jerome Darbon, Stanley Osher, Aarti Singh and Ryan Tibshirani for pointing them to works that are relevant to the grouped least squares framework, and also an anonymous reviewer of a previous submission of this paper for pointing out an alternate view of the proof of Theorem 5.1.

References

Appendix A Proofs about Graph Problems as Minimizing LASSO Objectives

In this section we give formal proofs that show the shortest path problem is an instance of LASSO and the minimum cut problem is an instance of fused-LASSO.

It’s worth noting that our proofs do not guarantee that the answers returned are a single path or cut. In fact, when multiple solutions have the same value it’s possible for our algorithm to return a linear combination of them. However, we can ensure that the optimum solution is unique using the Isolation Lemma of Mulmuley, Vazarani and Vazarani [MVV87] while only incurring polynomial increase in edge lengths/weights. This analysis is similar to the one in Section 3.5 of [DS08] for finding a unique minimum cost flow, and is omitted here.

We prove the two claims in Fact 1.1 about shortest path and minimum cut separately in Lemmas A.1 and A.2.

Given a $s$-$t$ shortest path instance in an undirected graph where edge lengths $\textbf{l}:E\rightarrow\Re^{+}$ are integers between $1$ and $n^{d}$. There is a LASSO minimization instance where all entries are bounded by $n^{O(d)}$ such that the value of the LASSO minimizer is within $1$ of the optimum answer.

Proof Our reductions rely crucially on the edge-vertex incidence matrix, which we denote using B. Entries of this matrix are defined as follows:

We first show the reduction for shortest path. Then a path from $s$ to $t$ corresponds to a flow value assigned to all edges, $\textbf{f}:E\rightarrow\Re$ such that $\textbf{B}^{T}\textbf{f}=\chi_{s,t}$. If we have another flow $\textbf{f}^{\prime}$ corresponding to any path from $s$ to $t$, then this constraint can be written as:

The first constraint is closer to the classical LASSO problem while the last one is within our definition of grouped least squares problem. The length of the path can then be written as $\sum_{e}l_{e}|f_{e}|$. Weighting these two terms together gives:

Where the maximum entry is bounded by $\max\{n^{d},n^{2}\lambda\}$. Clearly its objective is less than the length of the shortest path, let this solution be $\bar{\textbf{f}}$. Then since the total objective is at most $n^{d+1}$, we have that the maximum deviation between $\textbf{B}^{T}\bar{\textbf{f}}$ and $\chi_{s,t}$ is at most $n^{d+1}/\lambda$. Then given a spanning tree, each of these deviations can be routed to $s$ or $t$ at a cost of at most $n^{d+1}$ per unit of flow. Therefore we can obtain $\textbf{f}^{\prime}$ such that $\textbf{B}^{T}\textbf{f}^{\prime}=\chi_{s,t}$ whose objective is bigger by at most $n^{2d+2}/\lambda$. Therefore setting $\lambda=n^{2d+2}$ guarantees that our objective is within $1$ of the length of the shortest path, while the maximum entry in the problem is bounded by $n^{O(d)}$. $\blacksquare$

We now turn our attention to the minimum cut problem, which can be formulated as finding a vertex labeling $\textbf{x}^{(vert)}$ where $\textbf{x}^{(vert)}_{s}=0$, $\textbf{x}^{(vert)}_{t}=1$ and the size of the cut, $\sum_{uv\in E}|\textbf{x}^{(vert)}_{u}-\textbf{x}^{(vert)}_{v}|$. Since the $L_{1}$ term in the objective can incorporate single variables, we use an additional vector $\textbf{x}^{(edge)}$ to indicate differences along edges. The minimum cut problem then becomes minimizing $|\textbf{x}^{(edge)}|$ subject to the constraint that $\textbf{x}^{(edge)}=\textbf{B}^{\prime}\textbf{x}^{(vert)^{\prime}}+\textbf{B}\chi_{t}$, where $\textbf{B}^{\prime}$ and $\textbf{x}^{(vert)^{\prime}}$ are restricted to vertices other than $s$ and $t$ and $\chi_{t}$ is the indicator vector that’s $1$ on $t$. The equality constraints can be handled similar to the shortest path problem by increasing the weights on the first term. One other issue is that $|\textbf{x}^{(vert)^{\prime}}|$ also appears in the objective term, and we handle this by scaling down $\textbf{x}^{(vert)^{\prime}}$, or equivalently scaling up $\textbf{B}^{\prime}$.

Given a $s$-$t$ minimum cut instance in an undirected graph where edge weights $\textbf{w}:E\rightarrow\Re^{+}$ are integers between $1$ and $n^{d}$. There is a LASSO minimization instance where all entries are bounded by $n^{O(d)}$ such that the value of the LASSO minimizer is within $1$ of the minimum cut.

Consider the following objective, where $\lambda_{1}$ and $\lambda_{2}$ are set to $n^{d+3}$ and $n^{2}$:

Let $\bar{\textbf{x}}$ be an optimum vertex labelling, then setting $\textbf{x}^{(vert^{\prime})}$ to the restriction of $n^{-2}\bar{\textbf{x}}$ on vertices other than $s$ and $t$ and $\textbf{x}^{(edge)}$ to $\textbf{B}\textbf{x}^{(vert)}$ makes the first term . Since each entry of $\bar{\textbf{x}}$ is between and $1$, the additive increase caused by $|\textbf{x}^{(vert^{\prime})}|_{1}$ is at most $1/n$. Therefore this objective’s optimum is at most $1/n$ more than the size of the minimum cut.

For the other direction, consider any solution $\textbf{x}^{(vert)^{\prime}},\textbf{x}^{(edge)^{\prime}}$ whose objective is at most $1/n$ more than the size of the minimum cut. Since the edge weights are at most $n^{d}$ and $s$ has degree at most $n$, the total objective is at most $n^{d+1}+1$. This gives:

Therefore changing $\textbf{x}^{(edge)}$ to $\lambda_{2}\textbf{B}^{\prime}\textbf{x}^{(vert)^{\prime}}+\textbf{B}\chi_{t}$ increases the objective by at most $n^{-1/2}<1$. This gives a cut with weight within $1$ of the objective value and completes the proof. $\blacksquare$

We can also show a more direct connection between the minimum cut problem and the fused LASSO objective, where each absolute value term may contain a linear combination of variables. This formulation is closer to the total variation objective, and is also an instance of the problem formulated in Definition 2.1 with each edge in a group.

The minimum cut problem in undirected graphs can be written as an instance of the fused LASSO objective.

Proof Given a graph $G=(V,E)$ and edge weights cost, the problem can be formulated as labeling the vertices of $V\setminus\{s,t\}$ in order to minimize:

Appendix B Multiplicative Weights Based Approximate Algorithm

In this section we show that the approximate algorithm described in Section 5.1 finds a solution close to the optimum. We first show that if $\lambda^{(t)}$ as defined on Line 6 of Algorithm 1 is an upper bound for $\mathcal{OBJ}(\textbf{x}^{(t)}$. This is crucial in our use of it as the normalizing factor in our update step on Line 7.

Proof By the Cauchy-Schwarz inequality we have:

Taking square roots of both sides completes the proof. $\blacksquare$

At a high level, the algorithm assigns weights $\textbf{w}_{i}$ for each group, and iteratively reweighs them for $N$ iterations. Recall that our key potential functions are $\mu^{(t)}$ which is the sum of weights of all groups, and:

Where $\bar{\textbf{x}}$ is a solution such that $\mathcal{OBJ}(\bar{\textbf{x}})=\text{OPT}$. We will show that if $\mathcal{OBJ}(\textbf{x}^{(t)})$, or in turn $\lambda^{(t)}$ is large, then $\nu^{(t)}$ increases at a rate substantially faster than $\mu^{(t)}$. These bounds, and the relation between $\mu^{(t)}$ and $\nu^{(t)}$ is summarized below:

If in iteration $t$, $\lambda^{(t)}\geq(1+10\epsilon)\text{OPT}$ and $||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}\leq\rho\frac{\textbf{w}_{i}^{(t-1)}}{\mu^{(t-1)}}\lambda^{(t)}$ for all groups $i$, then:

The relationship between the upper and lower potentials can be established using the fact that $\textbf{w}_{i}$ is non-negative:

Part 2 follows directly from the local behavior of the $\log$ function:

Proof of Lemma B.2, Part 2: The update rules gives:

Using the fact that $1+x\leq\exp(x)$ when $x\geq 0$ we get:

Taking logs of both sides gives Equation 2.27.

This upper bound on the value of $\mu^{t}$ also allows us to show that the balancing rule keeps the $w^{t}_{i}$s reasonably balanced within a factor of $k$ of each other. The following corollary can also be obtained.

The weights at iteration $t$ satisfy $\textbf{w}^{(t)}_{i}\geq\frac{\epsilon}{k}\mu^{(t)}$.

The proof is by induction on $t$. When $t=0$ we have $\textbf{w}^{(0)}_{i}=1$, $\mu^{(0)}=k$ and the claim follows from $\frac{\epsilon}{k}k=\epsilon<1$. When $t>1$, we have:

The proof of Part 3 is the key part of our analysis. The first order change of $\nu^{(t)}$ is written as a sum of products of $L_{2}$ norms, which we analyze via. the fact that $\textbf{x}^{(t)}$ is the solution of a linear system from the quadratic minimization problem.

We make use of the following known fact about the behavior of the log function around $1$:

If $0\leq x\leq\epsilon$, then $\log(1+x)\geq(1-\epsilon)x$.

Since $\textbf{L}_{i}$ forms a P.S.D norm, by the Cauchy-Schwarz inequality we have:

Recall from Lemma 2.2 that since $\textbf{x}^{(t)}$ is the minimizer to $\mathcal{OBJ}2(\textbf{w}^{(t-1)})$, we have

Substituting this into Equation 2.32 gives:

Since the iteration count largely depends on $\rho$, it suffices to provide bounds for $\rho$ over all the iterations. The proof makes use of the following lemma about the properties of electrical flows, which describes the behavior of modifying the weights of a group $S_{i}$ that has a large contribution to the total energy. It can be viewed as a multiple-edge version of Lemma 2.6 of [CKM+11].

Assume that $\epsilon^{2}\rho^{2}<1/10k$ and $\epsilon<0.01$ and let $\textbf{x}^{(t-1)}$ be the minimizer for $\mathcal{OBJ}2(\textbf{w}^{(t-1)})$. Suppose there is a group $i$ such that $||\textbf{x}^{(t-1)}-\textbf{s}_{i}||_{\textbf{L}_{i}}\geq\rho\frac{\textbf{w}_{i}^{(t-1)}}{\mu^{(t-1)}}\lambda^{(t)}$, then

We first show that group $i$ contributes a significant portion to $\mathcal{OBJ}2(\textbf{w}^{(t-1)},\textbf{x}^{(t-1)})$. Squaring both sides of the given condition gives:

Also, by the update rule we have $\textbf{w}^{(t)}_{i}\geq(1+\epsilon)\textbf{w}^{(t-1)}_{i}$ and $\textbf{w}^{(t)}_{j}\geq\textbf{w}^{(t-1)}_{j}$ for all $1\leq j\leq k$. So we have:

This means the value of the quadratic minimization problem can be used as a second potential function. We first show that it’s monotonic and establish rough bounds for it.

$\text{OPT2}(\textbf{w}^{(0)})\leq n^{3d}$ and $\text{OPT2}(\textbf{w}^{(t)})$ is monotonically decreasing in $t$.

Proof By the assumption that the input is polynomially bounded we have that all entries of s are at most $n^{d}$ and $\textbf{L}_{i}\preceq n^{d}\textbf{I}$. Setting $\textbf{x}_{u}=0$ gives $||\textbf{x}-\textbf{s}_{i}||_{2}\leq n^{d+1}$. Combining this with the spectrum bound then gives $||\textbf{x}-\textbf{s}_{i}||_{\textbf{L}_{i}}\leq n^{2d+1}$. Summing over all the groups gives the upper bound.

The monotonicity of $\text{OPT2}(\textbf{w}^{(t)})$ follows from the fact that all weights are decreasing. $\blacksquare$

Combining this with the fact that $\text{OPT2}(\textbf{w}^{(N)})$ is not low enough for termination gives our bound on the total iteration count.

Proof of Theorem 5.1: The proof is by contradiction. Suppose otherwise, since $\mathcal{OBJ}(\textbf{x}^{(N)})>\epsilon$ we have:

Which combined with $\text{OPT2}(\textbf{w}^{(0)})\leq n^{3d}$ from Lemma B.6 gives:

Combining with Lemma B.5 implies that the number of iterations where $||\textbf{x}^{(t-1)}-\textbf{s}_{i}||_{\textbf{L}_{i}}\geq\rho\frac{\textbf{w}_{i}^{(t-1)}}{\mu^{(t-1)}}\lambda^{(t)}$ for $i$ is at most:

This means that we have $||\textbf{x}^{(t-1)}-\textbf{s}_{i}||_{\textbf{L}_{i}}\leq\rho\frac{\textbf{w}_{i}^{(t-1)}}{\mu^{(t-1)}}\lambda^{(t)}$ for all $1\leq i\leq k$ for at least $(1-\epsilon)N$ iterations and therefore by Lemma B.2 Part 3:

Appendix C Solving Linear Systems from Interior Point Algorithms

In this section we take a closer look at the linear system solves required by Lemma 5.2, specifically the Hessians of the objective given in Equation 5.10. We show that for small values $k$, having access to an efficient subroutines for solving the quadratic minimization problems gives improvements over general interior-point algorithms.

We first consider the barrier function corresponding to each group, $\phi(\textbf{x},y_{i})$. Its gradient is:

Since the variable $y_{i}$ only appears in $\phi(\textbf{x},y_{i})$, we may use partial Cholesky factorization to arrive at a linear system without it. The $n\times n$ matrix that we obtain is:

Note that by the Cauchy-Schwarz inequality this is a PSD matrix.

Since $\phi(\textbf{x},\textbf{y})=\sum_{i}\phi(\textbf{x},y_{i})$, the pivoted version of $\nabla^{2}\phi(\textbf{x},\textbf{y})$ can be written as:

Where $\textbf{u}_{i}=\textbf{L}_{i}(\textbf{x}-\textbf{s}_{i})$. To solve this linear system we invoke the Sherman-Morrison-Woodbury formula.

If A, U, C, V are $n\times n$, $n\times k$, $k\times k$ and $k\times n$ matrices respectively, then:

In our case we have $\textbf{A}=\sum_{i}\alpha_{i}\textbf{L}_{i}$, $\textbf{C}=-\textbf{I}$, $\textbf{U}=[\sqrt{\beta_{1}}\textbf{u}_{1}\ldots\sqrt{\beta_{k}}\textbf{u}_{k}]$ and $\textbf{V}=\textbf{U}^{T}$. So the linear system that we need to evaluate becomes:

The system $\textbf{A}^{-1}\textbf{U}$ can be found using $k$ solves in $\textbf{A}=\sum_{i}\alpha_{i}\textbf{L}_{i}$, which is equivalent to the quadratic minimization problem. Multiplying this by U can be done in $O(k^{2}n)$ time and gives us $\textbf{U}^{T}\textbf{A}^{-1}\textbf{U}$. This $k\times k$ system can in turn be solved in $k^{\omega}$ time. The other terms can be applied to vectors in either solves in A or matrix multiples in U, taking $O(T(n,m)+kn)$ time. $\blacksquare$

Appendix D Other Variants

Although our formulation of $\mathcal{OBJ}$ as a sum of $L_{2}$ objectives differs syntactically from some common formulations, we show below that the more common formulation involving a $L_{2}$-squared fidelity term can be reduced to finding exact solutions to $\mathcal{OBJ}$ using 2 iterations of ternary search. Most other formulations differs from our formulation in the fidelity term, but more commonly have $L_{1}$ smoothness terms as well. Since the anisotropic smoothness term is a special case of the isotropic one, our discussion of the variations will assume anisotropic objectives.

The most common form of the total variation objective used in practice is one with $L_{2}^{2}$ fidelity term. This term can be written as $||\textbf{x}-\textbf{s}_{0}||_{2}^{2}$, which corresponds to the norm defined by $\textbf{I}=\textbf{L}_{0}$. This gives:

We can establish the value of $||\textbf{x}-\textbf{s}_{0}||_{\textbf{L}_{0}}^{2}$ separately by guessing it as a constraint. Since the $t^{2}$ is convex in $t$, the following optimization problem is convex in $t$ as well:

Also, due to the convexity of $t^{2}$, ternary searching on the minimizer of this plus $t^{2}$ would allow us to find the optimum solution by solving $O(\log{n})$ instances of the above problem. Taking square root of both sides of the $||\textbf{x}-\textbf{s}_{0}||_{\textbf{L}_{0}}^{2}\leq t^{2}$ condition and taking its Lagrangian relaxation gives:

Which by the min-max theorem is equivalent to:

The term being minimized is identical to our formulation and its objective is convex in $\lambda$ when $\lambda\geq 0$. Since $-\lambda t$ is linear, their sum is convex and another ternary search on $\lambda$ suffices to optimize the overall objective.

Another common objective function to minimize is where the fidelity term is also under $L_{1}$ norm. In this case the objective function becomes:

This can be rewritten as a special case of $\mathcal{OBJ}$ as:

Which gives, instead, a grouped least squares problem with $m+k$ groups.