Path Following in the Exact Penalty Method of Convex Programming

1.3 MB·English

by  Hua Zhou

Path Following in the Exact Penalty Method
of Convex Programming
Hua Zhou Kenneth Lange
Department of Statistics Departments of Biomathematics,
2
1 2311 Stinson Drive Human Genetics, and Statistics
0 North Carolina State University University of California
2
Raleigh, NC 27695-8203 Los Angeles, CA 90095-1766

n
E-mail: hua [email protected] E-mail: [email protected]
a
J
January 18, 2012
7
1

]
A
Abstract
N
Classical penalty methods solve a sequence of unconstrained problems that put greater and greater
.
h stress on meeting the constraints. In the limit as the penalty constant tends to , one recovers the
t constrained solution. In the exact penalty method, squared penalties are replace∞d by absolute value
a
m penalties, and the solution is recovered for a finite value of the penalty constant. In practice, the kinks
inthepenaltyandtheunknownmagnitudeofthepenaltyconstantpreventwideapplicationoftheexact
[
penaltymethod in nonlinear programming. Inthis article, we examineastrategy of pathfollowing con-

1 sistent with theexact penalty method. Instead of performing optimization at a single penalty constant,
v wetracethesolutionasacontinuousfunctionofthepenaltyconstant. Thus,pathfollowingstartsatthe
3 unconstrained solution and follows the solution path as the penalty constant increases. In the process,
9 thesolution path hits, slides along, and exitsfrom thevarious constraints. For quadratic programming,
5
the solution path is piecewise linear and takes large jumps from constraint to constraint. For a gen-
3
eral convexprogram, thesolution path is piecewise smooth, and path following operates bynumerically
.
1 solving an ordinary differential equation segment by segment. Our diverse applications to a) projection
0 onto a convex set, b) nonnegative least squares, c) quadratically constrained quadratic programming,
2 d) geometric programming, and e) semidefinite programming illustrate the mechanics and potential of
1 path following. The final detour to image denoising demonstrates the relevance of path following to
:
v regularized estimation in inverse problems. In regularized estimation, one follows the solution path as
i thepenalty constant decreases from a large value.
X
r Keywords: constrained convex optimization, exact penalty, geometric programming, ordinary differ-
a
entialequation,quadratically constrainedquadraticprogramming,regularization, semidefiniteprogram-
ming
1 Introduction
Penaltiesandbarriersarebothpotentdevicesforsolvingconstrainedoptimizationproblems(Boyd and Vandenberghe,
2004;Forsgrenet al.,2002;Luenberger and Ye,2008;Nocedal and Wright,2006;Ruszczyn´ski,2006;Zangwill,
1967). The generalideaistoreplacehardconstraintsbypenalties orbarriersandthenexploitthe well-oiled
machineryforsolvingunconstrainedproblems. Penaltymethodsoperateontheexteriorofthefeasibleregion
1Research supported inpart byUSPHS grants GM53275 and MH59490 to KL,R01 HG006139 to KLand HZ,and NCSU
FRPDgranttoHZ.
1
andbarriermethodsontheinterior. Thestrengthofapenaltyorbarrierisdeterminedbyatuningconstant.
In classical penalty methods, a single global tuning constant is gradually sent to ; in barrier methods, it
∞
is gradually sent to 0. Either strategy generates a sequence of solutions that converges in practice to the
solution of the original constrained optimization problem.
Barriermethods are now generally conceded to offer a better approachto solving convex programs than
penalty methods. Application of log barriers and carefully controlled versions of Newton’s method make it
possibletofollowthecentralpathreliablyandquicklytotheconstrainedminimum(Boyd and Vandenberghe,
2004). Nonetheless, penalty methods should not be ruled out. Augmented Lagrangian methods (Hestenes,
1975)andexactpenaltymethods(Nocedal and Wright,2006)arepotentiallycompetitivewithinteriorpoint
methods for smooth convex programming problems. Both methods have the advantage that the solution
of the constrained problem kicks in for a finite value of the penalty constant. This avoids problems of ill
conditioning as the penalty constant tends to .
∞
The disadvantage of exact penalties over traditional quadratic penalties is lack of differentiability of
the penalized objective function. In the current paper, we argue that this impediment can be finessed by
path following. Our path following method starts at the unconstrained solution and follows the solution
path as the penalty constant increases. In the process, the solution path hits, exits, and slides along the
various constraint boundaries. The path itself is piecewise smooth with kinks at the boundary hitting and
escape times. One advances along the path by numerically solving a differential equation for the Lagrange
multipliers of the penalized problem. In the special case of quadratic programming with affine constraints,
the solution path is piecewise linear, and one can easily anticipate entire path segments (Zhou and Lange,
2011b). This special case is intimately related to the linear complementarity problem (Cottle et al., 1992)
in optimization theory.
Homotopy (continuation) methods for the solution of nonlinear equations and optimization problems
have been pursued for many years and enjoyed a variety of successes (Nocedal and Wright, 2006; Watson,
1986, 0001; Zangwill and Garcia, 1981). To our knowledge, however, there has been no exploration of path
followingasanimplementationoftheexactpenaltymethod. Ourmodestgoalhereistoassessthefeasibility
and versatility of exact path following for constrained optimization. Comparing its performance to existing
methods, particularly the interior point method, is probably best left for later, more practically oriented
papers. In our experience, coding the algorithm is straightforwardin Matlab. The rich numerical resources
of Matlab include differential equation solvers that alert the user when certain events such as constraint
hitting and escape occur.
The rest of the paper is organized as follows. Section 2 briefly reviews the exact penalty method for
optimizationandinvestigatessufficientconditionsforuniquenessandcontinuityofthesolutionpath. Section
2
3 derives the path following strategy for general convex programs, with particular attention to the special
cases of quadratic programming and convex optimization with affine constraints. Section 4 presents various
applicationsofthepathalgorithm. Ourmostelaborateexampledemonstratestherelevanceofpathfollowing
to regularized estimation. The particular problem treated, image denoising, is typical of many inverse
problems in applied mathematics and statistics Zhou and Wu (2011). In such problems one follows the
solution path as the penalty constant decreases. Finally, Section 5 discusses the limitations of the path
algorithm and hints at future generalizations.
2 Exact Penalty Methods
In this paper we consider the convex programming problem of minimizing the convex objective function
f(x) subject to r affine equality constraints g (x) = 0 and s convex inequality constraints h (x) 0. We
i j
≤
will further assumethat f(x)and the h (x) aretwice differentiable. The differential df(x)is the rowvector
j
of partial derivatives of f(x); the gradient f(x) is the transpose of df(x). The second differential d2f(x)
∇
is the Hessian matrix of second partial derivatives of f(x). Similar conventions hold for the differentials of
the constraint functions.
Exact penalty methods (Nocedal and Wright, 2006; Ruszczyn´ski, 2006) minimize the surrogate function
r s
(x) = f(x)+ρ g (x) +ρ max 0,h (x) . (1)
ρ i j
E | | { }
i=1 j=1
X X
This definition of (x) is meaningful regardless of whether the contributing functions are convex. If the
ρ
E
programis convex,then (x)is itself convex. Itis interestingto compare (x) to the Lagrangianfunction
ρ ρ
E E
r s
(x) = f(x)+ λ g (x)+ µ h (x),
i i j j
L
i=1 j=1
X X
which captures the behavior of f(x) near a constrained local minimum y. The Lagrangiansatisfies the sta-
tionarity condition (y)=0; its inequality multipliers µ are nonnegative and satisfy the complementary
j
∇L
slackness conditions µ h (y)=0. In an exact penalty method one takes
j j
ρ > max λ ,..., λ ,µ ,...,µ . (2)
1 r 1 s
{| | | | }
This choice creates the favorable circumstances
(x) (x) for allx
ρ
L ≤ E
(z) f(z) = (z) for all feasiblez
ρ
L ≤ E
(y) = f(y) = (y) fory optimal
ρ
L E
with profound consequences. As the next proposition proves, minimizing (x) is effective in minimizing
ρ
E
f(x) subject to the constraints.
3
Proposition 2.1. Suppose the objective function f(x) and the constraint functions are twice differentiable
and satisfy the Lagrange multiplier rule at the local minimum y. If inequality (2) holds and v∗d2 (y)v >0
L
for every vector v = 0 satisfying dg (y)v = 0 and dh (y)v 0 for all active inequality constraints, then
i j
6 ≤
y furnishes an unconstrained local minimum of (x). For a convex program satisfying Slater’s constraint
ρ
E
qualification and inequality (2), y is a minimum of (x) if and only if y is a minimum of f(x) subject to
ρ
E
the constraints. No differentiability assumptions are required for convex programs.
Proof. The conditions imposed on the quadratic form v∗d2 (y)v are well-known sufficient conditions for
L
a local minimum. Theorems 6.9 and 7.21 of the reference Ruszczyn´ski (2006) prove all of the foregoing
assertions.
Aspreviouslystressed,theexactpenaltymethodturnsaconstrainedoptimizationproblemintoanuncon-
strainedminimizationproblem. Furthermore,incontrasttothequadraticpenaltymethod(Nocedal and Wright,
2006, Section 17.1), the constrained solution in the exact method is achieved for a finite value of ρ. Despite
these advantages, minimizing the surrogate function (x) is complicated. For one thing, it is no longer
ρ
E
globally differentiable. For another, one must minimize (x) along an increasing sequence ρ because the
ρ n
E
Lagrangemultipliers (2) are usually unknownin advance. These hurdles have preventedwide applicationof
exact penalty methods in convex programming.
As a prelude to our derivation of the path following algorithm for convex programs, we record several
properties of (x) that mitigate the failure of differentiability.
ρ
E
Proposition 2.2. The surrogate function (x) is increasing in ρ. Furthermore, (x) is strictly convex
ρ ρ
E E
for one ρ > 0 if and only if it is strictly convex for all ρ > 0. Likewise, it is coercive for one ρ > 0 if and
only if is coercive for all ρ > 0. Finally, if f(x) is strictly convex (or coercive), then all (x) are strictly
ρ
E
convex (or coercive).
Proof. The first assertion is obvious. For the second assertion, consider more generally a finite family
u (x),...,u (x) of convex functions, and suppose a linear combination q c u (x) with positive coeffi-
1 q k=1 k k
cients is strictly convex. It suffices to prove that any other linear combinPation q b u (x) with positive
k=1 k k
coefficients is strictly convex. For any two points x=y and any scalar α (0,1P), we have
6 ∈
u [αx+(1 α)y] αu (x)+(1 α)u (y). (3)
k k k
− ≤ −
Since q c u (x) is strictly convex, strict inequality must hold for at least one k. Hence, multiplying
k=1 k k
inequaPlity (3) by b and adding gives
k
q q q
b u [αx+(1 α)y] < α b u (x)+(1 α) b u (y).
k k k k k k
− −
k=1 k=1 k=1
X X X
4
The third assertion follows from the fact that a convex function is coercive if and only if its restriction to
each half-line is coercive (Bertsekas, 2003, Proposition 3.2.2). Given this result, suppose (x) is coercive,
ρ
E
but ρ∗(x) is not coercive. Then there exists a point x, a direction v, and a sequence of scalars tn tending
E
to such that ρ∗(x+tnv) is bounded above. This requires the sequence f(x+tnv) and each of the
∞ E
sequences g (x+t v) and max 0,h (x+t v) to remain bounded above. But in this circumstance the
i n j n
| | {
sequence (x+t v) also remains bounded above. The final two assertions are obvious.
ρ n
E
3 The Path Following Algorithm
In this section, we take a different point of view. Insteadof minimizing (x) for anincreasing sequence ρ ,
ρ n
E
we study how the solution x(ρ) changes continuously with ρ and devise a path following strategy starting
from ρ = 0. For some finite value of ρ, the path locks in on the solution of the original convex program.
In regularized statistical estimation and inverse problems, the primary goal is to select relevant predictors
rather than to find a constrainedsolution. Thus,the entire solutionpath commands more interestthan any
singlepointalongit(Efron et al.,2004;Osborne et al.,2000;Tibshirani and Taylor,2011;Zhou and Lange,
2011b; Zhou and Wu, 2011). Although our theory will focus on constrained estimation, readers should bear
in mind this second application area of path following.
Thepathalgorithmreliescriticallyonthefirstorderoptimalityconditionthatcharacterizestheoptimum
point of the convex function (y).
ρ
E
Proposition 3.1. For a convex program, a point x = x(ρ) minimizes the function (y) if and only if x
ρ
E
satisfies the stationarity condition
r s
0 = f(x)+ρ s g (x)+ρ t h (x) (4)
i i j j
∇ ∇ ∇
i=1 j=1
X X
for coefficient sets s r and t s . These sets can be characterized as
{ i}i=1 { j}j=1
1 g (x)<0 0 h (x)<0
i j
{− } { }
s [ 1,1] g (x)=0 and t [0,1] h (x)=0. (5)
i ∈− i j ∈ j
 1 gi(x)>0  1 hj(x)>0
{ } { }
At most one point achievesthe minimum of (y) for a given ρ when (y) is strictly convex.
ρ ρ
E E
Proof. According to Fermat’s rule, x minimizes (y) if and only if 0 belongs to the subdifferential ∂ (x)
ρ ρ
E E
of (y). To derivethe subdifferentialdisplayedinequations (4) and(5), one applies the additionandchain
ρ
E
rules of the convex calculus. The sets defining the possible values of s and t are the subdifferentials of the
i j
functions s andt =max t,0 ,respectively. FormoredetailsseeTheorem3.5andancillarymaterialinthe
+
| | { }
book (Ruszczyn´ski, 2006). Finally, it is well knownthat strict convexity guaranteesa unique minimum.
5
To speak coherently of solution paths, one must validate the existence, uniqueness, and continuity of
the solution x(ρ) to the system of equations (1). Uniqueness follows from strict convexity as already noted.
Existence and continuity are more subtle.
Proposition 3.2. If (y) is strictly convex and coercive, then the solution path x(ρ) of equation (1) exists
ρ
E
and is continuous in ρ. If the gradient vectors g (x) : g (x) = 0 h (x) : h (x) = 0 of the active
i i j j
{∇ }∪{∇ }
constraints are linearly independent at x(ρ) for ρ > 0, then the coefficients s (ρ) and t (ρ) are unique and
i j
continuous near ρ as well.
Proof. In accordwith Proposition2.2, we assume that either f(x) is strictly convex and coercive or restrict
our attention to the open interval (0, ). Consider a subinterval [a,b] and fix a point x in the common
∞
domain of the functions (y). The coercivity of (y) and the inequalities
ρ a
E E
[x(ρ)] [x(ρ)] (x) (x)
a ρ ρ b
E ≤ E ≤ E ≤ E
demonstrate that the solution vector x(ρ) is bounded over [a,b]. To prove continuity, suppose that it fails
foragivenρ [a,b]. Thenthereexistsanǫ>0andasequenceρ tendingtoρsuch x(ρ ) x(ρ) ǫfor
n n 2
∈ k − k ≥
all n. Since x(ρ ) is bounded, we can pass to a subsequence if necessary and assume that x(ρ ) converges
n n
to some point y. Taking limits in the inequality [x(ρ )] (x) demonstrates that (y) (x) for
Eρn n ≤ Eρn Eρ ≤ Eρ
all x. Because x(ρ) is unique, we reach the contradictory conclusions y x(ρ) ǫ and y =x(ρ).
2
k − k ≥
Verificationofthesecondclaimisdeferredtopermitfurtherdiscussionofpathfollowing. Theclaimsays
that an active constraint (g (x)=0 or h (x)=0) remains active until its coefficient hits an endpoint of its
i j
subdifferential. Because the solution path is, in fact, piecewise smooth, one can follow the coefficient path
by numerically solving an ordinary differential equation (ODE).
Our path following algorithm works segment-by-segment. Along the path we keep track of the following
index sets
= i:g (x)<0 = j :h (x)<0
E i I j
N { } N { }
= i:g (x)=0 = j :h (x)=0 (6)
E i I j
Z { } Z { }
= i:g (x)>0 = j :h (x)>0
E i I j
P { } P { }
determinedby the signsofthe constraintfunctions. For the sakeofsimplicity, assumethat atthe beginning
of the current segment s does not equal 1 or 1 when i and t does not equal 0 or 1 when j . In
i E j I
− ∈Z ∈Z
other words, the coefficients of the active constraints occur on the interior of their subdifferentials. Let us
showinthis circumstancethatthe solutionpathcanbe extended ina smoothfashion. Ourplanofattackis
to reparameterize by the Lagrange multipliers for the active constraints. Thus, set λ =ρs for i and
i i E
∈Z
6
ω = ρt for j . The multipliers satisfy ρ < λ < ρ and 0 < ω < ρ. The stationarity condition now
j j I i j
∈ Z −
reads
0 = f(x) ρ g (x)+ρ g (x)+ρ h (x)
i i j
∇ − ∇ ∇ ∇
iX∈NE iX∈PE jX∈PI
+ λ g (x)+ ω h (x).
i i j j
∇ ∇
iX∈ZE jX∈ZI
To this we concatenate the constraint equations 0=g (x) for i and 0=h (x) for j .
i E j I
∈Z ∈Z
For convenience now define
dg (x)
UZ(x)=(cid:20) dhZZEI(x) (cid:21), uZ¯(x)=−iX∈NE∇gi(x)+iX∈PE∇gi(x)+jX∈PI∇hj(x).
In this notation the stationarity equation can be recast as
λ
0 = ∇f(x)+ρuZ¯(x)+UtZ(x) ω .
(cid:20) (cid:21)
Under the assumption that the matrix U (x) has full row rank, one can solve for the Lagrange multipliers
Z
in the form
λ
(cid:20)ωZZEI(cid:21) = −[UZ(x)UtZ(x)]−1UZ(x)[∇f(x)+ρuZ¯(x)]. (7)
Hence, the multipliers are unique. Continuity of the multipliers is a consequence of the continuity of the
solution vector x(ρ) and all functions in sight on the right-hand side of equation (7). This observation
completes the proof of Proposition 3.2.
Collectively the stationarity and active constraint equations can be written as the vector equation 0 =
k(x,λ,ω,ρ). To solve for x, λ and ω in terms of ρ, we apply the implicit function theorem (Lange, 2004;
Magnus and Neudecker, 1999). This requires calculating the differential of k(x,λ,ω,ρ) with respect to the
underlyingdependentvariablesx,λ,andω andtheindependentvariableρ. Becausetheequalityconstraints
are affine, a brief calculation gives
d2f(x)+ρ d2h (x)+ ω d2h (x) Ut (x)
∂ k(x,λ,ω,ρ) = j∈PI j j∈ZI j j Z
x,λ,ω U (x) 0
(cid:20) P Z P (cid:21)
∂ k(x,λ,ω,ρ) = uZ¯(x) .
ρ 0
(cid:18) (cid:19)
The matrix ∂ k(x,λ,ω,ρ) is nonsingular when its upper-left block is positive definite and its lower-left
x,λ,ω
block hasfull rowrank (Lange,2010, Proposition11.3.2). Giventhat itis nonsingular,the implicit function
theorem applies, and we can in principle solve for x, λ and ω in terms of ρ. More importantly, the implicit
function theorem supplies the derivative
x
d
λ = ∂ k(x,λ,ω,ρ)−1∂ k(x,λ,ω,ρ), (8)
dρ ZE − x,λ,ω ρ
ω
ZI
 
which is the key to path following. We summarize our findings in the next proposition.
7
Proposition 3.3. Suppose the surrogate function (y) is strictly convex and coercive. If at the point
ρ
E
x(ρ ) the matrix ∂ k(x,λ,ω,ρ) is nonsingular and the coefficient of each active constraints occurs on
0 x,λ,ω
the interior of its subdifferential, then the solution path x(ρ) and Lagrange multipliers λ(ρ) and ω(ρ) satisfy
the differential equation (8) in the vicinity of x(ρ ).
0
Inpracticeonetracesthesolutionpathalongthecurrenttimesegmentuntileitheraninactiveconstraint
becomes activeor the coefficientof anactiveconstrainthits the boundary ofits subdifferential. The earliest
hitting time or escape time over all constraints determines the duration of the current segment. When the
hitting time for an inactive constraint occurs first, we move the constraint to the appropriate active set
E
Z
or andkeeptheotherconstraintsinplace. Similarly,whentheescapetimeforanactiveconstraintoccurs
I
Z
first, we move the constraint to the appropriate inactive set and keep the other constraints in place. In the
second scenario,if s hits the value 1, then we move i to ; If s hits the value 1, then we move i to .
i E i E
− N P
Similar comments apply when a coefficient t hits 0 or 1. Once this move is executed, we commence path
j
followingalongthenewsegment. Pathfollowingcontinuesuntilforsufficientlylargeρ,thesets , ,and
E E
N P
PI are exhausted, uZ¯ = 0, and the solution vector x(ρ) stabilizes. Our previous paper (Zhou and Lange,
2011b) suggests remedies in the very rare situations where escape times coincide.
Pathfollowingsimplifiesconsiderablyintwospecialcases. Considerconvexquadraticprogrammingwith
objectivefunctionf(x)= 1xtAx+btxandequalityconstraintsVx=dandinequalityconstraintsWx e,
2 ≤
where A is positive semi-definite. The exact penalized objective function becomes
s t
1
(x) = xtAx+btx+ρ vtx d +ρ (wtx e ) .
Eρ 2 | i − i| j − j +
i=1 j=1
X X
Since both the equality and inequality constraints are affine, their second derivatives vanish. Both U and
Z
uZ¯ are constant on the current path segment, and the path x(ρ) satisfies
d x A Ut −1 u
λ = Z Z¯ . (9)
dρωZZEI −(cid:18)UZ 0 (cid:19) (cid:18) 0 (cid:19)
 
This implies that the solution path x(ρ) is piecewise linear. Our previous paper (Zhou and Lange, 2011b)
is devoted entirely to this special class of problems and highlights many statistical applications.
On the next rung on the ladder of generality are convex programs with affine constraints. For the exact
surrogate
s t
(x) = f(x)+ρ vtx d +ρ (wtx e ) ,
Eρ | i − i| j − j +
i=1 j=1
X X
the matrix UZ and vector uZ¯ are still constant along a path segment. The relevant differential equation
8
becomes
d λx = d2f(x) UtZ −1 uZ¯ . (10)
dρωZZEI −(cid:18) UZ 0 (cid:19) (cid:18) 0 (cid:19)
 
There are two approachesfor computing the right-handside of equation(10). When A=d2f(x) is positive
definite and B =U has full row rank, the relevant inverse amounts to
Z
A Bt −1 A−1 A−1Bt[BA−1Bt]−1BA−1 A−1Bt[BA−1Bt]−1
= − .
B 0 [BA−1Bt]−1BA−1 [BA−1Bt]−1
(cid:18) (cid:19) (cid:18) − (cid:19)
ThenumericalcostofcomputingtheinversescalesasO(n3)+O( 3). Whend2f(x)isaconstant,theinverse
|Z|
is computed once. Sequentially updating it for different active sets is then conveniently organizedaround
Z
thesweepoperatorofcomputationalstatistics(Zhou and Lange,2011b). Forageneralconvexfunctionf(x),
every time x changes, the inverse must be recomputed. This burden plus the cost of computing the entries
of d2f(x) slow the path algorithm for general convex problems.
Inmanyapplicationsf(x)isconvexbutnotnecessarilystrictlyconvex. Onecancircumventproblemsin
inverting d2f(x) by reparameterizing (Nocedal and Wright, 2006). For the sake of simplicity, suppose that
all of the constraints are affine and that U has full row rank. The set of points x satisfying the active
Z
constraints can be written as x = w +Yy, where w is a particular solution, y is free to vary, and the
columns of Y Rn×(n−|Z|) span the null space of U and hence are orthogonal to the rows of U . Under
Z Z
∈
the null space reparameterization, dx =Y dy. Furthermore,
dρ dρ
Ytd2f(x)Y = d2f(w+Yy)
y
YtuZ¯(x) = ∇y −ρ gi(w+Yy)+ρ gi(w+Yy)+ρ hj(w+Yy) .
h iX∈NE iX∈PE jX∈PI i
It follows that equation (10) becomes
d
dρy = −[Ytd2f(x)Y]−1YtuZ¯
d
dρx = −Y[Ytd2f(x)Y]−1YtuZ¯. (11)
Differentiating equation (7) gives the multiplier derivatives
d λ dx
dρ(cid:20)ωZZEI(cid:21) = −(UZUtZ)−1UZ(cid:18)d2f(x)dρ +uZ¯(cid:19). (12)
The obvious advantage of using equation (11) is that the matrix Ytd2f(x)Y can be nonsingular when
d2f(x) is singular. The computational cost of evaluating the right-hand sides of equations (11) and (12)
is O([n ]3)+O( 3). When n and are small compared to n, this is an improvement over
−|Z| |Z| −|Z| |Z|
the cost O(n3)+O( 3) of computing the right-hand side of equation (8). Balanced against this gain is
|Z|
9
the requirement of finding a basis of the null space of U . Fortunately, the matrix Y is constant over each
Z
path segment and in practice can be computed by taking the QR decomposition of the active constraint
matrix U . At each kink of the solution path, either one constraint enters or one leaves. Therefore,
Z
Z
Y can be sequentially computed by standard updating and downdating formulas (Lawson and Hanson,
1987; Nocedal and Wright, 2006). Which ODE (8) or (11) is preferable depends on the specific application.
When the loss function f(x) is not strictly convex, for example when the number of parameters exceeds
the number of cases in regression, path following requires the ODE (11). Interested readers are referred
to the book (Nocedal and Wright, 2006) for a more extended discussion of range-space versus null-space
optimization methods.
For a general convex program, one can employ Euler’s update
x(ρ+∆ρ) x(ρ) x(ρ)
d
λ(ρ+∆ρ) = λ(ρ) +∆ρ λ(ρ)
    dρ 
ω(ρ+∆ρ) ω(ρ) ω(ρ)
     
to advance the solution of the ODE (8). Euler’s formula may be inaccurate for ∆ρ large. One cancorrectit
byfixingρandperformingonestepofNewton’smethodtore-connectwiththe solutionpath. Thisamounts
to replacing the position-multiplier vector by
x
λ ∂ k(x,λ,ω,ρ)−1k(x,λ,ω,ρ).
x,λ,ω
 −
ω
In practice, it is certainly easierand probably safer to rely on ODE packages such as the ODE45 function in
Matlab to advance the solution of the ODE.
4 Examples of Path Following
Our examplesareintended to illuminate the mechanicsofpath followingandshowcaseits versatility. As we
emphasized in the introduction, we forgo comparisonswith other methods. Comparisons depend heavily on
programming details and problem choices, so a premature study might well be misleading.
Example 4.1. Projection onto the Feasible Region
Finding a feasible point is the initial stage in many convex programs. Dykstra’s algorithm (Dykstra,
1983; Deutsch, 2001) was designed precisely to solve the problem of projecting an exterior point onto the
intersection of a finite number of closed convex sets. The projection problem also yields to our generic path
following algorithm. Consider the toy example of projecting a point b R2 onto the intersection of the
∈
closed unit ball and the closed half space x 0 (Lange, 2004). This is equivalent to solving
1
≥
1
minimize f(x)= x b 2
2k − k
1 1
subject to h (x)= x 2 0, h (x)= x 0.
1 2 1
2k k − 2 ≤ − ≤
10