Fast Newton methods for the group fused lasso

Matt Wytock
Machine Learning Dept.
Carnegie Mellon University
Pittsburgh, PA

Suvrit Sra
Machine Learning Dept.
Carnegie Mellon University
Pittsburgh, PA

Abstract
We present a new algorithmic approach to thegroup fused lasso, a convex model that approximates a multi-dimensional signal via an approximately piecewise-constant signal. Thismodel has found many applications in multiple change point detection, signal compression, and total variation denoising, though
existing algorithms typically using first-order
or alternating minimization schemes. In thispaper we instead develop a specialized projected Newton method, combined with a primal active set approach, which we show to
be substantially faster that existing methods.
Furthermore, we present two applications that
use this algorithm as a fast subroutine for amore complex outer loop: segmenting linear
regression models for time series data, andcolor image denoising. We show that on these
problems the proposed method performs very
well, solving the problems faster than stateof-the-art methods and to higher accuracy.

1

Introduction

Given a multivariate signal y1 , y2 , . . . , yT , with yt ∈Rn , the (weighted) group fused lasso (GFL) estimator(Bleakley and Vert, 2011; Alaı́z et al., 2013) attempts
to find a roughly “piecewise-constant” approximation
to this signal. It determines this approximation bysolving the optimization problem
minimize

x1 ,x2 ,...,xT

T
T
−1
X
1X
wi kxt −yt k22 +λt kxt −xt+1 k2 (1)
2 t=1
t=1

where x1 , x2 , . . . , xT are the optimization variables,w ∈ RT+ are weights for each time point, λ ∈ RT+−1
are regularization parameters, and k · k2 denotes the
Euclidean norm. Intuitively, the `2 norm on the difference between consecutive points encourages sparsity in

J. Zico Kolter
Machine Learning Dept.
Carnegie Mellon University
Pittsburgh, PA

this difference: each difference xt − xt+1 will typicallybe either full or identically zero at the solution, i.e., thesignal x will be approximately piecewise-constant. Thisapproach generalizes the 1D total variation norm (Tibshirani et al., 2005; Barbero and Sra, 2011), which considers only univariate signals. Owing to the piecewiseconstant nature of the approximate signals formed by
the group fused lasso, the approach has found applications in signal compression, multiple change-point
detection, and total variation denoising. Though several algorithms have been proposed to solve (1), tothe best of our knowledge these have involved, at theirfoundation, first-order methods such as projected gradient, block coordinate descent, or splitting methods.
Although such algorithms can sometimes obtain reasonable performance, they often fail to quickly find accurate solutions, especially when one wants to solve (1)to high precision as a “subroutine” (or prox-operator)
in a larger algorithm (Barbero and Sra, 2011).
In this paper, we develop a fast algorithm for solving
the optimization problem (1), based upon a projected
Newton approach. Our method can solve group fusedlasso problems to high numerical precision, often several
orders of magnitude faster than existing state-of-the-art
approaches. At its heart, our method involves dualizing the optimization problem (1) twice, in a particularmanner, to eliminate the non-differentiable `2 norm
and replace it by simple nonnegativity constraints; wesolve the reformulated problem to high accuracy via a
projected Newton approach. In order to fully exploit
the sparsity of large-scale instances, we combine theabove ideas with a primal active-set method that iteratively solves reduced-size problems to find the final setof non-zero differences for the original GFL problem.
Although our fast fused group lasso method is valuable
in its own right, its real power comes when used asa proximal subroutine in a more complex algorithm,
an operation that often needs to be solved thousandsof times. With this motivation in mind, we apply ourapproach to two applications: segmenting linear regression models, and color total variation image denoising.

We demonstrate the power of our approach in experiments with real and synthetic data, both for the basic
group fused lasso and these applications, and show
substantial improvement over the state of the art.

2

A fast Newton method for the GFL

We begin by adopting slightly more compact notation,and rewrite (1) (the primal problem) as
minimize
X

yT

(P)



i

where ai is the ith column of A; and D ∈ RT,T −1
denotes the first order differencing operator
1
0 0 ···
 −1 1 0 · · · 
D =  0 −1 1 · · · 
(4)


..
..
.. . .
.
.
.
.
so that XD takes the difference of the columns of X.
Dual problems

To solve (P), it is useful to look at its dual and (for ouralgorithm) a modified dual of this dual. To derive these
problems, we transform (P) slightly by introducing theconstraint V = XD, and corresponding dual variables
U ∈ Rn×T −1 . The Lagrangian is then given by
LP (X, U, V ) := 12 k(X − Y )W 1/2 k2F + kV Λk1,2+ tr U T (V − XD).

(5)

?

T

X = Y − UD W

−1?

, V = 0 iff kut k2 ≤ λt(6)where ut is the t-th column of U ; this leads to the dual
maximize − 12 kU DT W −1/2 k2F + tr U DT Y T
Usubject to kut k2 ≤ λt ,

t=1

(D)

t = 1, . . . , T − 1.

Indeed, several past algorithmic approaches have solved(D) directly using projected gradient methods, seee.g., (Alaı́z et al., 2013).
The basis of our algorithm is to form the dual of (D),but in a manner that leads to a different problem than(8)

where Z := diag(z), and leads to the dual problem (thedual of the dual of (P))
z≥0

T
−1
1
D + Z)−1 DT Y T + 12 (λ2 )T z,
2 Y D(D W(DD)
where λ2 denotes squaring λ elementwise. This procedure, taking the dual of the dual of the original
optimization problem, has transformed the original,
non-smooth problem into a smooth optimization problem subject to a non-negativity constraint, a setting
for which there are several efficient algorithms. Although (DD) is not easily solved via a standard form
semidefinite program—it involves a matrix fractionalterm, for which the standard semidefinite programming
form is computationally unattractive—it can be solved
efficiently by a number of methods for smooth, boundconstrained optimization. However, as we will see
below, the Hessian for this problem is typically poorly
conditioned, so the choice of algorithm for minimizing(DD) has a large impact in practice. Furthermore, because the z dual variables are non-zero only for thechange points of the original X variables, we expectthat for many regimes we will have very few non-zero z
values. These points motivate the use of projected Newton methods (Bertsekas, 1982), which perform Newton
updates on the variables not bound (z 6= 0).
2.2

A projected Newton method for (DD)

Denote the objective of (DD) as f (z); the gradient andHessian of f are given by∇z f (z) = − 21 (U 2 )T 1 + 12 λ2 ,Minimizing (5) analytically over X and V gives(7)

Minimizing over U analytically yields

min
; (2)

W := diag(w) and Λ := diag(λ); k · kF denotes the
Frobenius norm; k · k1,2 denotes the mixed `1,2 -norm
X
kAk1,2 :=kai k2 ,(3)

2.1

LD (U, z) = − 12 kU DT W −1/2 k2F + tr U DT Y T
XT −1+
zt (kut k22 − λ2t ).
U ? = Y D(DT W −1 D + Z)−1 ,1/2 2
1
kF + kXDΛk1,2
2 k(X − Y )W

where X, Y ∈ Rn×T denote the matricesX = x1 · · · xT , Y = y1 · · ·the original primal. In particular, noting that theconstraint kut k2 ≤ λt is equivalent to the constraintthat kut k22 ≤ λ2t , we can remove the non-differentiable`2 norm, and form the Lagrangian∇2z f (z) = U T U ◦ (DT W −1 D + Z)−1 ,(9)where as above U = Y D(DT W −1 D + Z)−1 , U 2 denotes elementwise squaring of U , and ◦ denotes theelementwise (Hadamard) product. The projected Newton method proceeds as follows: at each iteration, weconstruct the set of bound variables
I := {i : zi = 0 and (∇z f (z))i > 0}.

(10)

We then perform a Newton update only on those variables that are not bound (Ī, referred to as the freeset), and project back onto the feasible setzĪ ← zĪ − α(∇2z f (z))−1(∇z f (z))Ī + ,
(11)
Ī,Ī

Algorithm 1 Projected Newton for GFL
input signal Y ∈ Rn×T ; weights w ∈ RT+ ; regularization parameters λ ∈ RT+−1 ; tolerance 
output: optimized signal X ∈ Rn×T
initialization: z ← 0
repeat1. Form dual variables and gradientU ← Y D(DW−1−1

D + Z)∇z f (z) ← − 12 (U 2 )T 1 + 12 λ2
2. Compute active constraints
I ← {i : zi = 0 and (∇z f (z))i > 0}
3. Compute reduced Hessian and Newton direction
H ← UĪT UĪ ◦ (DT W −1 D + Z)−1Ī,Ī∆zĪ ← −H −1 (∇z f (z))Ī
4. Update variableszĪ ← [zĪ + α∆zĪ ]+where α is chosen by line search
until k(∇z f (z))Ī k2 ≤ 

where α is a step size (chosen by backtracking, interpolation, or other line search), and [·]+ denotes projection onto the non-negative orthant. The full method isshown in Algorithm 1. Although the projected Newton
method is conceptually simple, it involves inverting
several (possibly T × T matrices), which is impractical
if these were to be computed as general matrix operations. Fortunately, there is a great amount of structurethat can be exploited in this problem.
Efficiently solving Y D(DT W −1 D+Z)−1 . One key
operation for the weighted GFL problem is to solve
linear systems of the form DT W −1 D + Z, where W
and Z are diagonal. Fortunately, the first matrix ishighly structured: it is a symmetric tridiagonal matrix
11
w1 + w2
1
− w2

− w12
1
1
w2 + w3
− w13
..
.


···

··· 


DT W −1 D = 
,
0
··· 


..
..
.
.
(12)
and adding Z to it only affects the diagonal. LAPACK
has customized routines for solving problems of thisform: dpttrf (which computes the LDLT factorization of the matrix) and dptts2 (which computes thesolution to LDLT X = B via backsubstitution). For


0
− w13
1
1
w3 + w4..
.

our work, we modified this latter code slightly to solve
systems with the unknown on the left hand size, as isrequired for our setting; this lends a slight speedup by
exploiting the memory locality of column-based matrices. The methods factor T − 1 × T − 1 matrix in O(T )time, and solve n left hand sides in time O(T n).
Computing entries of (DT W −1 D+Z)−1 . The projected Newton method also requires more than justsolving equations of the form above: to compute theHessian, we must actually also compute entries of theinverse (DT W −1 D + Z)−1 — we need to compute theentries with rows and columns in Ī. Naively, this would
require solving k = |Ī| left hand sides, corresponding
to the unit bases for the entries in Ī; even using thefast solver above, this takes time O(T k). To speed
up this operation, we instead use a fast method forcomputing the actual entries of the inverse of this tridiagonal, using an approach based upon (Usmani, 1994);
this ultimately lets us compute the k 2 entries in O(k 2 )time, which can be much faster for small free sets.
Specifically, let a ∈ RT −1 and b ∈ RT −2 denotethe diagonal and the negative off-diagonal entries of
1
DT W −1 D + Z respectively (that is, ai = w1i + wi+1
+ zi
1
and bi = wi+1 ), we can compute individual entries of(DT W −1 D + Z)−1 as follows (the following adapts thealgorithm in (Usmani, 1994), but has enough simplifications for our case that we state it explicitly here).
Define θ, φ ∈ RT via the recursions
θi+1 = ai θi − b2i−1 θi−1 , i = 2, . . . , T − 1
θ1 = 1, θ2 = a1 ,φi = ai φi+1 − b2i φi+2 , i = T − 2, . . . , 1
φT = 1, φT −1 = aT −1 .

(13)Then, the (i, j) entry of (DT W −1 D + Z)−1 for j ≤ i
is given by
1 Yj−1 (DT W −1 D + Z)−1bi θi φj+1 . (14)
ij =
k=i
θT
Finally,
we can compute all the needed running prodQj−1
ucts k=i bi by computing a single cumulative sum
Piof the logs of the bi terms ci = j=1 log bi and thenQj−1using the equality k=i bi = exp(cj − ci ).
2.3

A primal active set approach

Using the two optimizations mentioned above, the projected Newton method can very quickly find a solution
accurate to numerical precision for medium sized problems (T and n on the order of thousands). However,for problems with substantially larger T , which are
precisely those we are most interested in for many GFL
applications, the approach above begins to break down.
There are two reasons for this: 1) The size of the freeset k = |I|, though often small at the final solution, canbe significantly larger at intermediate iterations; since
the Newton method ultimately does involve an O(k 3 )time to invert the Hessian restricted to the free set, thiscan quickly render the algorithm impractical. 2) Even
with small free sets, the basic O(T n) cost required fora single pass over the data at each Newton iteration
starts to dominate, especially since a significant number of iterations to find the correct free set may berequired (only after finding the correct free set does
one obtain quadratic convergence rates).

term in the objective is also equivalent, we have thatk(X − Y )W 1/2 k2F
=

k+1
X

1/2

k(x0i 1T − YJi0 )WJ 0 ,J 0 k2F
i

i

i=1

=

k+1
X

T

T

1/2

T
0
0
0
(wJ
YJi0 wJi0 + kYJi0 WJ 0 ,J 0 k2F0 1)x i xi − 2x
i



i

i=1

=

k+1
X

wi0 kx0i − yi0 k22 + c.

i=1

To overcome these problems, we consider a further layer
to the algorithm, which wraps our fast projected Newton solver inside a primal active-set method. The basic
intuition is that, at the optimal solution to the original
GFL problem, there will typically be very few change
points in the solution X ? (these correspond exactly
to those z variables that are non-zero). If we knew
these changes points ahead of time, we could solve a
substantially reduced (weighted) GFL problem thatwas equivalent to the original problem. Specifically, letJ ⊆ {1, . . . , T − 1} denote the optimal set of change
point locations for the primal problem. By the relationship of dual problems, this will be identical to theset of free variables Ī at the optimal solution, but since
we treat these differently in the algorithmic design weuse different notation. Then the original problem
minimize k(X − Y )W 1/2 k2F + kXDΛk1,2 ,X

where X ∈ R

n×T

(15)

, is equivalent to the reduced problem

minimizek(X 0 − Y 0 )W 0
0
X

1/2 2
kF + kX 0 D0 ΛJ ,J k1,2

(16)

with optimization variable X 0 ∈ Rn×k+1 for k = |J |,where D0 ∈ Rk+1×k denotes the same first order differences matrix but now over only k + 1-sized vectors,
and where Y 0 and W 0 = diag(w0 ) are defined by
wi0 =

X
j∈Ji0

wj ,

yi0 =

1 X
wj yj ,
wi00This equivalence motivates a primal active set method
where we iteratively guess the active set J (with somefixed limit on its allowable size), use the projected Newton algorithm to solve the reduced problem, and thenuse the updated solution to re-estimate the active set.
This is essentially equivalent to a common “block pivoting” strategy for non-negative least squares (Portugal
et al., 1994) or `1 methods (Lee et al., 2007), and has
been shown to be very efficient in practice (Kim andPark, 2010). The full algorithm, which we refer to as
Active Set Projected Newton (ASPN, pronounced “aspen”), is shown in Algorithm 2. In total, the algorithm
is extremely competitive compared to past approaches
to GFL, as we show in Section 4, often outperforming
the existing state of the art by orders of magnitude.

3

Applications

Although the ASPN algorithm for the group fusedlasso is a useful algorithm in its own right, part of
the appeal of a fast solver for this type of problem
is the possibility of using it as a “subroutine” withinsolvers for more complex problems. In this section we
derive such algorithms for two instances: segmentationof time-varying linear regression models and multichannel total variance image denoising. Both models
have been considered in the literature previously, andthe method presented here offers a way of solving these
optimization problems to a relatively high degree of
accuracy using simple methods.

(17)

j∈Ji

where we define Ji0 = {Ji−1 + 1, . . . , Ji } for i =1, . . . , k + 1 (i.e., Ji0 denotes the list of indices within
the ith segment, there being k+1 segments for k change
points). Furthermore, all these terms can be computedin time O(nk) via cumulative sums similar to the cumulative sum used for b above (which take O(T n) tocompute once, but which thereafter only require O(kn)to form the reduced problem).
To see this equivalence, note first that since X onlychanges at the points |J |, it immediately holds thatkXDΛk1,2 = kX 0 D0 ΛJ ,J k1,2 . To show that the other

3.1

Linear model segmentation

In this setting, we observe a sequence of input/output
pairs (at ∈ Rn , yt ∈ R) over time and the goal is
to find model parameters xt such that yt ≈ aTt xt (itis more common to denote the input itself as xt andmodel parameters θt , but the notation here is more in
keeping with the rest of this paper). Naturally, if xt
is allowed to vary arbitrarily, we can always find (aninfinite number of) xt ’s that fit the output perfectly,
but if we constrain the sum of norms kxt −xt−1 k2 , thenwe will instead look for piecewise constant segments inthe parameter space; this model was apparently first
proposed in Ohlsson et al. (2010).

Algorithm 2 Active Set Projected Newton (ASPN)
n×T

input signal Y ∈ R; weights w ∈ RT+ ; regularization parameters λ ∈ RT+−1 ; maximum active set
size kmax ; tolerance 
output: optimized signal X ∈ Rn×T
initialization: z ← 0
repeat1. Form dual variables and gradientalternating direction method of multipliers (ADMM),
an algorithm that has attracted great attention recently (e.g. Boyd et al. (2011)). Briefly, ADMM solves
problems of the form
minimize f (x) + g(z), subject to Ax + Bz = c, (20)x,zvia a sequence of alternating minimizations over x andz and dual variable updates.

U ← Y D(DW −1 D + Z)−11
1
∇z f (z) ← − (U 2 )T 1 + λ2
22
2. Compute active set, containing all non-zero zi ’s
and additional element with negative gradients,
up to size kmax

The “standard” ADMM algorithm. The simplest way to apply ADMM to (18), considered for thepure group fused lasso e.g., in Wahlberg et al. (2012),
is to introduce variables Z = XD, and formulate theproblem as
minimize kA vec X − Y k22 + kZΛk1,2
X,Z

0

J ← {i : zi > 0}J 1 ← {i : zi = 0, ∇z f (z) < 0}

After some derivations, this leads to the updates

1
J ← J 0 ∪ J1:k0
max −|J |X k+1 ← argmin kA vec X − yk22 + ρ2 kXD − Z k + U k k2FX

3. Form reduced problem (Y 0 , w0 ) for J using (17)and solve using projected Newton

Z k+1 ← argmin kZΛk1,2 + ρ2 kX k+1 D − Z + U k k2F
Z

U

k+1

k

← U + X k+1 D − Z k+1 ,(22)

zJ ← Projected-Newton(Y 0 , w0 , λJ )until k(∇z f (z))J k2 ≤ This model may be cast as the optimization problem
minimize kA vec X − yk22 + kXDΛk1,2 ,X(21)subject to XD = Z.

(18)

where X ∈ Rn×T is the same optimization variable
as previously, y ∈ RT denotes the vector of outputs,
vec denotes the vectorization of a matrix (stacking its
columns into a single column vector), and A ∈ RT ×T n
is the block diagonal matrix
 T

a1
0
0 ···
 0 aT2
0 ··· 


T
A= 0
.
(19)
0
a
··· 
3


..
..
.. . .
.
.
.
.
While this problem looks very similar to the ordinary
GFL setting, the introduction of the additional matrix A renders it substantially more complex. While
it is possible to adapt the Newton methods above tosolve the problem directly, much of the special problem
structure is lost, and it requires, for examples, forming
T n × T n block tridiagonal matrices, which is substantially more computationally intensive, especially forlarge n (the methods scale like O(n3 )). While optimization may still be possible with such approaches,
we instead adopt a different approach that builds on thewhere ρ acts effectively like a stepsize for the problem.
This set of updates is particularly appealing because
minimization over X and Z can both be computed inclosed form: the minimization over X is unconstrained
quadratic optimization, and has the solution
−1 T
AT A + ρF T F
A y + ρF T vec(Z k − U k ) (23)where F = (DT ⊗ I). Furthermore, these updates
can be computed very efficiently, since AT A + ρF T F
is block tridiagonal, and since this matrix does notchange at each iteration, we can precompute its (sparse)
Cholesky decomposition once, and use it for all iterations; using these optimizations, the X update takes
time O(T n2 ). Similarly, the Z update is a proximal
operator that can be solved by soft thresholding thecolumns of X k+1 D + U k (an O(T n) operation). Although these elements make the algorithm appealing,
they hide a subtle issue: the matrix (AT A + ρF T F X)
is poorly conditioned (owing to the poor conditioningof DT D), even for large ρ. Because of this, ADMM
needs a large number of iterations for converging toa reasonable solution; even if each iteration is quite
efficient, the overall algorithm can still be impractical.
ADMM using the GFL proximal operator. Alternatively, we can derive a different ADMM algorithm
by considering instead the formulation
minimize kA vec X − Y k22 + kZDΛk1,2
X,Z

subject to X = Z,

(24)

which leads to the iterative updates
X k+1 ← argmin kA vec X − yk22 + ρ2 kX − Z k + U k k2FX

Z k+1 ← argmin kZDΛk1,2 + ρ2 kX k+1 − Z + U k k2F
Z

U

k+1

k

← U + X k+1 − Z k+1 .
(25)

The X update can still be computed in closed form
−1 T
vec X k+1 = AT A + ρI
A y + ρ vec(Z k − U k ) ,which is even simpler to compute than in the previous
case, since AT A+ρI is block diagonal with blocks ai aTi +ρI, which can be solved for in O(n) time; thus the entire
X update takes times O(T n). The downside is that
the Z update, of course, can no longer be solved with
soft-thresholding. But the Z update here is precisely
in the form of the group fused lasso; thus, we can use
ASPN directly to perform the Z update. The main
advantage here is that the matrix AT A + ρI is muchbetter conditioned, which translates into many fewer
iterations of ADMM. Indeed, as we show below, thisapproach can be many orders of magnitude faster thanstraight ADMM, which is already a very competitivealgorithm for solving these problems.
3.2

Color total variation denoising

Unfortunately, this optimization problem cannot be
solved directly via the group fused lasso, as the difference penalties on the rows and columns for the same
matrix X render the problem quite different from thebasic GFL. We can, however, adopt an approach similar to the one above, and create separate variables
corresponding to the row and column slices, plus a
constraint that they be equal; formally, we solve
minimize
X,Z

+ λ

m−1
n
XX

m n−1
X
X

kX:,i,j − X:,i,j+1 k2

i=1 j=1

kX:,i,j − X:,i+1,j k2 ,kX:,i,: − Y:,i,: k2F + λkX:,i,: Dk1,2



i=1

+

n
X

kZ:,:,j − Y:,:,j k2F + λkZ:,:,j Dk1,2



j=1

subject to X = Z.
(27)
The major advantage of this approach is that it decomposes the problem into m + n independent GFL tasks,
plus a meta-algorithm that adjusts each sub-problem
to make the rows and columns agree. Several suchalgorithms are possible, including ADMM; we present
here a slightly simpler scheme known as the “proximal Dykstra” method (Combettes and Pesquet, 2011),
which has been previously applied to the case of (single
channel, i.e., black and white) total variation denoising (Barbero and Sra, 2011). Starting with X 0 = Y ,P 0 = 0, Q0 = 0, the algorithm iterates as follows:

Next, we consider color total variation denoising example. Given an m × n RGB image represented as a thirdorder tensor, Y ∈ R3×m×n , total variation image denoising (Rudin et al., 1992; Blomgren and Chan, 1998)
attempts to find an approximation X ∈ R3×m×n suchthat differences between pixels in X favor being zero.
It does this by solving the optimization problem
minimize 12 kX − Y k2F + λ
X

m
X

k+1
k
k
Z:,:,j← GFL(X:,:,j+ P:,:,j, λ), j = 1, . . . , n

P k+1 ← P k + X k − Z k+1
k+1
k+1
X:,i,:← GFL(Z:,i,:+ Qk:,i,: , λ), i = 1, . . . , m(28)

Qk+1 ← Qk + Z k+1 − X k+1 .
Typically, very few iterations (on the order of 10) of
this outer loop are need to converge to high accuracy.
Furthermore, because each of the m or n GFL problems
solved in the first and third steps are independent, they
can be trivially parallelized.

i=1 j=1(26)
corresponding to an `2 norm penalty on the difference
between all adjacent pixels, where each pixel X:,i,j isrepresented as a 3 dimensional vector. We can write
this as a sum of m + n group fused lasso problems
minimizeXm
X

kX:,i,: − Y:,i,: k2F + λkX:,i,: Dk1,2



i=1

+

n
X


kX:,:,j − Y:,:,j k2F + λkX:,:,j Dk1,2 ,j=1

where X:,i,: ∈ R3×n denotes the slice of a single rowof the image and X:,:,j ∈ R3×n denotes the slice of asingle column.

4

Experimental results

We present experimental results for our approaches,both on the basic group fused lasso problem, where wecompare to several other potential approaches, and on
the two applications of linear model segmentation andcolor total variation denoising. C++ and MATLAB
code implementing our methods is available at http://www.cs.cmu.edu/~mwytock/gfl/.
4.1

Group fused lasso

Here we evaluate the ASPN algorithm versus several
alternatives to solving the group fused lasso problem,
evaluated on both synthetic and real data. Figure 1

3
10

ASPN
GFL
PG
APG
DR
LBFGS−Bf − f*

1
10
−1
10
−3
10
−5
10
0

1

2
3
Time (seconds)

4

5

5
ASPN
GFL
PG
APG
DR
LBFGS−Bf − f*

4
10

0

1
10
−2
10
−5
10

−5

0

5

10

1000 2000 3000 4000 5000 6000 7000 8000 9000 10000shows a synthetic time series with T = 10, 000, n = 100,
and 10 discrete change points in the data; the data was
generated by uniformly sampling the change points,
sampling the mean of each segment from N (0, I), andthen additional Gaussian noise. Figure 1 shows the recovered signal using the group fused lasso with wt = 1,
λt = 20. In Figure 2, we show timing results for thisproblem as well as a smaller problem with T = 1000
and n = 10; we compare ASPN to GFLseg (Bleakley and Vert, 2011) (which uses coordinate descenton the primal problem), an accelerated projected gradient on the dual (i.e., the FISTA algorithm) (Beckand Teboulle, 2009), Douglas-Rachford splitting (Combettes and Pesquet, 2007) (a generalization of ADMM
that performs slightly better here), a projected gradienton the dual (Alaı́z et al., 2013), and LBFGS-B (Byrdet al., 1995) applied to the dual of the dual. In allcases, ASPN performs as well as (often much betterthan) the alternatives.
Next, we evaluate how the ASPN algorithm scales as afunction of the number of time points T and the numberof change points at the solution, k. In Figure 3, thefirst set of experiments shows that when the numberof change points at the solution is fixed (k = 10), theamount of time required for a highly accurate solution
remains small even for large T , agreeing with analysis
that shows the number of operations required is O(T ).
In particular, a solution accurate to 10−6 is found in
4.8 seconds on a problem with T = 106 time points.

25

30

Figure 2: Above: timing results on synthetic problem
with T = 1000, n = 10. Below: timing results on
synthetic problem with T = 10000, n = 100.

1
0.5
X

Figure 1: Above: synthetic change point data, with
T = 10000, n = 100, and 10 true change points. Below:
recovered signal.

15
20Time (seconds)

0
−0.5
−1
0.5

1

1.5
T

2

2.5

3
4
x 10

Figure 4: Above: Lung data from (Bleakley and Vert,2011). Below: recovered signal using group fused lasso.

However, in the next set of experiments, we see thatcompute time grows rapidly as a function of k due tothe O(k 3 ) operations required to compute the Newton
step, suggesting that the proposed method is mostappropriate for problems with sparse solutions.

ActualO(T)110ActualO(k3)210Time (seconds)

TIme (seconds)

2
10

0
10
−1
10
−2
10

0
10

−2
10

2
10

3
10

4
10
Number of time points (T)

5
10

6
10

0
10

1
10

2
10
Number of change points (k)

3
10

4
10

Figure 3: Left: timing results vs. number of change points at solution for synthetic problem with T = 10000 and
n = 10. Right: timing results for varying T , n = 10, and sparse solution with 10 change points.

0.5

4
10

f − f*

1
10

y

ASPN
GFL
PG
APG
DR
LBFGS−B

0

−2
10

−5
10

−0.5
0

5

10

15
20
Time (seconds)

25

0

30

2000

4000

6000

8000

10000

6000

8000

10000

T

6
10

0.2

ASPN
GFL
PG
APG
DR
LBFGS−B

0.1
Xf − f*

3
10

0
10

0
−0.1
−0.2

−3
10
0

50

100
Time (seconds)

150

200

0

2000

4000
T

Figure 5: Above: Timing results on bladder problem,
T = 2143, n = 57. Below: Timing results on lungproblem, T = 31708, n = 18.

Figure 6: Above: Observed autoregressive signal zt .
Below: true autoregressive model parameters.

Finally, we evaluate the algorithm on two real timeseries previously used with the group fused lasso (Bleakley and Vert, 2011), from DNA profiles of bladder andlung cancer sequences. Figure 4 shows one of these
two series, along with the approximation produced by
the group fused lasso. Figure 5 shows timing results
for the above methods again on this problem: here weobserve the same overall behavior, that ASPN typicallydominates the other approaches.

regressive time series models. In particular, we observe
some time series z1 , . . . , zT , and we fit a linear model tothis data zt ≈ aTt xt where at = (zt−1 , zt−2 , . . . , zt−n ).
Figure 6 shows an example time series generated bythis process, as well as the true underlying model thatgenerated the data (with additional noise). This is
the rough setting used in (Ohlsson et al., 2010), which
was the first example we are aware of that uses such
regularization techniques within an linear regression
framework. Figure 7 shows the model parameters recovered using the method from Section 3.1, which here
match the ground truth closely.

4.2

Linear regression segmentation

Here we apply the two different ADMM methods discussed in Section 3.1 to the task of segmenting auto-

Of more importance, though, is the comparison between the two different ADMM approaches. Figure 8

0.2

X

0.1
0
−0.1
−0.2
0

2000

4000

6000

8000

10000

T

Figure 9: Left: original image. Middle: image corrupted with Gaussian noise. Right: imaged recovered
with total variation using proximal Dykstra and ASPN.
0
10

0.2

ADMM
Dykstra+ASPN
−1
10

0

(f−f*)/f*

X

0.1

−0.1

−2
10

−0.2
0

2000

4000

6000

8000

10000

−3
10

T

0

Figure 7: Above: Autoregressive parameters recovered
with “simple” ADMM algorithm. Below: parameters
recovered using alternative ADMM w/ ASPN.

ADMM w/ TV
ADMM

f − f*
−5
10

0

10

20

30

40

50

60

Time (seconds)

Figure 8: Convergence of simple ADMM versus alternative ADMM w/ ASPN.

shows convergence versus running time and here the“simple” ADMM approach, which encodes the difference operator in the constraints (and thus has simpler
updates), converges significantly slower than our alternative. Importantly, the X axis in this figure ismeasured in time, and we emphasize that even thoughthe “simple” ADMM updates are individually slightly
faster (they do not involve GFL subproblems), theiroverall performance is much poorer. Further, as illustrated in Figure 7, the “simple” ADMM approachnever actually obtains a piecewise constant X except
at the optimum, which is never reached in practice.

40

60
Iteration

80

100

120

Figure 10: Comparison of proximal Dykstra method
to ADMM for TV denoising of color image.
4.3

0
10

20

Color total variation denoising

Finally, as described in Section 3.2, we apply the proximal Dykstra algorithm, using ASPN as a fast subroutine, to color image denoising. Figure 9 shows a
256x256 image generated by combining various solidcolored shapes, corrupted with per-RGB-component
noise of N (0, 0.1), and then recovered with total variation denoising. There has been enormous work on
total variation denoising, and while a full comparison
is beyond the scope of this paper, ADMM or methods such as those used by the FTVd routines in Yanget al. (2009), for instance, are considered to be someof the fastest for this problem. In Figure 10, we show
the performance of our approach and ADMM versus
iteration number, and as expected observe better convergence; for single-core systems, ADMM is ultimately
a better solution for this problem, since each iterationof ADMM takes about 0.767 seconds in our implementation whereas 512 calls to ASPN take 20.4 seconds.
However, the advantage to the ASPN approach is that
all these calls can be trivially parallelized for a 256X
speedup (the calls are independent and all code is CPUbound), whereas parallelizing a generic sparse matrix
solve, as needed for ADMM-based approaches, is much
more challenging and thus per-iteration performancehighlights the potential benefits of the ASPN approach.
Acknowledgements This work was supported in partby the National Science Foundation under award IIS1320402, and by support from Google.

References
Alaı́z, C. M., Jiménez, Á. B., and Dorronsoro, J. R.
(2013). Group fused lasso. In International Conference on Artificial Neural Networks and Machine
Learning (ICANN), pages 66–73.
Barbero, Á. and Sra, S. (2011). Fast newton-type methods for total variation regularization. In Proceedingsof the 28th International Conference on Machine
Learning (ICML), pages 313–320.
Beck, A. and Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences,2(1):183–202.
Bertsekas, D. P. (1982). Projected newton methods for optimization problems with simple constraints. SIAM Journal on control and Optimization,20(2):221–246.
Bleakley, K. and Vert, J.-P. (2011). The group fusedlasso for multiple change-point detection. arXiv
preprint arXiv:1106.4199.
Blomgren, P. and Chan, T. F. (1998). Color TV: total
variation methods for restoration of vector-valued
images. Image Processing, IEEE Transactions on,
7(3):304–309.
Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2011). Distributed optimization and statistical learning via the alternating direction method ofmultipliers. Foundations and Trends R in Machine
Learning, 3(1):1–122.
Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C. (1995). A
limited memory algorithm for bound constrained optimization. SIAM Journal on Scientific Computing,16(5):1190–1208.
Combettes, P. L. and Pesquet, J.-C. (2007). A DouglasRachford splitting approach to nonsmooth convex
variational signal recovery. Selected Topics in Signal
Processing, IEEE Journal of, 1(4):564–574.
Combettes, P. L. and Pesquet, J.-C. (2011). Proximal
splitting methods in signal processing. In Fixedpoint algorithms for inverse problems in science andengineering, pages 185–212. Springer.
Kim, J. and Park, H. (2010). Fast active-set-type
algorithms for L1-regularized linear regression. In
International Conference on Artificial Intelligenceand Statistics, pages 397–404.
Lee, H., Battle, A., Raina, R., and Ng, A. Y. (2007).
Efficient sparse coding algorithms. In Neural Information Processing Systems.
Ohlsson, H., Ljung, L., and Boyd, S. (2010). Segmentation of ARX-models using sum-of-norms regularization. Automatica, 46(6):1107–1111.

Portugal, L. F., Judice, J. J., and Vicente, L. N. (1994).
A comparison of block pivoting and interior-point
algorithms for linear least squares problems with
nonnegative variables. Mathematics of Computation,
63(208):625–643.
Rudin, L. I., Osher, S., and Fatemi, E. (1992). Nonlinear total variation based noise removal algorithms.
Physica D: Nonlinear Phenomena, 60(1):259–268.
Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., andKnight, K. (2005). Sparsity and smoothness via thefused lasso. Journal of the Royal Statistical Society:
Series B (Statistical Methodology), 67(1):91–108.
Usmani, R. A. (1994). Inversion of a tridiagonal jacobi
matrix. Linear Algebra and Its Applications, 212:413–
414.
Wahlberg, B., Boyd, S., Annergren, M., and Wang,
Y. (2012). An ADMM algorithm for a class of total
variation regularized estimation problems. arXiv
preprint arXiv:1203.1828.
Yang, J., Yin, W., Zhang, Y., and Wang, Y. (2009). A
fast algorithm for edge-preserving variational multichannel image restoration. SIAM Journal on Imaging
Sciences, 2(2):569–592.