Path of Solutions for Fused Lasso Problems

Torpong Nitayanont

, Cheng Lu

and Dorit S. Hochbaum

Department of Industrial Engineering and Operations Research, University of California, Berkeley, Berkeley, CA, U.S.A.

Keywords:

Fused Lasso, Path of Solutions, Minimum Cut, Hyperparameter Selection, Signal Processing.

Abstract:

In a fused lasso problem on sequential data, the objective consists of two competing terms: the ﬁdelity term

and the regularization term. The two terms are often balanced with a tradeoff parameter, the value of which

affects the solution, yet the extent of the effect is not a priori known. To address this, there is an interest in

generating the path of solutions which maps values of this parameter to a solution. Even though there are

inﬁnite values of the parameter, we show that for the fused lasso problem with convex piecewise linear ﬁdelity

functions, the number of different solutions is bounded by n

q where n is the number of variables and q is the

number of breakpoints in the ﬁdelity functions. Our path of solutions algorithm, PoS, is based on an efﬁcient

minimum cut technique. We compare our PoS algorithm with a state-of-the-art solver, Gurobi, on synthetic

data. The results show that PoS generates all solutions whereas Gurobi identiﬁes less than 22% of the number

of solutions, on comparable running time. Even allowing for hundreds of times factor increase in time limit,

compared with PoS, Gurobi still cannot generate all the solutions.

1 INTRODUCTION

We consider here the class of convex piecewise linear

fused lasso problems:

(PFL) min

,...,x

∑

i=1

) + λ

n−1

∑

i=1

− x

i+1

(1)

where each ﬁdelity function f

) is a convex piece-

wise linear function and λ is a parameter that controls

the tradeoff between the ﬁdelity term and the regular-

ization term, which penalizes the differences between

consecutive variables.

Applications of problems in the class of PFL in-

clude the study of DNA copy number gains and losses

on chromosomes using Comparative Genomic Hy-

bridization (CGH) (Eilers and De Menezes, 2005). In

CGH, each piecewise linear function takes the form

of τ(x

− a

)

+ (1 −τ)(a

− x

)

for a given τ ∈ [0,1]

and data point a

, and the problem is called quantile

regression. In signal processing (Storath et al., 2016),

each piecewise linear function is f

) = w

− a

for a given weight w

and data point a

. Each piece-

wise linear function in these examples contains ex-

actly one breakpoint.

https://orcid.org/0009-0002-6976-1951

https://orcid.org/0000-0001-5137-7199

https://orcid.org/0000-0002-2498-0512

Clearly, the optimal solution of PFL (1) varies

with regard to the tradeoff parameter λ. In the works

mentioned above, the choice of λ is either handpicked

(Storath et al., 2016) or selected via cross validation

(Eilers and De Menezes, 2005). However, there could

still be λ values outside of the search range that lead

to desirable results, which go unnoticed, or the reso-

lution of the cross validation grid search could still be

made ﬁner. Our goal is to ﬁnd the path of solutions to

PFL for all λ ∈ [0,∞), that is we want to identify a set

of ranges or intervals of λ (we call them λ-ranges) that

span [0, ∞), i.e. Λ = {[0,λ

),[λ

,λ

),...,[λ

,∞)},

such that the solutions to PFL for two values of λ from

the same range are similar.

In addition to the identiﬁcation of the set of λ-

ranges, we also ﬁnd the optimal solution for each

λ-range, {x

∗

([0,λ

)),x

∗

([λ

,λ

)),...,x

∗

([λ

,∞))}.

Note that we use the notation of x

∗

(λ) to refer to the

solution of PFL for a speciﬁc value of λ. The notation

of x

∗

([λ

,λ

)), for example, refers to the optimal so-

lution to PFL for any value of λ in the range [λ

,λ

It also implies that the solution for all values of λ in

this range are identical.

The relevant literature on ﬁnding the path of solu-

tions include the work on fused-lasso signal approx-

imator (FLSA) by Hoeﬂing (2010). The problem is

Nitayanont, T., Lu, C. and Hochbaum, D.

Path of Solutions for Fused Lasso Problems.

DOI: 10.5220/0012433200003654

Paper published under CC license (CC BY-NC-ND 4.0)

In Proceedings of the 13th International Conference on Pattern Recognition Applications and Methods (ICPRAM 2024), pages 107-118

ISBN: 978-989-758-684-2; ISSN: 2184-4313

107

deﬁned as follows:

min

,...,x

∑

i=1

− a

)

+ λ

∑

i=1

+ λ

n−1

∑

i=1

− x

i+1

Hoeﬂing (2010) solves for the optimal solution

∗

(λ

,λ

) for all λ

and λ

by ﬁrst ﬁnding the

solution for λ

= 0 in terms of λ

, x

∗

(0,λ

), in

O(nlogn). Then, x

∗

(λ

,λ

) is computed as a func-

tion of x

∗

(0,λ

). Tibshirani (2011) presents path of

solution algorithms to a generalized lasso problem,

which is formulated as

min

x∈R

∥b − Ax∥

+ λ∥Dx∥

for different conﬁgurations of A and D such as when

D is an arbitrary matrix and A = I. Note that the reg-

ularization term of this problem can be made equiva-

lent to our regularization term,

∑

n−1

i=1

− x

i+1

|, when

D is set to an appropriate matrix.

Another work on the path of solutions is the work

by Wang et al. (2006), which solves the regularized

least absolute deviation regression problem (RLAD):

min

x∈R

||y − Ax||

+ λ||x||

given A ∈ R

n×p

,y ∈ R

for every possible value of λ

in O(log (np) min(p,n)

). They applied their method

on the image reconstruction problem.

Generalized isotonic median regression or GIMR

problem, which has the following formulation:

(GIMR) min

,...,x

∑

i=1

) +

n−1

∑

i=1

i,i+1

· (x

− x

i+1

)

n−1

∑

i=1

i+1,i

· (x

i+1

− x

)

(2)

for a given set of {d

i,i+1

i+1,i

}

n−1

i=1

, is solved by an

algorithm, called HL-algorithm hereafter, in the work

of Hochbaum and Lu (2017). HL-algorithm solves

GIMR by formulating the problem into a sequence

of minimum cut problems, which can be solved ef-

ﬁciently in O(q log n) where q is the total number of

breakpoints in ﬁdelity functions. This is the best com-

plexity for solving GIMR to date. PFL is a special

case of GIMR and therefore can be solved, for a ﬁxed

λ, using the HL-algorithm. The path of solutions al-

gorithm proposed here that solves PFL for all λ ≥ 0 is

an extension of the HL algorithm. We introduce here

the notations and concepts that will be used through-

out the paper.

Notation and Preliminaries. GIMR (2) can be

viewed as deﬁned on a bi-directional path graph G =

(V,A) with node set V = {1,2,...,n} and arc set

A = {(i,i + 1),(i + 1,i)}

i=1,...,n−1

. Each node i ∈ V

corresponds to the variable x

. Let the node interval

[i, j] in the graph G for i ≤ j be the subset of consec-

utive nodes in V, {i, i + 1, . .., j − 1, j}.

We deﬁne an associated graph G

with the set of

vertices V

= V ∪ {s,t} and the set of arcs A

= A ∪

∪A

. The two appended nodes s and t are called the

source and the sink node, respectively. A

= {(s,i) :

i ∈ V } and A

= {(i,t) : i ∈ V } are the sets of source

adjacent arcs and sink adjacent arcs. Each arc (i, j) ∈

has an associated nonnegative capacity c

i, j

An s,t-cut is a partition of V

,({s} ∪ S,T ∪ {t}),

where T =

S = V \S. For simplicity, we refer to an s,t-

cut partition as (S,T ). We refer to S as the source set

of the cut, excluding s, and T the sink set, excluding

t. For each node i ∈ V , we deﬁne its status in graph

as status(i) = s if i ∈ S (referred as an s-node),

otherwise status(i) = t)(i ∈ T ) (referred as a t-node).

The capacity of a cut (S,T ) is deﬁned as C({s} ∪

S,T ∪ {t}) where C (V

) =

∑

i∈V

, j∈V

i, j

. A mini-

mum s,t-cut in G

is an s,t-cut (S,T ) that minimizes

C({s} ∪S, T ∪{t}). Hereafter, any reference to a min-

imum cut is to the unique minimum s,t-cut with the

maximal source set, that is, the source set that is not

contained in any other source set of a minimum cut in

case that there are multiple minimum cuts.

A convex piecewise linear function f

) is

speciﬁed by its ascending list of q

breakpoints, a

i,1

i,2

< . . . < a

i,q

, and the slopes of the q

+ 1 linear

pieces between every two adjacent breakpoints, de-

noted by w

i,0

< w

i,1

< . . . < w

i,q

. To deﬁne a break-

point more formally, each convex piecewise linear

function f

) can be viewed as a maximum of mul-

tiple afﬁne functions, which are sorted according to

their slopes. We call the x-coordinate of an intersec-

tion of two consecutive afﬁne functions a breakpoint.

Let the sorted list of the union of q breakpoints of

all the n convex piecewise linear functions be a

, j

< . . . < a

, j

(w.l.o.g. we may assume that the

n sets of breakpoints are disjoint (Hochbaum and Lu,

2017)), where a

, j

, the k-th breakpoint in the sorted

list, is the breakpoint between the ( j

− 1)-th and the

-th linear pieces of f





Overview of Paper. Section 2 describes the algo-

rithm that solves PFL (1) for a ﬁxed value of λ, the

HL-algorithm. In Section 3, we give a description of

the path of solutions algorithm, or PoS, for ℓ

-ﬁdelity

fused lasso problem. In Section 4, we generalize it

to PFL with convex piecewise linear ﬁdelity function.

We provide the bound of the number of different so-

lutions of PFL across all nonnegative λ values as well

as the time complexity of the algorithm in Section 5.

Lastly, we conclude with experimental results in Sec-

tion 6.

ICPRAM 2024 - 13th International Conference on Pattern Recognition Applications and Methods

108

Figure 1: G

(α): For node i, if its right sub-gradient at

α, ( f

)

′

(α), is negative, s is connected to i with an arc of

capacity c

s,i

= −( f

)

′

(α). Otherwise, i is connected to t

with an arc of capacity c

i,t

= ( f

)

′

(α).

2 ALGORITHM THAT SOLVES

FUSED LASSO FOR A FIXED

VALUE OF λ

The PFL problem (1) that we consider in this work

can be solved for a ﬁxed λ using the HL-algorithm,

introduced by Hochbaum and Lu (2017). The HL-

algorithm solves problems in the class of generalized

isotonic median regression or GIMR (2). PFL is a

special case of GIMR (2) when d

i,i+1

= d

i+1,i

= λ

for all i ∈ {1,2,...,n − 1}. In this section, we give

a description of how the HL-algorithm solves GIMR,

which also explains for the case of PFL.

We construct a parametric graph G

(α) =

), shown in Figure 1, that is associated with

the bi-directional path graph G = (V,A), for any scalar

value α. The capacities of arcs (i, i + 1),(i + 1,i) ∈ A

are c

i,i+1

= d

i,i+1

and c

i+1,i

= d

i+1,i

respectively. Arc

(s,i) ∈ A

has capacity c

s,i

= max{0,−( f

)

′

(α)} and

arc (i,t) ∈ A

has capacity c

i,t

= max{0, ( f

)

′

(α)},

where ( f

)

′

(α) is the right sub-gradient of function

(·) at argument α. (One can select instead the left

sub-gradient.) For any given value of α, we have ei-

ther c

s,i

= 0 or c

i,t

= 0, that is, i is never connected to

both s and t for the same α.

The link between the minimum cut for any given

value of α and the optimal solution to GIMR (2) is

characterized in the following threshold theorem:

Theorem 2.1 (Threshold theorem (Hochbaum,

2001)). For any given α, let S

∗

be the maximal source

set of the minimum cut in graph G

(α). Then, there is

an optimal solution x

∗

to GIMR (2) satisfying x

∗

≥ α

if i ∈ S

∗

and x

∗

< α if i ∈ T

∗

An important property of G

(α) is that the capac-

ities of source adjacent arcs and sink adjacent arcs are

nonincreasing and nondecreasing functions of α, re-

spectively. The capacities of all the other arcs are con-

stants. This implies the following nested cut property:

Lemma 2.2 (Nested cut property (Gallo et al., 1989;

Hochbaum, 2001, 2008)). For any two parameter val-

ues α

≤ α

, let S

and S

be the respective max-

imal source set of the minimum cuts of G

(α

) and

(α

), then S

⊇ S

We remark that the above threshold theorem and

nested cut property both work not only for GIMR

(2) deﬁned on a bi-path graph, but also for a gen-

eralization of GIMR that is deﬁned on arbitrary (di-

rected) graphs, where the regularization term may pe-

nalize the difference between arbitrary pairs of vari-

ables rather than just consecutive variables.

Based on the threshold theorem, it is sufﬁcient

to solve the minimum cuts in the parametric graph

(α) for all values of α, in order to solve GIMR (2).

In piecewise linear functions, the right sub-gradients

for α values between any two adjacent breakpoints

are constant. Thus, the source and sink adjacent arc

capacities remain constant for α between any two ad-

jacent breakpoint values in the sorted list of break-

points over all the n convex piecewise linear func-

tions. This result leads to the following lemma given

by Hochbaum and Lu (2017).

Lemma 2.3. The minimum cuts in G

(α) remain

unchanged for α assuming any value between any

two adjacent breakpoints in the sorted list of break-

points of all the n convex piecewise linear functions,

{ f

)}

i=1,...,n

Thus, the values of α to be considered can be re-

stricted to the set of breakpoints of the n convex piece-

wise linear functions, { f

)}

i=1,...,n

. The HL algo-

rithm solves GIMR (2) by efﬁciently computing the

minimum cuts of G

(α) for subsequent values of α

in the ascending list of breakpoints, a

, j

< a

, j

... < a

, j

Let G

, for k ≥ 1, denote the parametric graph

(α) for α equal to a

, j

, i.e., G

= G

, j

). For

k = 0, we let G

= G

, j

− ε) for a small value

of ε > 0. Let (S

) be the minimum cut in G

for k ≥ 0. Recall that S

is the maximal source set.

The nested cut property (Lemma 2.2) implies that

⊇ S

k+1

for k ≥ 0. Based on the threshold theorem

and the nested cut property, we know that for each

node j ∈ {1, . ..,n}, x

∗

= a

, j

for the index k such

that j ∈ S

k−1

and j ∈ T

The HL-algorithm generates the respective mini-

mum cuts of G

in increasing order of k. It is shown

by Hochbaum and Lu (2017) that (S

) can be com-

puted from (S

k−1

) in time O(log n). Hence the

total complexity of the algorithm is O(q log n). The

efﬁciency of updating (S

) from (S

k−1

) is

based on the following key results.

The update of the graph from G

k−1

to G

is simple

as it only involves a change in the capacities of the

source and sink adjacent arcs of i

,(s,i

) and (i

,t).

Recall that from G

k−1

to G

, the right sub-gradient

Path of Solutions for Fused Lasso Problems

109

( f

)

′

changes from w

, j

−1

, the slope of the k-th lin-

ear piece of f

, to w

, j

, the slope of the (k + 1)-th

linear piece of f

. Thus, the changes of c

s,i

and c

from G

k−1

to G

depend on the signs of w

, j

−1

and

, j

, which have three possible cases.

Case 1. w

, j

−1

≤ 0, w

, j

≤ 0 : c

s,i

changes from

−w

, j

−1

to −w

, j

while c

remains zero.

Case 2. w

, j

−1

≤ 0, w

, j

≥ 0 : c

s,i

changes from

−w

, j

−1

to 0 and c

changes from 0 to w

, j

Case 3. w

, j

−1

≥ 0,w

, j

≥ 0 : c

changes from

, j

−1

to w

, j

while c

s,i

remains zero.

Note that the update from G

k−1

to G

does not

involve the values of d

i,i+1

and d

i+1,i

in GIMR (2), of

which both take the values of λ in PFL (1).

Based on the nested cut property, for any node

i, if i ∈ T

k−1

, then i remains in T

and the sink set

for all subsequent cuts. Hence, an update of the

minimum cut in G

from the minimum cut in G

k−1

can only involve shifting some nodes from source set

k−1

to sink set T

. Formally, the relation between

k−1

) and (S

) is characterized in Lemma

2.4 and 2.5 given by Hochbaum and Lu (2017):

Lemma 2.4. If i

∈ T

k−1

, then (S

) = (S

k−1

We deﬁne s-interval in graph G to be the interval

of consecutive s-nodes, that is the set of consecutive

nodes that are in the source set S of the minimum cut

of G. An s-interval containing i refers to an interval of

s-nodes in G that contains node i. If i is a t-node, then

such s-interval is an empty interval. The maximal s-

interval containing i refers to the s-interval containing

i that is not a subset of any other such s-intervals.

Lemma 2.5. If i

∈ S

k−1

, then all the nodes that

change their status from s in G

k−1

to t in G

must

form a (possibly empty) s-interval of i

in G

k−1

Both lemmas imply that the derivation of the min-

imum cut in G

from the minimum cut in G

k−1

can

be done by ﬁnding an optimal s-interval of nodes

containing i

that change their status from s to t.

If this interval of nodes, referred to as node sta-

tus change interval, is [i

] then it follows that

) = (S

k−1

\[i

],T

k−1

∪ [i

]). Combining

this with the threshold theorem (Theorem 2.1), all x

for i ∈ [i

] have their optimal solutions in GIMR

(2) equal to a

, j

To ﬁnd the node status change interval from G

k−1

to G

that gives the minimum cut in G

, we rely on

the following lemma from Hochbaum and Lu (2017).

Lemma 2.6. Given the maximal s-interval contain-

ing i

in G

k−1

, [i

], the node status change inter-

val [i

] is the optimal solution to the following

problem deﬁned on G

min

]

C({s},[i

]) +C([i

]\[i

],{t}) (3)

+ C([i

]\[i

],[i

k2]

] ∪ {i

− 1,i

+ 1})

s.t. i

∈[i

] ⊆ [i

] or [i

] =

Solving problem (3) in Lemma 2.6 can be easier

than directly solving the minimum cut problem on

. Our path of solutions algorithm rely on all of

the stated theorems and lemmas, but particularly on

Lemma 2.5 and 2.6.

3 PATH OF SOLUTIONS

ALGORITHM FOR ℓ

FIDELITY

FUNCTIONS

The proposed path of solutions algorithm, or PoS,

ﬁnds the optimal solution x

∗

for all nonnegative λ.

We ﬁrst demonstrate in this section how PoS solves a

speciﬁc PFL problem, the ℓ

-ﬁdelity fused lasso:

(ℓ

-PFL) min

,...,x

∑

i=1

− a

| + λ

n−1

∑

i=1

− x

i+1

| (4)

for all nonnegative λ, given a set of breakpoints

}

i=1

and a set of positive penalty weights {w

}

i=1

As we consider ℓ

-PFL, which is a special case

of both GIMR and PFL, we substitute some notations

introduced in Section 2 by simpler terms.

For ℓ

-PFL, f

= w

−a

| consists of two linear

pieces and a single breakpoint a

. w

i,0

= −w

i,1

and a

i,0

= a

. Moreover, both d

i,i+1

and d

i+1,i

are

equal to λ. Similar to the HL-algorithm, we let the

indices i

,...,i

be the sorted indices of n break-

points such that a

< a

i,2

< . .. < a

. In the graph

or G

), the weights of the source adjacent arc

(s,i) ∈ A

and the sink adjacent arc (i,t) ∈ A

are

s,i

(

−w

if a

< a

0 otherwise

, c

i,t

(

0 if a

< a

otherwise

The arcs (i,i + 1) and (i + 1, i) have weights equal to

λ, for i ∈ {1,2,...,n − 1}. The HL-algorithm com-

putes the optimal solution of problem (4) by solving

for the minimum cuts of graphs G

,...,G

At iteration k of the HL-algorithm, we update the

graph from G

k−1

to G

and compute the minimum cut

solution (S

), which allows us to discover the op-

timal solutions of variables whose nodes are in T

but

not in T

k−1

, as explained last section. In PoS where

we consider all λ ∈ [0,∞), the minimum cut solutions

of G

for different values of λ may be different. To

reﬂect the dependence on λ, we use the notations of

ICPRAM 2024 - 13th International Conference on Pattern Recognition Applications and Methods

110

(λ) and (S

(λ),T

(λ)) for the graph and its mini-

mum cut solutions at step k of the algorithm. We may

revert to the notations of G

and (S

) when we dis-

cuss the graph for any λ in general.

The overall idea of PoS is as follow: at the ﬁrst

iteration of the algorithm, we start with the graph

and the range of all possible values of λ, [0, ∞),

denoted by Λ

, on which all values of λ result in

the same minimum cut (S

(λ),T

(λ)) where S

(λ) =

{1,...,n} and T

(λ) =

0. We then update the graph

to G

and partition Λ

into different ranges of λ such

that λ in the same range has the same minimum cut

solution (S

(λ),T

(λ)). The resulting set of ranges of

λ for G

is denoted Λ

. For conciseness, we call each

range of values of λ a λ-range.

At iteration k, we start with a set of ranges of

λ-values, denoted by Λ

k−1

, spanning [0,∞), such

that each λ-range Λ ∈ Λ

k−1

corresponds to a par-

ticular minimum cut solution (S

k−1

(Λ),T

k−1

(Λ)) of

k−1

(Λ). We update G

k−1

to G

, and partition each

λ-range in Λ

k−1

into smaller λ-ranges according to

the solutions of (S

). These λ-ranges that are off-

spring of ranges in Λ

k−1

together form the set Λ

The set Λ

is then passed onto the next iteration.

We explain the ﬁrst iteration of PoS in Subsection

3.1. Subsection 3.2 describes the iteration k. In Sub-

section 3.2.1, we show the computation of the mini-

mum cuts of G

given the minimum cut solutions of

k−1

. Subsection 3.2.2 completes the iteration k by

using the minimum cut solutions to partition λ-ranges

from the previous iteration. We conclude with the out-

puts upon the completion of PoS in Subsection 3.3.

3.1 First Iteration of the Path of

Solutions Algorithm

At the beginning of the ﬁrst iteration, it is clear that

the minimum cut (S

(λ),T

(λ)) of G

(λ), shown in

Figure 2, is always ({1,2,...,n},{}) for any λ ∈ [0,∞)

since all nodes are connected to s and none is con-

nected to t. We let Λ

denote the set of λ-ranges

such that for λ

and λ

that are in the same λ-range,

we have (S

(λ

),T

(λ

)) = (S

(λ

),T

(λ

)). Since

there is only one solution of (S

(λ),T

(λ)) across all

λ ∈ [0, ∞), we have Λ

= {[0, ∞)}.

The ﬁrst step of the algorithm is to determine

the minimum cuts (S

(λ),T

(λ)) of G

(λ), in Figure

3. We consider the outcomes of the cut capacity in

terms of λ, for each possible case of T

(λ) accord-

ing to Lemma 2.5, which states that T

(λ)\T

(λ) must

form an s-interval containing i

in G

(λ). Since all

nodes in G

(λ) are s-nodes for all λ, possible T

(λ)

are node intervals with the form of [i

] where

∈ [i

] ⊆ [1,n]. We group them into cases such

that the coefﬁcient of λ in the cut capacity of each

case is different from other cases.

Case 1. T

(λ) =

0. None of the nodes changes its

status from s to t. C(S

(λ),T

(λ)) = w

Case 2. T

(λ) = [1, n].

C(S

(λ),T

(λ)) = (

∑

i=1

) − w

Case 3. T

(λ) consists of some, but not all, nodes.

There are three subcases.

Case 3.1. T

(λ) = [1, i

] where i

≤ i

< n.

C(S

(λ),T

(λ)) = λ + min

≤i

(

∑

i=1

) − w

The optimal T

(λ) under Case 3.1 is [1,i

], when

= i

, with C(S

(λ),T

(λ)) = λ +

∑

−1

i=1

Case 3.2. T

(λ) = [i

,n] where 1 < i

≤ i

C(S

(λ),T

(λ)) = λ + min

1<i

≤i

(

∑

i=i

) − w

The optimal T

(λ) under Case 3.2 is [i

,n], when

= i

, with C(S

(λ),T

(λ)) = λ +

∑

i=i

Case 3.3. T

(λ) = [i

] where 1 < i

≤

and i

≤ i

< n. C(S

(λ),T

(λ)) = 2λ +

min

1<i

≤i

(

∑

i=i

) − w

. The optimal T

(λ) under

Case 3.3 is {i

} or [i

], when i

= 1 and i

= 1.

C(S

(λ),T

(λ)) = 2λ.

The computations of the optimal cuts in Case 3.1,

3.2 and 3.3. rely on the fact that {w

}

i=1

are positive.

Notice that the comparison between solutions with

similar coefﬁcients of λ in the cut capacity is indepen-

dent of λ. For instance, Case 3.1, with the cut capac-

ity λ+

∑

−1

i=1

, and Case 3.2, with the cut capacity of

λ +

∑

i=i

, have the same coefﬁcient of λ, which

1, in their cut capacities. Determining which one is

more optimal is independent of λ. We only need to

compare the constant terms of the two cases.

Let T

denote the cut with the smallest term con-

stant among the cuts whose coefﬁcients of λ are equal

to ρ, and c

denote the corresponding minimum con-

stant terms. The cut capacity due to the sink set T

then equal to ρλ + c

. If multiple cut solutions have

the smallest constant, we select the cut solution that

results in the largest source set possible, due to the

maximal source set requirement.

For ρ = 0, we compare Case 1 and Case 2. c

min(w

∑

i=1

) − w

). If w

≤ (

∑

i=1

) − w

)

then T

0 (Case 1), otherwise, T

= [1,n] (Case

2). Notice that when w

= (

∑

i=1

)−w

), we select

0 (Case 1) rather than [1,n] (Case 2) since T

0 always results in a larger source set.

For ρ = 1, we compare Case 3.1 and Case

3.2. c

= min(

∑

−1

i=1

∑

i=i

). If

∑

−1

i=1

∑

i=i

then T

= [1, i

] (Case 3.1). If

∑

−1

i=1

∑

i=i

then T

= [i

,n] (Case 3.2). However,

∑

−1

i=1

∑

i=i

, we select the one, between

Case 3.1 and 3.2, that results in a larger source set.

Path of Solutions for Fused Lasso Problems

111

Figure 2: G

(λ) for ℓ

-PFL (1).

Figure 3: G

(λ) for ℓ

-PFL (1).

For ρ = 2, there is only one case, Case 3.3, whose

coefﬁcients of λ in cut capacity is 2. The cut capacity

of Case 3.3, T

(λ) = {i

}, is 2λ. Hence, c

= 0 and

= {i

With the notations of T

and c

, we compare the

capacities between cases of T

(λ) listed above and

give the optimal solution of T

(λ) for different λ in-

tervals as follow:

If T

0 then

(λ) =











, if λ ≥

and λ ≥ c

− c

, if λ < c

− c

and λ > c

if λ ≤ c

and λ <

(5)

However, if T

= [1, n] then

(λ) =











, if λ >

and λ > c

− c

, if λ ≤ c

− c

and λ > c

if λ ≤ c

and λ ≤

(6)

The derivation of both (5) and (6) comes from the

comparison of three cut capacities: c

, λ + c

and

2λ, of which we select the smallest one for different

ranges of λ. The difference between (5) and (6) is

a result of the fact that we always select the solution

with the maximal source set. Hence, for λ such that

and T

are equally optimal, we select the one with

a larger source set. When T

0, T

is always pre-

ferred. When T

= [1,n], T

is always preferred. The

comparison for such situation between T

and T

, and

between T

and T

can be done similarly.

Note that when 2c

> c

, the second case in (5)

does not exist since c

> c

−c

. Hence, when T

and 2c

> c

, the solution of T

(λ) is

(λ) =

(

, if λ ≥

otherwise

(5.1)

Otherwise, when T

0 and 2c

≤ c

(λ) =











, if λ ≥ c

− c

, if c

< λ < c

− c

if λ ≤ c

(5.2)

For the case where T

= [1,n], we can write the so-

lutions of T

(λ) for when 2c

> c

and 2c

≤ c

a similar way as (5.1) and (5.2), which we omit here

due to the space limit.

Suppose the given weights {w

}

i=1

and break-

points {a

}

i=1

result in the case where T

0 and

≤ c

, in which the solution of T

(λ) is (5.2).

We divide the λ-range [0,∞) into three λ-ranges:

[0,c

],(c

− c

) and [c

− c

,∞). We denote this

set of λ-ranges by Λ

, i.e. Λ

= {[0, c

],(c

−

),[c

− c

,∞)}. Each λ-range in Λ

has a different

solution of T

(λ), according to (5.2).

When appropriate, we may use the notation of

(Λ) in place of T

(λ) for λ ∈ Λ ∈ Λ

to emphasize

the fact that T

(λ) for all λ ∈ Λ are identical.

The ﬁrst iteration of PoS ends here as we obtain

the set Λ

and the minimum cut solutions of G

for λ-

ranges in Λ

. In the next iteration of PoS, we solve for

the minimum cuts of the graph G

for each λ-range

in Λ

, which is then partitioned further according to

their minimum cut solutions. The subsequent itera-

tions of PoS will be elaborated in the next section.

Figure 4 illustrates Λ

,Λ

and Λ

which are

the sets of λ-ranges that lead to different minimum

cut solutions. Each λ-range in Λ

corresponds

to a particular sequence of minimum cut solution

(λ),T

(λ)}.

The process continues until we complete the it-

eration n of the algorithm and obtain the set of

λ-ranges Λ

. Each λ-range in Λ

corresponds

to a particular sequence of minimum cut solution

(λ),T

(λ),...,T

(λ)}, which corresponds to a par-

ticular solution of x

∗

(λ).

ICPRAM 2024 - 13th International Conference on Pattern Recognition Applications and Methods

112

Figure 4: λ-Ranges according to the minimum cut solutions of G

(λ), G

(λ) and G

(λ).

3.2 Iteration k of the Path of Solutions

Algorithm

Iteration k of the algorithm involves similar steps and

has the same goal as in the ﬁrst iteration, that is to

partition each λ-range from the previous iteration ac-

cording to their respective minimum cut solutions.

3.2.1 Finding the Minimum Cuts of G

(λ)

By the end of iteration k − 1, we have found the set

of λ-ranges, Λ

k−1

, and (S

k−1

(Λ),T

k−1

(Λ) for each

Λ ∈ Λ

k−1

. Our goal in the iteration k is to ﬁnd

the minimum cut solution (S

(λ),T

(λ)) of the graph

(λ) for λ in each Λ ∈ Λ

k−1

. This minimum cut so-

lution might contain several cases depending on the

value of λ, like the result in (5). These different cases

will then divide the λ-range Λ into smaller λ-ranges.

If i

∈ T

k−1

(λ), then according to Lemma 2.4, the

minimum cut solution (S

(λ),T

(λ)) is then equal to

k−1

(λ),T

k−1

(λ)). A more interesting and less trivial

case is when i

∈ S

k−1

(λ).

We discussed toward the end of Section 2 that,

based on Lemma 2.5, ﬁnding the minimum cut

(λ),T

(λ)) of G

(λ) is equivalent to ﬁnding the

node status change interval in S

k−1

(λ) that move to

(λ). According to Lemma 2.6, given the maxi-

mal s-interval containing i

in G

k−1

(λ), [i

], we

can ﬁnd the node status change interval containing i

], by solving problem (3) deﬁned on G

(λ), re-

stated here:

min

]

C({s},[i

]) +C([i

]\[i

],{t}) (3)

+ C([i

]\[i

],[i

k2]

] ∪ {i

− 1,i

+ 1})

s.t. i

∈[i

] ⊆ [i

] or [i

] =

For each Λ ∈ Λ

k−1

, we know its max-

imal s-interval containing i

, [i

], from

k−1

(Λ),T

k−1

(Λ)). To solve for the node sta-

tus change interval [i

], we write out the cost (3)

for different feasible solutions of [i

] ⊆ [i

for general λ ∈ [0, ∞). Similar to the computation

of T

(λ) in Subsection 3.1, the cost functions for

different possible cases of [i

] may have different

coefﬁcients of λ. These different coefﬁcients of λ

Figure 5: Example when 1 < i

< i

≤ i

< i

n. Here, we only show arcs of weight λ, which connect

consecutive nodes.

in (3) are contingent upon the positions of i

and

], as well as [i

Consider, for example, the most typical case of the

maximal s-interval containing i

, [i

], such that

∈ (1,i

) and i

∈ (i

,n). Different cases of [i

]

that are feasible solutions to problem (3) consist of:

Case 1. i

∈ (i

] and i

∈ [i

), illustrated

in Figure 5. The cost function (3) due to this case

of [i

] can be written as 4λ +

∑

i=i

· (a

∑

−1

i=i

· (a

≤ a

∑

i=i

· (a

≤ a

The ﬁrst term is the cost due to arcs of weight λ.

There are 4 such arcs of weights λ that connect nodes

in [i

]\[i

], which remain in the source set

(λ), to other nodes in the sink set T

(λ). The co-

efﬁcients of λ is then equal to 4. These 4 arcs are

displayed in bold in Figure 5.

The second term is the sum of the weights w for

nodes in [i

] that are still connected to the node s.

The sum of the third term and the fourth term is the

sum of the weights w for nodes in [i

− 1] ∪ [i

1,i

] that are connected to the node t.

The optimal node status change interval for Case

1 is argmin

]:i

∈(i

],i

∈[i

)

∑

i=i

· (a

∑

−1

i=i

· (a

≤ a

∑

i=i

· (a

≤ a

Case 2. There are two subcases.

Case 2.1. i

= i

∈ [i

). The cost is 2λ +

∑

i=i

· (a

> a

) +

∑

i=i

· (a

≤ a

Case 2.2. i

∈ (i

],i

= i

. The cost is equal to

2λ +

∑

i=i

· (a

> a

) +

∑

−1

i=i

· (a

≤ a

Case 3. [i

] is an empty interval, that is,

] ⊆ S

(λ). The cost function can be written as

2λ +

∑

i=i

· (a

≤ a

Case 4. i

= i

. The cost is

∑

i=i

· (a

Path of Solutions for Fused Lasso Problems

113

Let [i

]

denote the node status change in-

terval with the smallest cost among all s-intervals in

] whose coeffcients of λ in the cost are ρ, for

ρ = 0,2,4. This is a similar notation as in Subsection

3.1. For ρ = 0, there is only one case, that is Case 4.

For ρ = 2, we compare Case 2 and 3. For ρ = 4, there

is only one case, Case 1. We do not provide explicit

forms of [i

]

here due to the space limit.

With the same method as presented in Subsection

3.1, we can write out the optimal solution of [i

]

in the following form:

](λ) =











]

, if λ > τ

]

, if τ

≤ λ ≤ τ

]

if λ ≤ τ

(7)

(see the next paragraph for the solution when λ =

) where τ

and τ

are the threshold values that

split λ ∈ [0,∞) into ranges. τ

and τ

can be ob-

tained by comparing the objective functions (3) due

to [i

]

,[i

]

and [i

]

Remark here that when λ = τ

, we select between

]

and [i

]

the solution that results in a

larger source set, as it depends on the problem in-

stance. When λ = τ

, we always prefer [i

]

over

]

= [i

We see for this example of [i

], where i

∈

(1,i

) and i

∈ (i

,n), that there can be up to 3 differ-

ent solutions, as shown in (7), due to 3 different co-

efﬁcients of λ in the cost function across all possible

cases of [i

]. These 3 different solutions corre-

spond to 3 λ-ranges, [0,τ

],[τ

,τ

] and (τ

,∞).

For other types of the maximal s-interval contain-

ing i

, [i

], we list all possible coefﬁcients of λ

and the number of solutions of the node status change

interval, [i

], in Table 1.

We have shown here the procedure to ﬁnd the node

status change interval [i

] from G

k−1

(λ) to G

(λ),

for any λ ≥ 0, given a maximal s-interval containing

, [i

], which can be obtained from the minimum

cut (S

k−1

(λ),T

k−1

(λ)). Next, we show how we apply

this procedure to partition each λ-ranges in Λ

k−1

and

get the set Λ

to complete the iteration k of PoS.

3.2.2 Partitioning of λ-Ranges

The implication of the procedure in Subsection 3.2.1

is that λ-ranges in Λ

k−1

that have the same maximal

s-interval containing i

in G

k−1

also have the same

node status change interval solution. Therefore, the ﬁ-

nal step of the iteration k of PoS is to group λ-ranges

in Λ

k−1

based on their maximal s-intervals contain-

ing i

. For each group, with a particular maximal s-

intervals containing i

, we compute for the node status

change interval, like in Subsection 3.2.1, and apply

the corresponding solution to λ-ranges in that group.

Table 1: Different types of the maximal s-intervals [i

]

containing i

in the minimum cut of G

k−1

, considered in

the k-th iteration of PoS. For each type of [i

], we list

all possible coefﬁcients of λ in the objective function (3).

Types of [i

̸= 1, n and . . .

Coef.

of λ

# cases

∈ (1, i

),i

∈ (i

,n) 0, 2, 4 3

∈ (1, i

),i

= n 0, 1, 2, 3 4

= 1, i

∈ (i

,n) 0, 1, 2, 3 4

= 1, i

= n 0, 1, 2 3

= i

= n 0, 2 2

= 1, i

= i

0, 2 2

= i

∈ (i

,n) 0, 1 2

∈ (1, i

),i

= i

0, 1 2

= i

0, 2 2

Types of [i

∈ {1, n} and . . .

Coef.

of λ

# cases

∈ {1, n} 0, 1 2

= i

= 1, i

∈ (1, n) 0, 1, 2 3

∈ (1, n), i

= i

= n 0, 1, 2 3

For example, suppose the λ-ranges depicted by

rounded rectangular boxes in Figure 6(b) are λ-ranges

in Λ

k−1

that have the same maximal s-interval con-

taining i

, and their node status change interval solu-

tion, in the form similar to (7), is shown in Figure

6(a). τ

,τ

and τ

are threshold values that divide

[0,∞) into 4 segments, each corresponds to a partic-

ular node status change interval solution. Here, we

show the case with 4 solutions, separated by 3 thresh-

olds, which is the maximum number possible (see Ta-

ble 1). This is a different case from (7) where we have

3 solutions.

Then, we intersect these λ-ranges, which have the

same maximal s-interval containing i

(Figure 6(b)),

with the ranges of λ from the node status change in-

terval solution (Figure 6(a)). The result is shown in

Figure 6(c) where some λ-ranges are split into mul-

tiple λ-ranges by threshold values of λ that belong

to the ranges. Each of these resulting λ-ranges takes

the node status change interval solution from the seg-

ments of λ-values (Figure 6(a)) that it belongs to, as

written in the rectangular boxes in Figure 6(c). For

each of these λ-ranges, say a λ-range Λ, with the

node status change interval solution [i

](Λ), its

minimum cut in G

is (S

(Λ),T

(Λ)) where T

(Λ) is

k−1

(Λ) ∪ [i

](Λ).

After we perform this step on all groups of λ-

ranges that share the same maximal s-intervals con-

taining i

in G

k−1

, we obtain a set of λ-ranges span-

ning [0,∞), which together form the set Λ

ICPRAM 2024 - 13th International Conference on Pattern Recognition Applications and Methods

114

Figure 6: At iteration k, we group λ-ranges that have the same maximal s-interval containing i

in the minimum cut of G

k−1

together. Let λ-ranges in (b) be an example of such group of λ-ranges for a speciﬁc maximal s-interval containing i

. The

node status change interval solution can be computed accordingly for this s-interval. An example of such solution is illustrated

in (a). This solution (a) is then applied onto λ-ranges in (b), resulting in a set of λ-ranges in (c), with their node status change

intervals solutions written in their respective rectangular boxes.

3.3 Outputs of the PoS Algorithm

The Path of Solutions algorithm divides the range of

nonnegative λ values, [0,∞), into a set of smaller

ranges, which are then subsequently divided into

smaller ranges in each iteration.

At the end of PoS, we obtain the set of λ-ranges,

. We show in Section 5 that the number of λ-

ranges in Λ

is at most n

. For an interval Λ in

, we have an associated sequence of graph cuts

{(S

(Λ),T

(Λ))}

k=0

. We then obtain the solution

∗

(λ) for λ ∈ Λ based on the HL-algorithm by set-

ting x

∗

(λ) to a

, the k-th smallest breakpoint, for the

index k such that i ∈ S

k−1

(Λ) and i ∈ T

(Λ).

PoS solves ℓ

-PFL in O(n

qlogn). We show this

in Section 5 as well.

4 PATH OF SOLUTIONS

ALGORITHM FOR CONVEX

PIECEWISE LINEAR FIDELITY

FUNCTIONS

Here, we address the path of solutions for a general

fused lasso problem with convex piecewise linear ﬁ-

delity function, or PFL (1):

(PFL) min

,...,x

∑

i=1

) + λ

n−1

∑

i=1

− x

i+1

| (1)

The extension from ℓ

-PFL to PFL can be done in

a similar way as the extension of the HL-algorithm

on GIMR from the case of ℓ

-ﬁdelity function to

convex piecewise linear ﬁdelity functions, given by

Hochbaum and Lu (2017). We give a brief descrip-

tion here:

Let the number of breakpoints of f

) be q

and they are denoted by a

i,1

i,2

,...,a

i,q

. The to-

tal number of breakpoints is

∑

i=1

, denoted by q.

We sort the q breakpoints in the ascending order with

the indices (i

, j

),(i

, j

),...(i

, j

) such that a

, j

< . . . < a

, j

The PoS algorithm applied on PFL involves

the computation of the minimum cuts of graphs

,...,G

, for different λ-ranges, which is done in

a similar manner as when PoS is applied on ℓ

-PFL

in Section 3. The difference here comes from the

update of G

k−1

to G

, which is described in Section

2. The differences between the two graphs are the

weights of the source and sink adjacent arcs of node

. In G

, w

s,i

= max(−( f

)

′

),0) and w

max(( f

)

′

),0). The maximal s-interval of inter-

est when ﬁnding the node status change interval in it-

eration k is the maximal s-interval containing i

, sim-

ilar to when we solve ℓ

-PFL.

5 BOUND OF THE NUMBER OF

SOLUTIONS AND TIME

COMPLEXITY

We ﬁrst provide the following lemma and theorem

that lead to a conclusion that the number of different

solutions to PFL across all nonnegative λ is equal to

the number of λ-ranges in Λ

Lemma 5.1. The set of λ-ranges obtained from the

iteration k of PoS, Λ

, has the following property:

for two different nonnegative values of λ, λ

and λ

they belong to the same λ-range in Λ

if and only if

(λ

) = T

(λ

) for all j = 0,1,...,k.

Proof. Suppose T

(λ

) ̸= T

(λ

) for some j ∈

{0,1,...,k}. Let h be the smallest index such that

(λ

) ̸= T

(λ

). We know that h ≥ 1 since T

(λ

) =

(λ

) = {}. It follows that T

(λ

) = T

(λ

) for all

j ≤ h − 1. Hence, λ

and λ

are in the same λ-

range in Λ

h−1

. Since T

(λ

) ̸= T

(λ

) and T

h−1

(λ

) =

h−1

(λ

), the node status change interval going from

h−1

to G

for λ

, which is T

(λ

)\T

h−1

(λ

), and

Path of Solutions for Fused Lasso Problems

115

that for λ

, which is T

(λ

)\T

h−1

(λ

), must be dif-

ferent. After the partition of their common λ-range in

h−1

, they must be in separate λ-ranges in Λ

. Since

any two λ-ranges never merge, they must also be in

different λ-intervals in Λ

. Hence, we have proven

that if λ

and λ

are in the same λ-range in Λ

then

(λ

) = T

(λ

) for all j = 0,1,...,k.

For the opposite direction of the proof, suppose

and λ

are not in the same λ-range in Λ

. Let

h be the largest index such that λ

and λ

are in

the same λ-range in Λ

(there exists such h because

and λ

are in the same range of [0,∞) ∈ Λ

From the proof of the ﬁrst direction, it follows that

(λ

) = T

(λ

). Since λ

and λ

are not in the

same λ-range in Λ

h+1

, it implies that their node

status change intervals, denoted by [i

](λ

) and

](λ

) are different. Hence, T

h+1

(λ

), which is

equal to T

(λ

) ∪ [i

](λ

), and T

h+1

(λ

), which is

equal to T

(λ

) ∪ [i

](λ

), must be different.

Theorem 5.2. The solutions of PFL (1) for two dif-

ferent values of λ, x

∗

(λ

) and x

∗

(λ

), are equal if

and only if λ

and λ

are in the same λ-range in Λ

Proof. It follows from Lemma 5.1 that λ

and λ

are

in the same λ-range in Λ

if and only if T

(λ

) =

(λ

) for all j = 0,1,...,n.

The threshold theorem implies that, if T

(λ

) =

(λ

) for j = 0,...,n then x

∗

(λ

) = x

∗

(λ

If T

(λ

) ̸= T

(λ

) for some j ∈ {0, . . .,n}, then

there is a node i such that i ∈ T

(λ

) but i /∈ T

(λ

) (or

i /∈ T

(λ

) but i ∈ T

(λ

)). Suppose i ∈ T

(λ

) but i /∈

(λ

). It follows that x

∗

(λ

) ≤ a

< x

∗

(λ

). Hence,

∗

(λ

) ̸= x

∗

(λ

The bounds of the number of solutions to ℓ

-PFL

and PFL across all λ ≥ 0 are given in Theorem 5.3

and 5.4.

Theorem 5.3. The number of different optimal so-

lutions of ℓ

-PFL (4) across all λ values is at most

+3n

+2n+2

, which is smaller than n

when n ≥ 4.

Proof. Theorem 5.2 implies that the number of solu-

tions is equal to the number of λ-ranges in Λ

The iteration k of the algorithm involves dividing

λ-ranges in Λ

k−1

into groups based on the maximal s-

interval containing i

in G

k−1

and then solving for the

node status change interval solution for each maximal

s-interval containing i

The node status change interval solution for each

maximal s-interval containing i

may consist of mul-

tiple solutions, similar to how (7) consists of multiple

cases. If there are p cases of the node status change

interval solution, there will be p − 1 λ-threshold val-

ues that cut into the set of λ-ranges that have the same

maximal s-interval containing i

. In Figure 6, for ex-

ample, p = 4. In (7), p = 3. Each threshold con-

tributes to an increase of at most one additional λ-

ranges. In Figure 6(c), there is an increase of 3 λ-

ranges from Figure 6(b) after we use the λ-threshold

values to split the λ-ranges in Figure 6(b).

Our summary in Table 1 shows that the number of

cases of the node status change interval solution for

any type of the maximal s-interval is at most 4. This

implies that for a set of λ-ranges in Λ

k−1

that have

the same maximal s-interval containing i

, there will

be at most 3 threshold values of λ that cut through

them and the increase in the number of λ-ranges in

due to this maximal s-interval can be at most 3.

Since the number of s-intervals that contain i

, or

the number of [i

] such that 1 ≤ i

≤ i

≤ n is

(n − i

+ 1), the increase in the number of λ-ranges

from Λ

k−1

to Λ

is at most 3i

(n − i

+ 1).

= {[0,∞)} has only one λ-range. After n iter-

ations, in Λ

, the number of λ-ranges is at most 1 +

∑

(n −i

+1) = 1+

n(n+1)(n+2)

+3n

+2n+2

For n ≥ 4, this bound is smaller than n

Theorem 5.4. The number of different optimal solu-

tions of PFL (1), with q breakpoints, across all non-

negative λ is at most 1 +

3q(n+1)

, which is smaller

than n

q for n ≥ 7.

Proof. The proof is similar to that of Theorem 5.3. In

each graph update, or each iteration, say iteration k of

the algorithm, the increase in the number of λ-ranges

in Λ

from Λ

k−1

can be at most 3. When solving

PFL with q breakpoints, the number of iterations is q,

as we update the graph from G

to G

, and so forth,

until G

Suppose the k-th smallest breakpoint comes from

the i-th ﬁdelity function, f

. At iteration k, we do the

computation on the maximal s-interval that contains

node i. The increase in the number of λ-ranges in this

iteration is then at most 3i(n − i + 1).

There are q

iterations where we consider the max-

imal s-interval that contains the node i. Hence, this in-

crease of at most 3i(n − i + 1) happens q

times. This

implies that the bound of the number of λ-ranges in

is 1 +

∑

i=1

i(n − i + 1).

When n is an odd number, max

i=1,...,n

i(n − i +

1) =

(n+1)

. For an even number n, max

i=1,...,n

i(n −

i+1) =

+2n

(n+1)

. Hence, the bound of the num-

ber of λ-ranges in Λ

is at most 1 +

∑

i=1

(n+1)

1 +

3q(n+1)

, which is at most n

q for n ≥ 7.

Theorem 5.5. PFL (1) is solved by PoS in

ICPRAM 2024 - 13th International Conference on Pattern Recognition Applications and Methods

116

Table 2: The number of solutions found by our path of solutions algorithm and Gurobi as well as the computation times, taken

from one of the 10 experiments. For Gurobi, we vary the resolution of the list of values of λ that we provide to the solver.

Gurobi

PoS

Resolution (∆) 0.1 0.05 0.025 0.01 0.005 0.0025 0.001 0.0005

# Solutions 24 36 53 77 93 101 104 108 114

Time (s) 0.45 0.82 1.63 4.10 8.19 16.43 40.72 83.55 0.42

Table 3: Average percentage of solutions found by Gurobi compared to those found by our path of solutions algorithm as well

as the average relative computation time of Gurobi compared to the path of solutions algorithm, across 10 runs.

Gurobi

PoS

Resolution (∆) 0.1 0.05 0.025 0.01 0.005 0.0025 0.001 0.0005

Avg % solutions found 23% 34% 48% 69% 81% 88% 95% 97% 100%

Avg runtime compared to PoS 1.22 2.40 4.78 11.92 23.72 47.47 119.21 237.05 1

Figure 7: Percentages of different solutions of ℓ

-PFL

found by Gurobi with varying resolutions, averaged across

10 experiments. Error bars of 1 standard deviation are dis-

played atop the bar plots. PoS can ﬁnd all solutions across

all λ, indicated by the dashed line at the 100% level.

O(n

qlogn) time where q is the number of break-

points of n convex piecewise linear ﬁdelity functions.

Proof. Suppose the k-th smallest breakpoint comes

from the i-th ﬁdelity function, f

. At iteration k of

the algorithm, we compute for the node status change

interval for each group of λ-ranges in Λ

k−1

that have

the same maximal s-interval containing i. There are

i(n − i + 1) such maximal s-intervals. This is also the

number of the node status change interval computa-

tion at iteration k since we do one computation for

one maximal s-interval.

With the same reasoning as in the proof of Theo-

rem 5.4, the number of the computation for node sta-

tus change interval across q iterations is

∑

i=1

i(n −

i + 1), which is bounded by

as shown earlier.

Hochbaum and Lu (2017) provided a method that

performs a computation of the node status change in-

terval in O(log n). By relying on the same method

and the data structures used in that work, the total

time complexity due to the node status change interval

computation throughout the algorithm is O(n

qlogn).

Other steps such as the sorting of the breakpoints

as well as the updates in the data structures are domi-

nated by the node status change interval computation

time. Therefore, the time complexity of PoS in solv-

ing PFL for all λ ≥ 0 is O(n

qlogn).

6 EXPERIMENTS

To demonstrate the performance of PoS in ﬁnding dif-

ferent solutions for all λ ≥ 0, we compare it to Gurobi,

a state-of-the-art linear programming solver (Gurobi

Optimization, LLC, 2023). We implemented PoS in

Python and run both PoS and Gurobi on the same lap-

top with Apple M2 CPU, 16GB RAM.

We evaluate both PoS and Gurobi on ℓ

-PFL (4)

where the number of samples, n, is 100, the weights

,...,w

and the breakpoints a

,...,a

are sampled

from a uniform distribution between 0 and 1.

PoS generates the solutions for all λ. However,

Gurobi solves the problem only for a speciﬁc value of

λ, which needs to be speciﬁed, one at a time. We pro-

vide a list of values of λ to Gurobi: [0, ∆, 2∆, ...,U].

∆ is the resolution for the grid search on the values of

λ, and U is the upper bound of λ. We take the value

of U from the result of PoS. We vary the resolution ∆

and report the result for all resolutions that we tested.

We run the experiment 10 times. In each run, the

weights and breakpoints are resampled. We measure

the computation time of PoS and Gurobi as well as

the number of different solutions that they found. The

results for one particular run, as an example, are re-

ported in Table 2. In this table, we report the number

of different solutions found by both algorithms and

their runtimes. The reported number of solutions for

Path of Solutions for Fused Lasso Problems

117

PoS is the true number of different solutions since it

ﬁnds all solutions for all nonnegative λ.

As shown in Table 2, Gurobi, with a low resolu-

tion (large ∆), fails to ﬁnd some solutions that cor-

respond to λ not included in the list. With ∆ = 0.1,

both PoS and Gurobi have about the same compu-

tation time. However, Gurobi found only 24 differ-

ent solutions while PoS was able to produce all 114

different solutions. As ∆ decreases, the list of λ for

Gurobi becomes ﬁner. Gurobi is able to ﬁnd more

different solutions, but it also takes longer time. At

the resolution of 0.0005, Gurobi found 108 solutions,

which is close to the total number of different solu-

tions. However, it takes Gurobi more than 80 seconds

while PoS can achieve a better result in 0.42 seconds.

In Table 3 and Figure 7, we report the percent-

ages of solutions found by Gurobi compared to PoS,

averaged across 10 runs. A reported number of 100%

implies that we ﬁnd all possible solutions. We also re-

port the average relative runtime of Gurobi compared

to PoS, where the reported number of a for Gurobi

implies the runtime that is a times of that of PoS.

When we set the resolution for Gurobi at 0.1,

Gurobi and PoS take approximately the same amount

of time. However, with this resolution, Gurobi only

found about 23% of all possible solutions while PoS

can ﬁnd all of them. At the ﬁnest resolution included

in the experiment, which is 0.0005, Gurobi can ﬁnd

97% of all solutions. However, it takes more than 200

times longer than PoS.

7 CONCLUSIONS

We provide in this work an efﬁcient minimum cut-

based algorithm called Path of Solutions, or PoS, that

generates the solutions of the convex piecewise linear

fused lasso problem (1) for all values of the tradeoff

parameter, λ. PoS is a more efﬁcient alternative to

the traditional method of parameter tuning, in which

a single parameter value is evaluated one at a time. In

traditional parameter tuning, we might unknowingly

overlook some important values or test two redundant

values that lead to the same result. These challenges

are overcome by PoS.

In addition to the algorithm design, we provide

both the time complexity of the algorithm and the

bound of the number of different solutions across all

nonnegative λ in terms of the number of variables and

the number of breakpoints in the ﬁdelity functions.

We demonstrate the efﬁciency of PoS via a set of

experiments in comparison with Gurobi, a state-of-

the-art solver for mathematical programming. PoS

is capable of generating all solutions for all λ while

Gurobi found only a small fraction in the same

amount of time. To ﬁnd all solutions, Gurobi requires

more than a factor of 200 of the time that PoS needs.

Regarding future directions of this work, we are

interested in extending our algorithm to solve other

variants of the fused lasso problems as well as con-

ducting more experiments on synthetic and real data.

ACKNOWLEDGEMENTS

This research was supported in part by AI Institute

NSF Award 2112533.

REFERENCES

Eilers, P. H. and De Menezes, R. X. (2005). Quan-

tile smoothing of array cgh data. Bioinformatics,

21(7):1146–1153.

Gallo, G., Grigoriadis, M. D., and Tarjan, R. E. (1989). A

fast parametric maximum ﬂow algorithm and applica-

tions. SIAM Journal on Computing, 18(1):30–55.

Gurobi Optimization, LLC (2023). Gurobi Optimizer Ref-

erence Manual.

Hochbaum, D. S. (2001). An efﬁcient algorithm for image

segmentation, markov random ﬁelds and related prob-

lems. Journal of the ACM (JACM), 48(4):686–701.

Hochbaum, D. S. (2008). The pseudoﬂow algorithm: A new

algorithm for the maximum-ﬂow problem. Operations

research, 56(4):992–1009.

Hochbaum, D. S. and Lu, C. (2017). A faster algorithm

solving a generalization of isotonic median regression

and a class of fused lasso problems. SIAM Journal on

Optimization, 27(4):2563–2596.

Hoeﬂing, H. (2010). A path algorithm for the fused lasso

signal approximator. Journal of Computational and

Graphical Statistics, 19(4):984–1006.

Storath, M., Weinmann, A., and Unser, M. (2016). Ex-

act algorithms for lˆ1-tv regularization of real-valued

or circle-valued signals. SIAM Journal on Scientiﬁc

Computing, 38(1):A614–A630.

Tibshirani, R. J. (2011). The solution path of the generalized

lasso. Stanford University.

Wang, L., Gordon, M. D., and Zhu, J. (2006). Regularized

least absolute deviations regression and an efﬁcient al-

gorithm for parameter tuning. In Sixth International

Conference on Data Mining (ICDM’06), pages 690–

700. IEEE.

ICPRAM 2024 - 13th International Conference on Pattern Recognition Applications and Methods

118