Journal of Machine Learning Research 8 (2007) 1519-1555

Submitted 7/06; Revised 1/07; Published 7/07

An Interior-Point Method for Large-Scale `1 -Regularized
Logistic Regression
Kwangmoo Koh
Seung-Jean Kim
Stephen Boyd

DENEB 1@ STANFORD . EDU
SJKIM @ STANFORD . EDU
BOYD @ STANFORD . EDUInformation Systems Laboratory
Electrical Engineering Department
Stanford University
Stanford, CA 94305-9510, USA

Editor: Yi Lin

Abstract
Logistic regression with `1 regularization has been proposed as a promising method for feature
selection in classification problems. In this paper we describe an efficient interior-point method
for solving large-scale `1 -regularized logistic regression problems. Small problems with up to athousand or so features and examples can be solved in seconds on a PC; medium sized problems,with tens of thousands of features and examples, can be solved in tens of seconds (assuming somesparsity in the data). A variation on the basic method, that uses a preconditioned conjugate gradientmethod to compute the search step, can solve very large problems, with a million features and examples (e.g., the 20 Newsgroups data set), in a few minutes, on a PC. Using warm-start techniques,a good approximation of the entire regularization path can be computed much more efficiently thanby solving a family of problems independently.
Keywords: logistic regression, feature selection, `1 regularization, regularization path, interiorpoint methods.

1. Introduction
In this section we describe the basic logistic regression problem, the ` 2 - and `1 -regularized versions,and the regularization path. We set out our notation, and review existing methods and literature.
Finally, we give an outline of this paper.
1.1 Logistic Regression
Let x ∈ Rn denote a vector of explanatory or feature variables, and b ∈ {−1, +1} denote the associated binary output or outcome. The logistic model has the form
exp b(wT x + v)
1=,
Prob(b|x) =1 + exp (−b(wT x + v)) 1 + exp (b(wT x + v))
where Prob(b|x) is the conditional probability of b, given x ∈ R n . The logistic model has parametersv ∈ R (the intercept) and w ∈ Rn (the weight vector). When w 6= 0, wT x + v = 0 defines the neutral
hyperplane in feature space, on which the conditional probability of each outcome is 1/2. On theshifted parallel hyperplane wT x + v = 1, which is a distance 1/kwk2 from the neutral hyperplane,c 2007 Kwangmoo Koh, Seung-Jean Kim and Stephen Boyd.

KOH , K IM AND B OYD

the conditional probability of outcome b = 1 is 1/(1 + 1/e) ≈ 0.73, and the conditional probabilityof b = −1 is 1/(1 + e) ≈ 0.27. On the hyperplane wT x + v = −1, these conditional probabilities
are reversed. As wT x + v increases above one, the conditional probability of outcome b = 1 rapidly
approaches one; as wT x + v decreases below −1, the conditional probability of outcome b = −1
rapidly approaches one. The slab in feature space defined by |w T x + v| ≤ 1 defines the ambiguity
region, in which there is substantial probability of each outcome; outside this slab, one outcome ismuch more likely than the other.
Suppose we are given a set of (observed or training) examples,(xi , bi ) ∈ Rn × {−1, +1},

i = 1, . . . , m,assumed to be independent samples from a distribution. We use p log (v, w) ∈ Rm to denote the vector
of conditional probabilities, according to the logistic model,plog (v, w)i = Prob(bi |xi ) =exp(wT ai + vbi ),1 + exp(wT ai + vbi )i = 1, . . . , m,where ai = bi xi . The likelihood function associated with the samples is ∏m
i=1 plog (v, w)i , and the
log-likelihood function is given by
m

m

i=1

i=1

∑ log plog (v, w)i = − ∑ f (wT ai + vbi ),where f is the logistic loss function
f (z) = log(1 + exp(−z)).

(1)This loss function is convex, so the log-likelihood function is concave. The negative of the loglikelihood function is called the (empirical) logistic loss, and dividing by m we obtain the averagelogistic loss,m

lavg (v, w) = (1/m) ∑ f (wT ai + vbi ).
i=1

We can determine the model parameters w and v by maximum likelihood estimation from theobserved examples, by solving the convex optimization problem
minimize lavg (v, w),(2)

with variables v ∈ R and w ∈ Rn , and problem data A = [a1 · · · am ]T ∈ Rm×n and the vector of binary
outcomes b = [b1 · · · bm ]T ∈ Rm . The problem (2) is called the logistic regression problem (LRP).
The average logistic loss is always nonnegative, that is, lavg (v, w) ≥ 0, since f (z) ≥ 0 for any z.
For the choice w = 0, v = 0, we have lavg (0, 0) = log 2 ≈ 0.693, so the optimal value of the LRPlies between 0 and log 2. In particular, the optimal value can range (roughly) between 0 and 1. The
optimal value is 0 only when the original data are linearly separable, that is, there exist w and v suchthat wT xi + v > 0 when bi = 1, and wT xi + v < 0 when bi = −1. In this case the optimal value ofthe LRP (2) is not achieved (except in the limit with w and v growing arbitrarily large). The optimalmvalue is log 2, that is, w = 0, v = 0 are optimal, only if ∑m
i=1 bi = 0 and ∑i=1 ai = 0. (This follows
from the expression for ∇lavg , given in Section 2.1.) This occurs only when the number of positive
1520

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

examples (i.e., those for which bi = 1) is equal to the number of negative examples, and the averageof xi over the positive examples is the negative of the average value of x i over the negative examples.
The LRP (2) is a smooth convex optimization problem, and can be solved by a wide variety ofmethods, such as gradient descent, steepest descent, Newton, quasi-Newton, or conjugate-gradients(CG) methods (see, for example Hastie et al., 2001, § 4.4).
Once we find maximum likelihood values of v and w, that is, a solution of (2), we can predict theprobability of the two possible outcomes, given a new feature vector x ∈ R n , using the associatedlogistic model. For example, when w 6= 0, we can form the logistic classifier
φ(x) = sgn(wT x + v),
where
sgn(z) =

(3)+1 z > 0−1 z ≤ 0,

which picks the more likely outcome, given x, according to the logistic model. This classifier islinear, meaning that the boundary between the two decision outcomes is a hyperplane (defined bywT x + v = 0).
1.2 `2 -Regularized Logistic Regression
When m, the number of observations or training examples, is not large enough compared to n, thenumber of feature variables, simple logistic regression leads to over-fit. That is, the classifier foundby solving the LRP (2) performs perfectly (or very well) on the training examples, but may perform
poorly on unseen examples. Over-fitting tends to occur when the fitted model has many feature
variables with (relatively) large weights in magnitude, that is, w is large.
A standard technique to prevent over-fitting is regularization, in which an extra term that penalizes large weights is added to the average logistic loss function. The ` 2 -regularized logistic
regression problem isnT2minimize lavg (v, w) + λkwk22 = (1/m) ∑m
i=1 f (w ai + vbi ) + λ ∑i=1 wi .

(4)

Here λ > 0 is the regularization parameter, used to control the trade-off between the average logistic
loss and the size of the weight vector, as measured by the ` 2 -norm. No penalty term is imposedon the intercept, since it is a parameter for thresholding the weighted sum w T x in the linear classifier (3). The solution of the `2 -regularized regression problem (4) (which exists and is unique) canbe interpreted in a Bayesian framework, as the maximum a posteriori probability (MAP) estimate
of w and v, when w has a Gaussian prior distribution on Rn with zero mean and covariance λI andv has the (improper) uniform prior on R; see, for example, Chaloner and Larntz (1989) or Jaakkola
and Jordan (2000).
The objective function in the `2 -regularized LRP is smooth and convex, and so (like the ordinary, unregularized LRP) can be minimized by standard methods such as gradient descent, steepest
descent, Newton, quasi-Newton, truncated Newton, or CG methods; see, for example, Luenberger(1984), Lin et al. (2007), Minka (2003), Nocedal and Wright (1999), and Nash (2000). Other methods that have been used include optimization transfer (Krishnapuram and Hartemink, 2005; Zhang
et al., 2004) and iteratively re-weighted least squares (Komarek, 2004). Newton’s method and variants are very effective for small and medium sized problems, while conjugate-gradients and limited
memory Newton (or truncated Newton) methods can handle very large problems. In Minka (2003)
1521

KOH , K IM AND B OYD

the author compares several methods for `2 -regularized LRPs with large data sets. The fastest methods turn out to be conjugate-gradients and limited memory Newton methods, outperforming IRLS,
gradient descent, and steepest descent methods. Truncated Newton methods have been applied tolarge-scale problems in several other fields, for example, image restoration (Fu et al., 2006) andsupport vector machines (Keerthi and DeCoste, 2005). For large-scale iterative methods such as
truncated Newton or CG, the convergence typically improves as the regularization parameter λ isincreased, since (roughly speaking) this makes the objective more quadratic, and improves the conditioning of the problem.
1.3 `1 -Regularized Logistic Regression
More recently, `1 -regularized logistic regression has received much attention. The ` 1 -regularized
logistic regression problem isnTminimize lavg (v, w) + λkwk1 = (1/m) ∑m
i=1 f (w ai + vbi ) + λ ∑i=1 |wi |,(5)where λ > 0 is the regularization parameter. The only difference with ` 2 -regularized logistic regression is that we measure the size of w by its `1 -norm, instead of its `2 -norm. A solution of the`1 -regularized logistic regression must exist, but it need not be unique. Any solution of the ` 1 regularized logistic regression problem (5) can be interpreted in a Bayesian framework as a MAP
estimate of w and v, when w has a Laplacian prior distribution and v has the (improper) uniform
prior. The objective function in the `1 -regularized LRP (5) is convex, but not differentiable (specifically, when any of the weights is zero), so solving it is more of a computational challenge thansolving the `2 -regularized LRP (4).
Despite the additional computational challenge posed by ` 1 -regularized logistic regression,compared to `2 -regularized logistic regression, interest in its use has been growing. The main motivation is that `1 -regularized LR typically yields a sparse vector w, that is, w typically has relatively
few nonzero coefficients. (In contrast, `2 -regularized LR typically yields w with all coefficientsnonzero.) When w j = 0, the associated logistic model does not use the jth component of the feature
vector, so sparse w corresponds to a logistic model that uses only a few of the features, that is,components of the feature vector. Indeed, we can think of a sparse w as a selection of the relevant
or important features (i.e., those associated with nonzero w j ), as well as the choice of the intercept
value and weights (for the selected features). A logistic model with sparse w is, in a sense, simpler or more parsimonious than one with nonsparse w. It is not surprising that ` 1 -regularized LR
can outperform `2 -regularized LR, especially when the number of observations is smaller than thenumber of features (Ng, 2004; Wainwright et al., 2007).
We refer to the number of nonzero components in w as its cardinality, denoted card(w). Thus,`1 -regularized LR tends to yield w with card(w) small; the regularization parameter λ roughly controls card(w), with larger λ typically (but not always) yielding smaller card(w).
The general idea of `1 regularization for the purposes of model or feature selection (or justsparsity of solution) is quite old, and widely used in other contexts such as geophysics (Claerbout
and Muir, 1973; Taylor et al., 1979; Levy and Fullagar, 1981; Oldenburg et al., 1983). In statistics,
it is used in the well-known Lasso algorithm (Tibshirani, 1996) for ` 1 -regularized linear regression,and its extensions such as the fused Lasso (Tibshirani et al., 2005), the grouped Lasso (Kim et al.,2006; Yuan and Lin, 2006; Zhao et al., 2007), and the monotone Lasso (Hastie et al., 2007). The
idea also comes up in signal processing in basis pursuit (Chen and Donoho, 1994; Chen et al.,1522

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

2001), signal recovery from incomplete measurements (Candès et al., 2006, 2005; Donoho, 2006),
and wavelet thresholding (Donoho et al., 1995), decoding of linear codes (Cand ès and Tao, 2005),
portfolio optimization (Lobo et al., 2005), controller design (Hassibi et al., 1999), computer-aided
design of integrated circuits (Boyd et al., 2001), computer vision (Bhusnurmath and Taylor, 2007),sparse principal component analysis (d’Aspremont et al., 2005; Zou et al., 2006), graphical model
selection (Wainwright et al., 2007), maximum likelihood estimation of graphical models (Banerjee
et al., 2006; Dahl et al., 2005), boosting (Rosset et al., 2004), and ` 1 -norm support vector machines(Zhu et al., 2004). A recent survey of the idea can be found in Tropp (2006). Donoho and Elad(2003) and Tropp (2006) give some theoretical analysis of why ` 1 regularization leads to a sparsemodel in linear regression. Recently, theoretical properties of ` 1 -regularized linear regression havebeen studied by several researchers; see, for example, Zou (2006), Zhao and Yu (2006), and Zouet al. (2007).
To solve the `1 -regularized LRP (5), generic methods for nondifferentiable convex problems canbe used, such as the ellipsoid method or subgradient methods (Shor, 1985; Polyak, 1987). Thesemethods are usually very slow in practice, however. (Because ` 1 -regularized LR typically results ina weight vector with (many) zero components, we cannot simply ignore the nondifferentiability ofthe objective in the `1 -regularized LRP (5), hoping to not encounter points of nondifferentiability.)
Another approach is to transform the problem to one with differentiable objective and constraintfunctions. We can solve the `1 -regularized LRP (5), by solving an equivalent formulation, withlinear inequality constraints,
minimize lavg (v, w) + λ1T u
subject to −ui ≤ wi ≤ ui , i = 1, . . . , n,(6)where the variables are the original ones v ∈ R, w ∈ Rn , along with u ∈ Rn . Here 1 denotes thevector with all components one, so 1T u is the sum of the components of u. (To see the equivalencewith the `1 -regularized LRP (5), we note that at the optimal point for (6), we must have u i = |wi |,in which case the objectives in (6) and (5) are the same.) The reformulated problem (6) is a convex
optimization problem, with a smooth objective, and linear constraint functions, so it can be solved
by standard convex optimization methods such as SQP, augmented Lagrangian, interior-point, andother methods. High quality solvers that can directly handle the problem (6) (and therefore, carry out`1 -regularized LR) include for example LOQO (Vanderbei, 1997), LANCELOT (Conn et al., 1992),
MOSEK (MOSEK ApS, 2002), and NPSOL (Gill et al., 1986). These general purpose solvers cansolve small and medium scale `1 -regularized LRPs quite effectively.
Other recently developed computational methods for ` 1 -regularized logistic regression include
the IRLS method (Lee et al., 2006; Lokhorst, 1999), a generalized LASSO method (Roth, 2004)
that extends the LASSO method proposed in Osborne et al. (2000) to generalized linear models,
generalized iterative scaling (Goodman, 2004), bound optimization algorithms (Krishnapuram et al.,2005), online algorithms (Balakrishnan and Madigan, 2006; Perkins and Theiler, 2003), coordinate
descent methods (Friedman et al., 2007; Genkin et al., 2006), and the Gauss-Seidel method (Shevadeand Keerthi, 2003). Some of these methods can handle very large problems (assuming some sparsity
in the data) with modest accuracy. But the additional computational cost required for these methods
to achieve higher accuracy can be very large.
The main goal of this paper is to describe a specialized interior-point method for solving the ` 1 regularized logistic regression problem that is very efficient, for all size problems. In particular our
method handles very large problems, attains high accuracy, and is not much slower than the fastest
1523

KOH , K IM AND B OYD

large-scale methods (conjugate-gradients and limited memory Newton) applied to the ` 2 -regularized
LRP.
Numerical experiments show that our method is as fast as, or faster than, other methods, andreliably provides very accurate solutions. Compared with high-quality implementations of general
purpose primal-dual interior-point methods, our method is far faster, especially for large problems.
Compared with first-order methods such as coordinate descent methods, our method is comparable
in solving large problems with modest accuracy, but is able to solve them with high accuracy withrelatively small additional computational cost.
In this paper we focus on methods for solving the `1 -regularized LRP; we do not discuss thebenefits or advantages of `1 -regularized LR, compared to `2 -regularized LR or other methods formodeling or constructing classifiers for two-class data.
1.4 Regularization Path
Let (v?λ , w?λ ) be a solution for the `1 -regularized LRP with regularization parameter λ. The family ofsolutions, as λ varies over (0, ∞), is called the (`1 -) regularization path. In many applications, theregularization path (or some portion) needs to be computed, in order to determine an appropriate
value of λ. At the very least, the `1 -regularized LRP must be solved for multiple, and often many,
values of λ.
In `1 -regularized linear regression, which is the problem of minimizingkFz − gk22 + λkzk1
over the variable z, where λ > 0 is the regularization parameter, F ∈ R p×n is the covariate matrix,
and g ∈ R p is the vector of responses, it can be shown that the regularization path is piecewise
linear, with kinks at each point where any component of the variable z transitions from zero tononzero, or vice versa. Using this fact, the entire regularization path in a (small or medium size) ` 1 regularized linear regression problem can be computed efficiently (Hastie et al., 2004; Efron et al.,
2004; Rosset, 2005; Rosset and Zhu, 2007; Osborne et al., 2000). These methods are related tonumerical continuation techniques for following piecewise smooth curves, which have been wellstudied in the optimization literature (Allgower and Georg, 1993).
Path following methods have been applied to several problems (Hastie et al., 2004; Park andHastie, 2006a,b; Rosset, 2005). Park and Hastie (2006a) describe an algorithm for (approximately)
computing the entire regularization path for general linear models (GLMs) including logistic regression models. In Rosset (2005), a general path-following method based on a predictor-corrector
method is described for general regularized convex loss minimization problems. Path-following
methods can be slow for large-scale problems, where the number of kinks or events is very large (atleast n). When the number of kinks on the portion of the regularization path of interest is modest,
however, these path-following methods can be very fast, requiring just a small multiple of the effort
needed to solve one regularized problem to compute the whole path (or a portion).
In this paper we describe a fast method for computing a large number of points on the regularization path, using a warm-start technique and our interior-point method. Unlike the methods
mentioned above, our method does not attempt to track the path exactly (i.e., jumping from kinkto kink on the path); it remains efficient even when successive values of λ jump over many kinks.
This is essential when computing the regularization path in a large-scale problem. Our method allows us to compute a large number of points (but much fewer than n, when n is very large) on the1524

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

regularization path, much more efficiently than by solving a family of the problems independently.
Our method is far more efficient than path following methods in computing a good approximationof the regularization for a medium-sized or large data set.
1.5 Outline
In Section 2, we give (necessary and sufficient) optimality conditions, and a dual problem, for the`1 -regularized LRP. Using the dual problem, we show how to compute a lower bound on the suboptimality of any pair (v, w). We describe our basic interior-point method in Section 3, and demonstrate
its performance in Section 4 with small and medium scale synthetic and machine learning benchmark examples. We show that `1 -regularized LR can be carried out within around 35 or so iterations,where each iteration has the same complexity as solving an ` 2 -regularized linear regression problem.
In Section 5, we describe a variation on our basic method that uses a preconditioned conjugate
gradient approach to compute the search direction. This variation on our method can solve verylarge problems, with a million features and examples (e.g., the 20 Newsgroups data set), in under
an hour, on a PC, provided the data matrix is sufficiently sparse.
In Section 6, we consider the problem of computing the regularization path efficiently, at avariety of values of λ (but potentially far fewer than the number of kinks on the path). Using warmstart techniques, we show how this can done much more efficiently than by solving a family ofproblems independently. In Section 7, we compare the interior-point method with several existing
methods for `1 -regularized logistic regression. In Section 8, we describe generalizations of our
method to other `1 -regularized convex loss minimization problems.

2. Optimality Conditions and DualIn this section we derive a necessary and sufficient optimality condition for the ` 1 -regularized LRP,
as well as a Lagrange dual problem, from which we obtain a lower bound on the objective that wewill use in our algorithm.
2.1 Optimality Conditions
The objective function of the `1 -regularized LRP, lavg (v, w)+λkwk1 , is convex but not differentiable,
so we use a first-order optimality condition based on subdifferential calculus (see Bertsekas, 1999,
Prop. B.24 or Borwein and Lewis, 2000, §2). The average logistic loss is differentiable, withm∇v lavg (v, w) = (1/m) ∑ f 0 (wT ai + vbi )bi
i=1m= −(1/m) ∑ (1 − plog (v, w)i )bi
i=1= −(1/m)bT (1 − plog (v, w)),
1525

KOH , K IM AND B OYDandm∇w lavg (v, w) = (1/m) ∑ f 0 (wT ai + vbi )ai
i=1m= −(1/m) ∑ (1 − plog (v, w)i )ai
i=1= −(1/m)AT (1 − plog (v, w)).
The subdifferential of kwk1 is given bywi > 0, {1}{−1} wi < 0,(∂kwk1 )i =
[−1, 1] wi = 0.

The necessary and sufficient condition for (v, w) to be optimal for the ` 1 -regularized LRP (5) is∇v lavg (v, w) = 0,

0 ∈ ∇w lavg (v, w) + λ∂kwk1 ,which can be expressed asbT (1 − plog (v, w)) = 0,
and

(7)


wi > 0, {+λ}T{−λ}wi < 0,(1/m) A (1 − plog (v, w)) i ∈[−λ, λ] wi = 0,
i = 1, . . . , n.

(8)

Let us analyze when a pair of the form (v, 0) is optimal. This occurs if and only ifk(1/m)AT (1 − plog (v, 0))k∞ ≤ λ.

bT (1 − plog (v, 0)) = 0,

The first condition is equivalent to v = log(m+ /m− ), where m+ is the number of training examples
with outcome 1 (called positive) and m− is the number of training examples with outcome −1(called negative). Using this value of v, the second condition becomesλ ≥ λmax = k(1/m)AT (1 − plog (log(m+ /m− ), 0))k∞ .
The number λmax gives us an upper bound on the useful range of the regularization parameter λ:For any larger value of λ, the logistic model obtained from ` 1 -regularized LR has weight zero (andtherefore has no ability to classify). Put another way, for λ ≥ λ max , we get a maximally sparseweight vector, that is, one with card(w) = 0.
We can give a more explicit formula for λmax :λmax = (1/m)
where
b̃i =

m−
m+
ai +
ai = (1/m) X T b̃ ,∑m bi =1
m bi ∑
∞
∞
=−1


m− /m
bi = 1i = 1, . . . , m.
−m+ /m bi = −1,

Thus, λmax is a maximum correlation between the individual features and the (weighted) output
vector b̃. When the features have been standardized, we have ∑m
i=1 xi = 0, so we get the simplified
expression
λmax = (1/m)

∑ xi

bi =1

∞

= (1/m)

1526

∑ xi

bi =−1∞.

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

2.2 Dual Problem
To derive a Lagrange dual of the `1 -regularized LRP (5), we first introduce a new variable z ∈ Rm ,as well as new equality constraints zi = wT ai + vbi , i = 1, . . . , m, to obtain the equivalent problem
minimize (1/m) ∑m
i=1 f (zi ) + λkwk1
subject to zi = wT ai + vbi , i = 1, . . . , m.

(9)

Associating dual variables θi ∈ R with the equality constraints, the Lagrangian ismL(v, w, z, θ) = (1/m) ∑ f (zi ) + λkwk1 + θT (−z + Aw + bv).
i=1

The dual function ism
inf L(v, w, z, θ) = (1/m) inf ∑ ( f (zi ) − mθi zi ) + inf λkwk1 + θT Aw + inf θT bv

v,w,z

z

=



v

w

i=1

∗
T
−(1/m) ∑m
i=1 f (−mθi ) kA θk∞ ≤ λ,−∞otherwise,bT θ = 0,

where f ∗ is the conjugate of the logistic loss function f : (−y) log(−y) + (1 + y) log(1 + y), −1 < y < 0∗0
y = −1 or y = 0
f (y) = sup (yu − f (u)) =u∈R∞,otherwise.

For general background on convex duality and conjugates, see, for example, Boyd and Vandenberghe (2004, Chap. 5) or Borwein and Lewis (2000).
Thus, we have the following Lagrange dual of the `1 -regularized LRP (5):maximize G(θ)
subject to kAT θk∞ ≤ λ,wherebT θ = 0,(10)m

G(θ) = −(1/m) ∑ f ∗ (−mθi )i=1

is the dual objective. The dual problem (10) is a convex optimization problem with variable θ ∈ R m ,and has the form of an `∞ -norm constrained maximum generalized entropy problem. We say that
θ ∈ Rm is dual feasible if it satisfies kAT θk∞ ≤ λ, bT θ = 0.
From standard results in convex optimization we have the following.
• Weak duality. Any dual feasible point θ gives a lower bound on the optimal value p ? of the(primal) `1 -regularized LRP (5):G(θ) ≤ p? .
(11)
• Strong duality. The `1 -regularized LRP (5) satisfies a variation on Slater’s constraint qualification, so there is an optimal solution of the dual (10) θ? , which satisfies
G(θ? ) = p? .
In other words, the optimal values of the primal (5) and dual (10) are equal.
1527

KOH , K IM AND B OYD

We can relate a primal optimal point (v? , w? ) and a dual optimal point θ? to the optimality
conditions (7) and (8). They are related by
θ? = (1/m)(1 − plog (v? , w? )).
We also note that the dual problem (10) can be derived starting from the equivalent problem (6),by introducing new variables zi (as we did in (9)), and associating dual variables θ+ ≥ 0 for theinequalities w ≤ u, and θ− ≥ 0 for the inequalities −u ≤ w. By identifying θ = θ+ − θ− we obtain
the dual problem (10).
2.3 Suboptimality Bound
We now derive an easily computed bound on the suboptimality of a pair (v, w), by constructing adual feasible point θ̄ from an arbitrary w. Define v̄ asv̄ = arg min lavg (v, w),v(12)that is, v̄ is the optimal intercept for the weight vector w, characterized by b T (1 − plog (v̄, w)) = 0.
Now, we define θ̄ asθ̄ = (s/m)(1 − plog (v̄, w)),
(13)
where the scaling constant s is
s = min mλ/kAT (1 − plog (v̄, w))k∞ , 1 .

Evidently θ̄ is dual feasible, so G(θ̄) is a lower bound on p? , the optimal value of the `1 -regularized
LRP (5).
To compute the lower bound G(θ̄), we first compute v̄. This is a one-dimensional smooth convex
optimization problem, which can be solved very efficiently, for example, by a bisection method on
the optimality conditionbT (1 − plog (v, w)) = 0,
since the lefthand side is a monotone function of v. Newton’s method can be used to ensure extremely fast terminal convergence to v̄. From v̄, we compute θ̄ using (13), and then evaluate thelower bound G(θ̄).
The difference between the primal objective value of (v, w), and the associated lower boundG(θ̄), is called the duality gap, and denoted η:η(v, w) = lavg (v, w) + λkwk1 − G(θ̄)
T∗= (1/m) ∑m
i=1 f (w ai + vbi ) + f (−mθi ) + λkwk1 .

(14)

We always have η ≥ 0; and (by weak duality (11)) the point (v, w) is no more than η-suboptimal.
At the optimal point (v? , w? ), we have η = 0.

3. An Interior-Point Method
In this section we describe an interior-point method for solving the ` 1 -regularized LRP (5), in theequivalent formulation
minimize lavg (v, w) + λ1T u
subject to −ui ≤ wi ≤ ui , i = 1, . . . , n,with variables w, u ∈ Rn and v ∈ R.
1528

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

3.1 Logarithmic Barrier and Central Path
The logarithmic barrier for the bound constraints −ui ≤ wi ≤ ui is
n

n

n

i=1

i=1

i=1

Φ(w, u) = − ∑ log(ui + wi ) − ∑ log(ui − wi ) = − ∑ log(u2i − w2i ),with domain
dom Φ = {(w, u) ∈ Rn × Rn | |wi | < ui , i = 1, . . . , n}.
The logarithmic barrier function is smooth and convex. We augment the weighted objective function
by the logarithmic barrier, to obtainφt (v, w, u) = tlavg (v, w) + tλ1T u + Φ(w, u),where t > 0 is a parameter. This function is smooth, strictly convex, and bounded below, and so
has a unique minimizer which we denote (v? (t), w? (t), u? (t)). This defines a curve in R × Rn × Rn ,parametrized by t, called the central path. (See Boyd and Vandenberghe, 2004, Chap. 11 for moreon the central path and its properties.)
With the point (v? (t), w? (t), u? (t)) we associateθ? (t) = (1/m)(1 − plog (v? (t), w? (t))),which can be shown to be dual feasible. (Indeed, it coincides with the dual feasible point θ̄ constructed from w? (t) using the method of Section 2.3.) The associated duality gap satisfies
lavg (v? (t), w? (t)) + λkw? (t)k1 − G(θ? (t)) ≤ lavg (v? (t), w? (t)) + λ1T u? (t) − G(θ? (t)) = 2n/t.
In particular, (v? (t), w? (t)) is no more than 2n/t-suboptimal, so the central path leads to an optimal
solution.
In a primal interior-point method, we compute a sequence of points on the central path, for an
increasing sequence of values of t, using Newton’s method to minimize φt (v, w, u), starting from thepreviously computed central point. A typical method uses the sequence t = t 0 , µt0 , µ2t0 , . . ., whereµ is between 2 and 50 (see Boyd and Vandenberghe, 2004, §11.3). The method can be terminatedwhen 2n/t ≤ ε, since then we can guarantee ε-suboptimality of (v ? (t), w? (t)). The reader is referred
to Nesterov and Nemirovsky (1994), Wright (1997), and Ye (1997) for more on (primal) interiorpoint methods.
3.2 A Custom Interior-Point Method
Using our method for cheaply computing a dual feasible point and associated duality gap for any(v, w) (and not just for (v, w) on the central path, as in the general case), we can construct a custom
interior-point method that updates the parameter t at each iteration.
C USTOM I NTERIOR - POINT M ETHOD FOR `1 - REGULARIZED LR.
given tolerance ε > 0, line search parameters α ∈ (0, 1/2), β ∈ (0, 1)
Set initial values. t := 1/λ, v := log(m+ /m− ), w := 0, u := 1.
repeat1529

KOH , K IM AND B OYD

1. Compute search direction.
∆v
Solve the Newton system ∇2 φt (v, w, u)  ∆w  = −∇φt (v, w, u).
∆u
2. Backtracking line search. Find the smallest integer k ≥ 0 that satisfies∆v
φt (v + βk ∆v, w + βk ∆w, u + βk ∆u) ≤ φt (v, w, u) + αβk ∇φt (v, w, u)T  ∆w .
∆u
3. Update. (v, w, u) := (v, w, u) + βk (∆v, ∆w, ∆u).
4. Set v := v̄, the optimal value of the intercept, as in (12).
5. Construct dual feasible point θ from (13).
6. Evaluate duality gap η from (14).
7. quit if η ≤ ε.
8. Update t.
This description is complete, except for the rule for updating the parameter t, which will be
described below. Our choice of initial values for v, w, u, and t can be explained as follows. The
choice w = 0 and u = 1 seems to work very well, especially when the original data are standardized.
The choice v = log(m+ /m− ) is the optimal value of v when w = 0 and u = 1, and the choice t = 1/λminimizes k(1/t)∇φt (log(m+ /m− ), 0, 1)k2 . (In any case, the choice of the initial values does not
greatly affect performance.) The construction of a dual feasible point and duality gap, in steps 4–6,
is explained in Section 2.3. Typical values for the line search parameters are α = 0.01, β = 0.5,
but here too, these parameter values do not have a large effect on performance. The computational
effort per iteration is dominated by step 1, the search direction computation.
There are many possible update rules for the parameter t. In a classical primal barrier method, t
is held constant until φt is (approximately) minimized, that is, k∇φt k2 is small; when this occurs, t
is increased by a factor typically between 2 and 50. More sophisticated update rules can be foundin, for example, Nesterov and Nemirovsky (1994), Wright (1997), and Ye (1997).
The update rule we propose is(max {µ min{tˆ,t},t} , s ≥ smin
(15)
t :=
t,
s < sminwhere tˆ = 2n/η, and s = βk is the step length chosen in the line search. Here µ > 1 and s min ∈ (0, 1]
are algorithm parameters; we have found good performance with µ = 2 and s min = 0.5.
To explain the update rule (15), we first give an interpretation of tˆ. If (v, w, u) is on the central
path, that is, φt is minimized, the duality gap is η = 2n/t. Thus tˆ is the value of t for which theassociated central point has the same duality gap as the current point. Another interpretation is that
if t were held constant at t = tˆ, (v, w, u) would converge to (v? (tˆ), w? (tˆ), u? (tˆ)), at which point theduality gap would be exactly η.
We use the step length s as a crude measure of proximity to the central path. When the current
point is near the central path, that is, φt is nearly minimized, we have s = 1; far from the central
path, we typically have s  1. Now we can explain the update rule (15). When the current pointis near the central path, as judged by s ≥ smin and tˆ ≈ t, we increase t by a factor µ; otherwise, we
keep t at its current value.
1530

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

We can give an informal justification of convergence of the custom interior-point algorithm. (Aformal proof of convergence would be quite long.) Assume that the algorithm does not terminate.
Since t never decreases, it either increases without bound, or converges to some value t¯. In thefirst case, the duality gap η converges to zero, so the algorithm must exit. In the second case, thealgorithm reduces (roughly) to Newton’s method for minimizing φt¯. This must converge, whichmeans that (v, w, u) converges to (v? (t¯), w? (t¯), u? (t¯)). Therefore the duality gap converges to η̄ =2n/t¯. A basic property of Newton’s method is that near the solution, the step length is one. At thelimit, we therefore havet¯ = max {µ min{2n/η̄, t¯}, t¯} = µt¯,which is a contradiction since µ > 1.
3.3 Gradient and Hessian
In this section we give explicit formulas for the gradient and Hessian of φt . The gradient g =∇φt (v, w, u) is given by
g1
g =  g2  ∈ R2n+1 ,
g3
where

g1 = ∇v φt (v, w, u) = −(t/m)bT (1 − plog (v, w)) ∈ R,


g2

g3


2w1 /(u21 − w21 )..
n= ∇w φt (v, w, u) = −(t/m)AT (1 − plog (v, w)) + ∈R ,.
2wn /(u2n − w2n )
2u1 /(u21 − w21 )..
n= ∇u φt (v, w, u) = tλ1 − ∈R .
.
2un /(u2n − w2n )The Hessian H = ∇2 φt (v, w, u) is given bytbT D0 b
tbT D0 A
0
H =  tAT D0 b tAT D0 A + D1 D2  ∈ R(2n+1)×(2n+1) ,
0
D2
D1


where
D0 = (1/m) diag( f 00 (wT a1 + vb1 ), . . . , f 00 (wT am + vbm )),
D1 = diag 2(u21 + w21 )/(u21 − w21 )2 , . . . , 2(u2n + w2n )/(u2n − w2n )2 ,

D2 = diag −4u1 w1 /(u21 − w21 )2 , . . . , −4un wn /(u2n − w2n )2 .

Here, we use diag(z1 , . . . , zm ) to denote the diagonal matrix with diagonal entries z1 , . . . , zm , wherezi ∈ R, i = 1, . . . , m. The Hessian H is symmetric and positive definite.
1531

KOH , K IM AND B OYD

3.4 Computing the Search Direction
The search direction is defined by the linear equations (Newton system)




 T
∆v
g1
tb D0 btbT D0 A
0 tAT D0 b tAT D0 A + D1 D2   ∆w  = −  g2  .
g3
∆u
0
D2
D1
We first eliminate ∆u to obtain the reduced Newton system


∆v
Hred
= −gred ,
∆w

(16)

where
Hred =

tbT D0 b
tbT D0 A
T
T
tA D0 b tA D0 A + D3

,

gred =



g1
g2 − D2 D−1
1 g3



,

D3 = D1 − D2 D−1
1 D2 .

Once this reduced system is solved, ∆u can be recovered as∆u = −D−1
1 (g3 + D2 ∆w).
Several methods can be used to solve the reduced Newton system (16), depending on the relative
sizes of n and m and the sparsity of the data A.
3.4.1 M ORE E XAMPLES THAN F EATURES
We first consider the case when m ≥ n, that is, there are more examples than features. We form
Hred , at a cost of O(mn2 ) flops (floating-point operations), then solve the reduced system (16) by
Cholesky factorization of Hred , followed by back and forward substitution steps, at a cost of O(n 3 )flops. The total cost using this method is O(mn2 + n3 ) flops, which is the same as O(mn2 ) when
there are more examples than features.
When A is sufficiently sparse, the matrix tAT D0 A + D3 is sparse, so Hred is sparse, with a dense
first row and column. By exploiting sparsity in forming tAT D0 A + D3 , and using a sparse Cholesky
factorization to factor Hred , the complexity can be much smaller than O(mn2 ) flops (see Boyd andVandenberghe, 2004, App. C or George and Liu, 1981).
3.4.2 F EWER E XAMPLES THAN F EATURES
When m ≤ n, that is, there are fewer examples than features, the matrix Hred is a diagonal matrix
plus a rank m + 1 matrix, so we can use the Sherman-Morrison-Woodbury formula to solve thereduced Newton system (16) at a cost of O(m2 n) flops (see Boyd and Vandenberghe, 2004, §4.3).
We start by eliminating ∆w from (16) to obtain(tbT D0 b − t 2 bT D0 AS−1 AT D0 b)∆v = −g1 + tbT D0 AS−1 (g2 − D2 D−1
1 g3 ),
where S = tAT D0 A + D3 . By the Sherman-Morrison-Woodbury formula (Golub and Van Loan,
1996, p. 50), the inverse of S is given by
−1 T
−1 T
(1/t)D−1
S−1 = D−1
0 + AD3 A
3 − D3 A

1532

−1

AD−1
3 .

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION−1 T
We can now calculate ∆v via Cholesky factorization of the matrix (1/t)D−1
and two
0 + AD3 A
backsubstitutions (Boyd and Vandenberghe, 2004, App. C). Once we compute ∆v, we can computethe other components of the search direction asT∆w = −S−1 (g2 − D2 D−1
1 g3 + tA D0 b∆v),∆u = −D−1
1 (g3 + D2 ∆w).
The total cost of computing the search direction is O(m2 n) flops. We can exploit sparsity in the−1 T
Cholesky factorization, whenever (1/t)D−1
0 + AD3 A is sufficiently sparse, to reduce the complexity.
3.4.3 S UMMARY
In summary, the number of flops needed to compute the search direction isO(min(n, m)2 max(n, m)),using dense matrix methods. If m ≥ n and AT A is sparse, or m ≤ n and AAT is sparse, we can use(direct) sparse matrix methods to compute the search direction with less effort. In each of thesecases, the computational effort per iteration of the interior-point method is the same as the effort ofsolving one `2 -regularized linear regression problem.

4. Numerical Examples
In this section we give some numerical examples to illustrate the performance of the interior-point
method described in Section 3, using algorithm parameters
α = 0.01,

β = 0.5,

smin = 0.5,

µ = 2,

ε = 10−8 .

(The algorithm performs well for much smaller values of ε, but this accuracy is more than adequate for any practical use.) The algorithm was implemented in both Matlab and C, and run on a
3.2GHz Pentium IV under Linux. The C implementation, which is more efficient than the Matlab
implementation (especially for sparse problems), is available online (www.stanford.edu/˜boyd/
l1_logreg).
4.1 Benchmark Problems
The data are four small or medium standard data sets taken from the UCI machine learning benchmark repository (Newman et al., 1998) and other sources. The first data set is leukemia cancer gene
expression data (Golub et al., 1999), the second is colon tumor gene expression data (Alon et al.,1999), the third is ionosphere data (Newman et al., 1998), and the fourth is spambase data (Newman
et al., 1998).
For each data set, we considered four values of the regularization parameter: λ = 0.5λ max ,λ = 0.1λmax , λ = 0.05λmax , and λ = 0.01λmax . We discarded examples with missing data, andstandardized each data set. The dimensions of each problem, along with the number of interiorpoint method iterations (IP iterations) needed, and the execution time, are given in Table 1. In
reporting card(w), we consider a component wi to be zero when(1/m) AT (1 − plog (v, w)) i ≤ τλ,1533

KOH , K IM AND B OYD

Data
Leukemia(Golub et al., 1999)

Features n
7129

Examples m
38

Colon
(Alon et al., 1999)

2000

62

Ionosphere(Newman et al., 1998)

34

351

Spambase(Newman et al., 1998)

57

4061

λ/λmax
0.5
0.1
0.05
0.01
0.5
0.1
0.05
0.01
0.5
0.1
0.05
0.01
0.5
0.1
0.05
0.01

card(w)
6
14
14
18
7
22
25
28
3
11
14
24
8
28
38
52

IP iterations
37
38
39
37
35
32
33
32
30
29
30
33
31
32
33
36

Time (sec)
0.60
0.62
0.63
0.60
0.26
0.25
0.26
0.25
0.02
0.02
0.02
0.03
0.63
0.66
0.69
0.75

Table 1: Performance of the interior-point method on 4 data sets, each for 4 values of λ.
where τ = 0.9999. This rule is inspired by the optimality condition in (8).
In all sixteen examples, around 35 iterations were required. We have observed this behavior overa large number of other examples as well. The execution times are well predicted by the complexity
order min(m, n)2 max(m, n).
Figure 1 shows the progress of the interior-point method on the four data sets, for the same four
values of λ. The vertical axis shows duality gap, and the horizontal axis shows iteration number,which is the natural measure of computational effort when dense linear algebra methods are used.
The figures show that the algorithm has linear convergence, with duality gap decreasing by a factor
around 1.85 in each iteration.
4.2 Randomly Generated Problems
To examine the effect of problem size on the number of iterations required, we generate 100 random
problem instances for each of 20 values of n, ranging from n = 100 to n = 10000, with m = 0.1n,that is, 10 times more features than examples. Each problem has an equal number of positive
and negative examples, that is, m+ = m− = m/2. Features of positive (negative) examples are
independent and identically distributed, drawn from a normal distribution N (v, 1), where v is inturn drawn from a uniform distribution on [0, 1] ([−1, 0]).
For each of the 2000 data sets, we solve the `1 -regularized LRP for λ = 0.5λmax , λ = 0.1λmax ,and λ = 0.05λmax . The lefthand plot in Figure 2 shows the mean and standard deviation of thenumber of iterations required to solve the 100 problem instances associated with each value of nand λ. It can be seen that the number of iterations required is very near 35, for all 6000 problem
instances.
In the same way, we generate a family of data sets with m = 10n, that is, 10 times more examples
than features, with 100 problem instances for each of 20 values of n ranging from n = 10 to n = 1000,
and for the same 3 values of λ. The righthand plot in Figure 2 shows the mean and standard deviation
1534

duality gap

duality gap

duality gap

10
PSfrag replacements
20
30
0duality gap
40
10`1 -R EGULARIZED L OGISTIC R EGRESSION
PSfrag replacements A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE
50
20
10−10
30
10−8
40
10−6
50
10−4 0
10−10
0
10
10
10−2
10−8
10
500 −210−6 −2
10
−10
10 10
10−4
dualitygap
10−2 −4
10−4
10
0
10
50
iterations
10−10
duality gap
duality gap
10−6
10−6
10−8

10−8
0

10

20

30

40

10−2

50 −2
10
10−10

10−4
50

10−4

duality gap

100

duality gap

100

10−10
10−6

0

10

20

30

40

0

10

20
iterations

30

40

10−6
10−8

10−8
0

10

20
iterations

30

40

Figure 1: Progress of the interior-point method on 4 data sets, showing duality gap versus iteration
number. Top left: Leukemia cancer gene data set. Top right: Colon tumor gene data set.
Bottom left: Ionosphere data set. Bottom right: Spambase data set.
of the number of iterations required to solve the 100 problem instances associated with each valueof n and λ. The results are quite similar to the case with m = 0.1n.

5. Truncated Newton Interior-Point Method
In this section we describe a variation on our interior-point method that can handle very large problems, provided the data matrix A is sparse, at the cost of having a run time that is less predictable.
The basic idea is to compute the search direction approximately, using a preconditioned conjugate
gradients (PCG) method. When the search direction in Newton’s method is computed approximately, using an iterative method such as PCG, the overall algorithm is called a conjugate gradient
Newton method, or a truncated Newton method (Ruszczynski, 2006; Dembo and Steihaug, 1983).
Truncated Newton methods have been applied to interior-point methods (see, for example, Vandenberghe and Boyd, 1995 and Portugal et al., 2000).
5.1 Preconditioned Conjugate Gradients
The PCG algorithm (Demmel, 1997, §6.6) computes an approximate solution of the linear equations
Hx = −g, where H ∈ RN×N is symmetric positive definite. It uses a preconditioner P ∈ RN×N , also
symmetric positive definite.
1535

number of features n
iterations
102
103 B OYD
KOH , K IM AND
104

40
30

iterations

iterations

PSfrag replacements

10
20
30 40
40

20
10
0 2
10

30
20
10

103
number of features n

0 1
10

104

102
number of features n

103

Figure 2: Average number of iterations required to solve 100 randomly generated ` 1 -regularized
LRPs with different problem size and regularization parameter. Left: n = 10m. Right:
n = 0.1m. Error bars show standard deviation.

P RECONDITIONED C ONJUGATE G RADIENTS A LGORITHM
given relative tolerance εpcg > 0, iteration limit Npcg , and x0 ∈ Rk
k := 0, r0 := Hx0 − g, p1 := −P−1 g, y0 := P−1 r0 .
repeatk := k + 1
z := H pk
θk := yTk−1 rk−1 /pTk z
xk := xk−1 + θk pk
rk := rk−1 − θk z
yk := P−1 rkµk+1 := yTk rk /yTk−1 rk−1
pk+1 := yk + µk+1 pk
until krk k2 /kgk2 ≤ εpcg or k = Npcg .

Each iteration of the PCG algorithm involves a handful of inner products, the matrix-vector
product H pk and a solve step with P in computing P−1 rk . With exact arithmetic, and ignoring the
stopping condition, the PCG algorithm is guaranteed to compute the exact solution x = −H −1 g in
N steps. When P−1/2 HP−1/2 is well conditioned, or has just a few extreme eigenvalues, the PCG
algorithm can compute an approximate solution in a number of steps that can be far smaller than N.
Since P−1 rk is computed in each step, we need this computation to be efficient.
5.2 Truncated Newton Interior-Point Method
The truncated Newton interior-point method is the same as the interior-point algorithm described in
Section 3, with the search direction computed using the PCG algorithm.
1536

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

We can compute H pk in the PCG algorithm using


 T
pk1
tb D0 btbT D0 A
0
H pk =  tAT D0 b tAT D0 A + D1 D2   pk2 
0
D2
D1
pk3


bT u
T

A u + D1 pk2  ,=D2 pk2 + D1 pk3

where u = tD0 (bpk1 + Apk2 ) ∈ Rm . The cost of computing H pk is O(p) flops when A is sparsewith p nonzero elements. (We assume p ≥ n, which holds if each example has at least one nonzero
feature.)We now describe a simple choice for the preconditioner P. The Hessian can be written as
H = t∇2 lavg (v, w) + ∇2 Φ(w, u).
To obtain the preconditioner, we replace the first term with its diagonal part, to get


d0 0
0

P = diag t∇2 lavg (v, w) + ∇2 Φ(w, u) =  0 D3 D2  ,0 D 2 D1

(17)

where

d0 = tbT D0 b,

D3 = diag(tAT D0 A) + D1 .

(Here diag(S) is the diagonal matrix obtained by setting the off-diagonal entries of the matrix S
to zero.) This preconditioner approximates the Hessian of tl avg with its diagonal entries, while
retaining the Hessian of the logarithmic barrier. For this preconditioner, P −1 rk can be computed
cheaply as


−1 
d0 0
0
rk1
P−1 rk =  0 D3 D2   rk2 
rk3
0 D 2 D1
rk1 /d0=  (D1 D3 − D22 )−1 (D1 rk2 − D2 rk3 )  ,(D1 D3 − D22 )−1 (−D2 rk2 + D3 rk3 )which requires O(n) flops.
We can now explain how implicit standardization can be carried out. When using standardized
data, we work with the matrix Astd defined in (20), instead of A. As mentioned in Appendix A,
Astd is in general dense, so we should not form the matrix. In the truncated Newton interior-point
method we do not need to form the matrix Astd ; we only need a method for multiplying a vector by
Astd and a method for multiplying a vector by Astd T . But this is easily done efficiently, using thefact that Astd is a sparse matrix (i.e., A) times a diagonal matrix, plus a rank-one matrix; see (20) in
Appendix A.
There are several good choices for the initial point in the PCG algorithm (labeled x 0 in Section 5.1), such as the negative gradient, or the previous search direction. We have found good
performance with both, with a small advantage in using the previous search direction.
1537

KOH , K IM AND B OYDThe PCG relative tolerance parameter εpcg has to be carefully chosen to obtain good efficiency
in a truncated Newton method. If the tolerance is too small, too many PCG steps are needed tocompute each search direction; if the tolerance is too high, then the computed search directions do
not give adequate reduction in duality gap per iteration. We experimented with several methods ofadjusting the PCG relative tolerance, and found good results with the adaptive rule
εpcg = min {0.1, ξη/kgk2 } ,(18)

where g is the gradient and η is the duality gap at the current iterate. Here, ξ is an algorithm
parameter. We have found that ξ = 0.3 works well for a wide range of problems. In other words, wesolve the Newton system with low accuracy (but never worse than 10%) at early iterations, and solve
it more accurately as the duality gap decreases. This adaptive rule is similar in spirit to standard
methods used in inexact and truncated Newton methods (see Nocedal and Wright, 1999).
The computational effort of the truncated Newton interior-point algorithm is the product ofs, the total number of PCG steps required over all iterations, and the cost of a PCG step, whichis O(p), where p is the number of nonzero entries in A, that is, the total number of (nonzero)features appearing in all examples. In extensive testing, we found the truncated Newton interiorpoint method to be very efficient, requiring a total number of PCG steps ranging between a few
hundred (for medium size problems) and several thousand (for large problems). For medium size
(and sparse) problems it was faster than the basic interior-point method; moreover the truncatedNewton interior-point method was able to solve very large problems, for which forming the Hessian
H (let alone computing the search direction) would be prohibitively expensive.
While the total number of iterations in the basic interior-point method is around 35, and nearly
independent of the problem size and problem data, the total number of PCG iterations required by
the truncated Newton interior-point method can vary significantly with problem data and the valueof the regularization parameter λ. In particular, for small values of λ (which lead to large values ofcard(w)), the truncated Newton interior-point method requires a larger total number of PCG steps.
Algorithm performance that depends substantially on problem data, as well as problem dimension,is typical of all iterative (i.e., non direct) methods, and is the price paid for the ability to solve verylarge problems.
5.3 Numerical Examples
In this section we give some numerical examples to illustrate the performance of the truncatedNewton interior-point method. We use the same algorithm parameters for line search, update rule,and stopping criterion as those used in Section 4, and the PCG tolerance given in (18) with ξ =0.3. We chose the parameter Npcg to be large enough (5000) that the iteration limit was neverreached in our experiments; the typical number of PCG iterations was far smaller. The algorithm isimplemented in both Matlab and C, on a 3.2GHz Pentium IV running Linux, except for very large
problems. For very large problems whose data could not be handled on this computer, the method
was run on AMD Opteron 254 with 8GB main memory. The C implementation is available onlineat www.stanford.edu/˜boyd/l1_logreg.
5.3.1 A M EDIUM S PARSE P ROBLEM
We consider the Internet advertisements data set (Newman et al., 1998) with 1430 features and 2359
examples (discarding examples with missing data). The total number of nonzero entries in the data
1538

(d)

PSfrag replacements

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

10−10 −210

10−1010−2

10−4

10−4

duality gap

100

duality gap

100
50

10−6
(d)
(a),(b),(c)

10−8
0

10

20
iterations

30

10−6
10−8
0

40

(d)
(a) (b) (c)
200
400
600
700 1000
cumulative PCG iterations

Figure 3: Progress of the truncated Newton interior-point method on the Internet advertisements
data set with four regularization parameters: (a) λ = 0.5λmax , (b) λ = 0.1λmax , (c) λ =0.05λmax , and (d) λ = 0.01λmax .
matrix A is p = 39011. We standardized the data set using implicit standardization, as explainedin Section 5.2, solving four `1 -regularized LRPs, with λ = 0.5λmax , λ = 0.1λmax , λ = 0.05λmax ,and λ = 0.01λmax . Figure 3 shows the convergence behavior. The lefthand plot shows the duality
gap versus outer iterations; the righthand plot shows duality gap versus cumulative PCG iterations,which is the more accurate measure of computational effort.
The lefthand plot shows that the number of Newton iterations required to solve the problem is
not much more than in the basic interior-point method described in Section 3. The righthand plotshows that the total number of PCG steps is several hundred, and depends substantially on the valueof λ. Thus, the search directions are computed using on the order of ten PCG iterations.
To give a very rough comparison with the direct method applied to this sparse problem, thetruncated Newton interior-point method is much more efficient than the basic interior-point methodthat does not exploit the sparsity of the data. It is comparable to or faster than the basic interior-point
method that uses sparse linear algebra methods, when the regularization parameter is not too small.
5.3.2 A L ARGE S PARSE P ROBLEM
Our next example uses the 20 Newsgroups data set (Lang, 1995). We processed the data set in away similar to Keerthi and DeCoste (2005). The positive class consists of the 10 groups with namesof form sci.*, comp.*, and misc.forsale, and the negative class consists of the other 10 groups. We
used McCallum’s Rainbow program (McCallum, 1996) with the commandrainbow -g 3 -h -s -O 2 -i
to tokenize the (text) data set. These options specify trigrams, skip message headers, no stoplist,
and drop terms occurring fewer than two times. The resulting data set has n = 777811 features(trigrams) and m = 11314 examples (articles). Each example contains an average of 425 nonzero
features. The total number of nonzero entries in the data matrix A is p = 4802169. We standardized
the data set using implicit standardization, as explained in Section 5.2, solving three ` 1 -regularized
LRPs, with λ = 0.5λmax , λ = 0.1λmax , and λ = 0.05λmax . (For the value λ = 0.01λmax , the runtime
is on the order of one hour. This case is not of practical interest, and so not reported here, since1539

λ/λmax
0.5
0.1
0.05

PSfrag replacements

60
80
100
10−10
−8 B OYD
KOH , K IM10
AND
10−6
10−4
10−2
0
card(w) Iterations 10 PCG iterations
9
43
558
(a)
544
60
1036
(b)
2531
58
2090
(c)

Time (sec)
134
256
501

Table 2: Performance of truncated Newton interior-point method on the 20 newsgroup data set(n = 777811 features, m = 11314 examples) for 3 values of λ.

100

10−10
duality gap

duality gap

80
100
100
10−10
10−2
10−4

10−2
10−4
10−6

10−6
10−8
0

(a)
40

20

(c)(b)
60

10−80iterations

(b)
(a)
(c)
1000
500
1500
2000
cumulative PCG iterations

Figure 4: Progress of the truncated Newton interior-point method on the 20 Newsgroups data setfor (a) λ = 0.5λmax , (b) λ = 0.1λmax , and (c) λ = 0.05λmax . Left. Duality gap versus
iterations. Right. Duality gap versus cumulative PCG iterations.
the cardinality of the optimal solution is around 10000 and comparable to the number of examples.)The performance of the algorithm, and the cardinality of the weight vectors, is given in Table 2.
Figure 4 shows the progress of the algorithm, with duality gap versus iteration (lefthand plot), andduality gap versus cumulative PCG iteration (righthand plot).
The number of iterations required to solve the problems ranges between 43 and 60, depending on
λ. The more relevant measure of computational effort is the total number of PCG iterations, whichranges between around 500 and 2000, again, increasing with decreasing λ, which corresponds toincreasing card(w). The average number of PCG iterations, per iteration of the truncated Newton
interior-point method, is around 13 for λ = 0.5λmax , 17 for λ = 0.1λmax , and 36 for λ = 0.05λmax .
(The variance in the number of PCG iterations required per iteration, however, is large.) The running
time is consistent with a cost of around 0.24 seconds per PCG iteration. The increase in runningtime, for decreasing λ, is due primarily to an increase in the average number of PCG iterations
required per iteration, but also from an increase in the overall number of iterations required.
5.3.3 R ANDOMLY G ENERATED P ROBLEMS
We generated a family of 21 data sets, with the number of features n varying from one hundred
to ten million, and m = 0.1n examples. The data were generated using the same general method
described in Section 4.2, but with A sparse, with an average number of nonzero features per example
around 30. Thus, the total number of nonzero entries in A is p ≈ 30m. We standardized the data
1540

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION
PSfrag replacements

λ = 0.1λmax

104
λ = 0.05λmax

runtime in seconds

103
102
101

λ = 0.5λmax

100
10−1
10−2 2
10

103

104
105
number of features n

106

107

Figure 5: Runtime of the truncated Newton interior-point method, for randomly generated sparse
problems, with three values of λ.

set using implicit standardization, as explained in Section 5.2, solving each problem instance forthe three values λ = 0.5λmax , λ = 0.1λmax , and λ = 0.05λmax . The total runtime, for the 63 `1 regularized LRPs, is shown in Figure 5. The plot shows that runtime increases as λ decreases, andgrows approximately linearly with problem size.
We compare the runtimes of the truncated Newton interior-point and the basic interior-point
method using dense linear algebra methods to compute the search direction. Figure 6 shows the results for λ = 0.1λmax . The truncated Newton interior-point method is far more efficient for medium
problems. For large problems, the basic interior-point method fails due to memory limitations, or
extremely long computation times.
By fitting an exponent to the data over the range from n = 320 to the largest problem successfully
solved by each method, we find that the basic interior-point method scales as O(n 2.8 ) (which isconsistent with the basic flop count analysis, which predicts O(n 3 )). For the truncated Newton
interior-point method, the empirical complexity is O(n1.3 ).
When sparse matrix methods are used to compute the search direction in the basic interiorpoint method, we get an empirical complexity of O(n2.2 ) for the Matlab implementation of thebasic interior-point method that uses sparse matrix methods, showing a good efficiency gain overdense methods, for medium scale problems. The C implementation would have the same empiricalcomplexity as the Matlab one with a smaller constant hidden in the O(·) notation.
5.3.4 P RECONDITIONER P ERFORMANCETo examine the effect of the preconditioner (17) on the efficiency of the approximate search direction
computation, we compare the eigenvalue distributions of the Hessian H and the preconditioned
Hessian P−1/2 HP−1/2 , for the colon gene tumor problem (n = 2000 features, m = 62 examples) at
the 15th iterate, in Figure 7. The eigenvalues of the preconditioned Hessian are tightly clustered,1541

KOH , K IM AND B OYD
PSfrag replacements

PSfrag replacements
index
eigenvalues

104

preconditioned(a) Hessian
Hessianruntime in seconds

103
102
101
100
10−1

PSfrag replacements

10−2 2
10

103

1
1000
2000
3000
4000
10−1
100
101
102
103
104
104 105 105
number of features n

(b)

106

107

Figure 6: Runtime of (a) the basic interior-point method and (b) the truncated Newton interior-point
method, for a family of randomly generated sparse problems.
108

104
103

106
Hessianeigenvalues

eigenvalues

Hessian102
101

preconditioned Hessian

preconditioned Hessian

102
100

100
10−1
105 1

104

1000

2000
index

3000

10−2
4000 1010 1

1000

2000
index

3000

4000

Figure 7: Eigenvalue distributions of Hessian and preconditionned Hessian, at the 15th iterate, forthe colon gene tumor problem, for λ = 0.5λmax (left) and λ = 0.05λmax (right).
with just a few extreme eigenvalues, which explains the good performance with relatively few PCG
iterations per iteration (Demmel, 1997, §6.6).

6. Computing the Regularization PathIn this section we consider the problem of solving the `1 -regularized LRP for M values of theregularization parameter λ,λmax = λ1 > λ2 > · · · > λM > 0.
This can be done by applying the methods described above, for each of the M problems. Thisis called a cold-start approach, since each problem is solved independently of the others. This
1542

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

is efficient when multiple processors are used, since the LRPs can be solved simultaneously, on
different processors. But when one processor is used, we can solve these M problems much more
efficiently by solving them sequentially, using the previously computed solution as a starting point
for the next computation. This is called a warm-start approach.
We first note that the solution for λ = λ1 = λmax is (log(m+ /m− ), 0, 0). Since this point does notsatisfy |wi | < ui , it cannot be used to initialize the computation for λ = λ2 . We modify it by adding
a small increment to u to get(v(1) , w(1) , u(1) ) = (log(m+ /m− ), 0, (εabs /(nλ))1),
which is strictly feasible. In fact, it is on the central path with parameter t = 2n/ε abs , and so isεabs -suboptimal. Note that so far we have expended no computational effort.
Now for k = 2, . . . , M we compute the solution (v(k) , w(k) , u(k) ) of the problem with λ = λk , byapplying the interior-point method, with starting point modified to be(vinit , winit , uinit ) = (v(k−1) , w(k−1) , u(k−1) ),and initial value of t set to t = 2n/εabs .
In the warm-start technique described above, the number of grid points, M, is fixed in advance.
The grid points (and M) can be chosen adaptively on the fly, while taking into account the curvatureof the regularization path trajectories, as described in Park and Hastie (2006a).
6.1 Numerical Results
Our first example is the leukemia cancer gene expression data, for M = 100 values of λ, uniformly
distributed on a logarithmic scale over the interval [0.001λ max , λmax ]. (For this example, λmax =0.37.) The left plot in Figure 8 shows the regularization path, that is, w (k) , versus regularization
parameter λ. The right plot shows the number of iterations required to solve each problem from awarm-start, and from a cold-start.
The number of cold-start iterations required is always near 36, while the number of warm-start
iterations varies, but is always smaller, and typically much smaller, with an average value of 3.1.
Thus the computational savings for this example is over 11 : 1.
Our second example is the 20 Newsgroups data set, with M = 100 values of λ uniformly spaced
on a logarithmic scale over [0.05λmax , λmax ]. For this problem we have λmax = 0.12. The top plot
in Figure 9 shows the regularization path. The bottom left plot shows the total number of PCG
iterations required to solve each problem, with the warm-start and cold-start methods. The bottom
right plot shows the cardinality of w as a function of λ.
Here too the warm-start method gives a substantial advantage over the cold-start method, at leastfor λ not too small, that is, as long as the optimal weight vector is relatively sparse. The total runtime
using the warm-start method is around 2.8 hours, and the total runtime using the cold-start method
is around 6.2 hours, so the warm-start methods gives a savings of around 2 : 1. If we consider only
the range from 0.1λmax to λmax , the savings increases to 5 : 1.
We note that for this example, the number of events (i.e., a weight transitioning between zeroand nonzero) along the regularization path is very large, so methods that attempt to track every event
will be very slow.
1543

10 λmax
λmax
−0.3
0
KOH , K IM AND B OYD
0.3
0.6
0.9

PSfrag replacements

cold start
warm start40
iterations required

weights

0.9
0.6
0.3
0
−0.3
10−3 λmax

10−2 λmax

λ10−1 λmax

λmax

cold start
30
20
warm start
10010−3 λmax

10−2 λmax

λ10−1 λmax

λmax

Figure 8: Left. Regularization path for leukemia cancer gene expression data. Right. Iterations
required for cold-start and warm-start methods.

7. Comparison
In this section we compare the performance of our basic and truncated Newton interior-point methods, implemented in C (called l1 logreg), with several existing methods for ` 1 -regularized logistic
regression, We make comparisons with MOSEK (MOSEK ApS, 2002), IRLS-LARS (Lee et al., 2006),
BBR (Genkin et al., 2006), and glmpath (Park and Hastie, 2006a).
MOSEK is a general purpose primal-dual interior-point solver, which is known to be quite efficient compared to other standard solvers. MOSEK can solve ` 1 -regularized LRPs using the separable
convex formulation (9), or by treating the problem as a geometric program (GP) (see Boyd et al.,2006). We used both formulations and report the better results here in each case. MOSEK uses a
stopping criterion based on the duality gap, like our method.
IRLS-LARS alternates between approximating the average logistic loss by a quadratic approximation at the current iterate, and solving the resulting `1 -regularized least squares problem using
the LARS method (Efron et al., 2004) to update the iterate. IRLS-LARS outperforms many existing
methods for `1 -regularized logistic regression including GenLASSO (Roth, 2004), SCGIS (Goodman,2004), Gl1ce (Lokhorst, 1999), and Grafting (Perkins and Theiler, 2003). IRLS-LARS used inour comparison is implemented in Matlab and C, with the LARS portion implemented in C. The
hybrid implementation is called IRLS-LARS-MC and available from http://ai.stanford.edu/
˜silee/softwares/irlslars.htm. We ran it until the primal objective is within tolerance fromthe optimal objective value, which is computed using l1 logreg, with small (10 −12 ) duality gap.
BBR, implemented in C, uses the cyclic coordinate descent method for Bayesian logistic regression. The C implementation is available from http://www.stat.rutgers.edu/˜madigan/BBR/.
The stopping criterion is based on lack of progress, and not on a suboptimality bound or duality gap.
Tightening the tolerance for BBR greatly increased its running time, and only had a minor effect on
the final accuracy.
glmpath uses a path-following method for generalized linear models, including logistic models,and computes a portion of the regularization path. It is implemented in the R environment and available from http://cran.r-project.org/src/contrib/Descriptions/glmpath.html. We compared glmpath to our warm-start method, described in Section 6.
1544

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

PSfrag replacements
λweights0.05λmax
0.1λmax
0.5λmax
λmax
−0.4
−0.2
0
0.2
0.4
λ
PCG iterations

0.4

weights

PSfrag replacements
0.2
PSfrag replacements
λweights0.05λmax
0.1λmax
0.5λmax
λmax
−0.4
−0.2
0
0.2
0.4

0

cold start
warm start

−0.2
λmax
−0.4
0.05λmax

0.05λmax
0.1λmax
0.5λmax
λmax
0
1000λ
2000
3000
3000

0.1λ0.5λmax maxcard(w)

PCG iterations

3000

2000
1000

2000

1000

cold start
warm start

00.05λmax0.1λmax

λ0.5λmax λmax

00.05λmax0.1λmax

λ0.5λmax λmax

Figure 9: Top. Regularization path for 20 newsgroup data. Bottom left. Total PCG iterations required by cold-start and warm-start methods. Bottom right. card(w) versus λ.

1545

KOH , K IM AND B OYD

We report the run times of the methods, using the different stopping criteria described above. We
also report the actual accuracy achieved, that is, the difference between the achieved primal objective
and the optimal objective value (as computed by l1 logreg with duality gap 10 −12 ). However, it
is important to point out that it is very difficult, if not impossible, to carry out a fair comparisonof solution methods, due to the issue of implementation (which can have a great influence on thealgorithm performance), the choice of algorithm parameters, and the different stopping criteria.
Therefore, the comparison results reported below should be interpreted with caution.
We report comparison results using four data sets: the leukemia cancer gene expression andspambase data sets, two dense benchmark data sets used in Section 4.1, the Internet advertisements
data set, the medium sparse data set used in Section 5.3, and the 20 Newsgroups data set, the large
sparse data set used in Section 5.3. When the large 20 Newsgroups data set was standardized, thethree existing solvers could not handle a data set of this size, so it was not standardized. The solvers
could handle the standardized Internet advertisements data set but do not exploit the sparse plus
rank-one structure of the standardized data matrix. Therefore, the Internet advertisements data set
was not standardized as well. For small problems (leukemia and spambase), the solvers were all run
on a 3.2GHz Pentium IV under Linux; for medium and large problem (Internet advertisements and
20 Newsgroups), the solvers were run on an AMD Opteron 254 (with 8GB RAM) under Linux.
The regularization parameter λ can strongly affect the runtime of the methods, including ours.
For each data set, we considered many values of the regularization parameter λ over the interval[0.0001λmax , λmax ] (which appears to cover the range of interest for many standardized problems).
We report here the results with λ = 0.001λmax for the two small benchmark data. The cardinalityof the optimal solution is 21 for the leukemia data and 54 for the spambase data, larger than half ofthe minimum of the number of features and the number of examples. For the two unstandardized
problems, `1 -regularized LR with a regularization parameter in the interval [0.0001λ max , λmax ] yields
a very sparse solution. We used λ = 0.0001λmax for the Internet advertisements and λ = 0.001λmax
for the 20 Newsgroups data. The cardinalities of optimal solutions are relatively very small (19for the Internet advertisements and 247 for the 20 Newsgroups) compared with the sizes of the
problems.
Table 3 summarizes the comparison results for the four problems. Here, the tolerance has a different meaning depending on the method, as described above. As shown in this table, our method
is as fast as, or faster than, existing methods for the small two data sets, but the discrepancy inperformance is not significant, since the problems are small. For the unstandardized Internet advertisements data set, our method is most efficient. MOSEK could not handle the unstandardized large
20 Newsgroups data set, so we compared l1 logreg with BBR and IRLS-LARS-MC on the unstandardized 20 Newsgroups data set, for the regularization parameter λ = 0.001λ max . The truncatedNewton interior-point method solves the problem to the accuracy 2 × 10 −6 in around 100 seconds
and to a higher accuracy with relatively small additional run time; BBR solves them to this ultimate
accuracy in a comparable time, but slows down when attempting to compute a solution with higher
accuracy.
Finally, we compare the runtimes of the warm-start method and glmpath. We consider theleukemia cancer gene expression data set and the Internet advertisements data set as benchmark
examples of small dense and medium sparse problems, respectively. For each data set, the warmstart method finds M points (w(k) , k = 1, . . . , M) on the regularization path, with λ uniformly spaced
on a logarithmic scale over the interval [0.001λmax , λmax ]. glmpath finds an approximation of theregularization path by choosing the kink points adaptively over the same interval. The results are1546

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

ProgramTol.

l1 logreg

10−4
10−8
10−4
10−8
10−4
10−8
10−4
10−11

IRLS-LARS-MC
MOSEK
BBR

Leukemiatime accuracy
0.37 2 × 10−6
0.57 1 × 10−100.75 6 × 10−60.81 6 × 10−11
9 8 × 10−6
10 3 × 10−9
15 1 × 10−7
73 4 × 10−8

Spambasetime accuracy
0.34 2 × 10−50.74 1 × 10−90.29 9 × 10−6
0.37 3 × 10−11
10 8 × 10−8
11 5 × 10−13
39 3 × 10−6
300 1 × 10−11

Internet adv.
time accuracy
0.14 4 × 10−50.27 5 × 10−9
1.2 1 × 10−5
2.5 1 × 10−10
1.3 3 × 10−9
1.4 4 × 10−130.44 1 × 10−7
1.1 1 × 10−12

Newsgroups
time accuracy
90 2 × 10−6
140 2 × 10−10
450 8 × 10−51200 7 × 10−9
- - 140 2 × 10−6
850 1 × 10−6

Table 3: Comparison results with two standardized data sets (leukemia and spambase) and two
unstandardized data sets (Internet advertisements and 20 Newsgroups). The regularization
parameter is taken as λ = 0.001λmax for leukemia and spambase, λ = 0.0001λmax forInternet advertisements and λ = 0.001λmax for 20 Newsgroups.

Data
Leukemia
Internetwarm-start (M = 25)
2.8
5.2

warm-start (M = 100)6.4
13glmpath
1.9
940

Table 4: Regularization path computation time (in seconds) of the warm-start method and glmpath
for standardized leukemia cancer gene expression data and Internet advertisements data.
shown in Table 4. For the small leukemia data set, glmpath is faster than the warm-start method.
The warm-start method is more efficient than glmpath for the Internet advertisements data set, amedium-sized sparse problem. The performance discrepancy is partially explained by the fact that
our warm-start method exploits the sparse plus rank-one structure of the standardized data matrix,
whereas glmpath does not.

8. Extensions and Variations
The basic interior-point method and the truncated Newton variation can be extended to general ` 1 regularized convex loss minimization problems, with twice differentiable loss functions, that havethe form
minimize (1/m) ∑m
i=1 φ(zi ) + λkwk1(19)
subject to zi = wT ai + vbi + ci , i = 1, . . . , m,with variables are v ∈ R, w ∈ Rn , and z ∈ Rm , and problem data ci ∈ R, ai ∈ Rn , and bi ∈ R
determined by a set of given (observed or training) examples(xi , yi ) ∈ Rn × R,

i = 1, . . . , m.

Here φ : R → R is a loss function which is convex and twice differentiable. Prior work related tothe extension includes Park and Hastie (2006a), Rosset (2005), and Tibshirani (1997).
In `1 -regularized (binary) classification, we have yi ∈ {−1, +1} (binary labels), and zi has theform zi = yi (wT xi + v), so we have the form (19) with ai = yi xi , bi = yi , and ci = 0. The associated classifier is given by y = sgn(wT x + v). The loss function φ is small for positive arguments,
1547

KOH , K IM AND B OYD

and grows for negative arguments. When φ is the logistic loss function f in (1), this general ` 1 regularized classification problem reduces to the `1 -regularized LRP. When φ is the convex loss
function φ(u) = − log Φ(u), where Φ is the cumulative distribution function of the standard normal
distribution, this `1 -regularized classification problem reduces to the `1 -regularized probit regression problem. More generally, `1 -regularized estimation problems that arise in generalized linear
models (McCullagh and Nelder, 1989; Hastie et al., 2001) for binary response variables (which include logistic and probit models) can be formulated problems of the form (19); see Park and Hastie(2006a) for the precise formulation.
In `1 -regularized linear regression, we have yi ∈ R, and zi has the form zi = wT xi + v − yi ,which is the difference between yi and its predicted value, wT xi + v. Thus `1 -regularized regression
problems have the form (19) with ai = xi , bi = 1, and ci = −yi . Typically φ is symmetric, with φ(0) =0. When the loss function is quadratic, that is, φ(u) = u2 , the convex loss minimization problem (19)is the `1 -regularized least squares regression problem studied extensively in the literature.
The dual of the `1 -regularized convex loss minimization problem (19) is∗Tmaximize −(1/m) ∑m
i=1 φ (−mθi ) + θ c
TT
subject to kA θk∞ ≤ λ, b θ = 0,

where A = [a1 · · · am ]T ∈ Rm×n , the variable is θ ∈ Rm , and φ∗ is the conjugate of the loss functionφ,φ∗ (y) = sup (yu − φ(u)) .
u∈R

As with `1 -regularized logistic regression, we can derive a bound on the suboptimality of (v, w), byconstructing a dual feasible point θ̄, from an arbitrary w,
 0 T

φ (w a1 + v̄b1 + c1 )..
θ̄ = (s/m)p(v̄, w),p(v̄, w) = ,.
φ0 (wT am + v̄bm + cm )where v̄ is the optimal intercept for the offset w,

mv̄ = arg min(1/m) ∑ φ(wT ai + vbi + ci ),v

i=1

and the scaling constant s is given by s = min mλ/kAT p(v̄, w))k∞ , 1 .
Using this method for cheaply computing a dual feasible point and associated duality gap forany (v, w), we can extend the custom interior-point method for ` 1 -regularized LRPs to general `1 regularized convex (twice differentiable) loss problems.
Other possible extensions include `1 -regularized Cox proportional hazards models (Cox, 1972).
The associated `1 -regularized problem does not have the form (19), but the idea behind the custom
interior-point method for `1 -regularized LRPs can be readily extended. The reader is referred toPark and Hastie (2006a) and Tibshirani (1997) for related work on computational methods for ` 1 regularized Cox proportional hazards models.


1548

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

Acknowledgments
This work is supported in part by the National Science Foundation under grants #0423905 and(through October 2005) #0140700, by the Air Force Office of Scientific Research under grant#F49620-01-1-0365, by MARCO Focus center for Circuit & System Solutions contract #2003CT-888, and by MIT DARPA contract #N00014-05-1-0700. The authors thank Michael Grant,
Trevor Hastie, Honglak Lee, Suin Lee, Mee Young Park, Robert Tibshirani, Yinyu Ye, and Sungroh
Yoon for helpful comments and suggestions. The authors thank anonymous reviewers for helpfulcomments and suggestions (especially on comparison with existing methods).

Appendix A. Standardization
Standardization is a widely used pre-processing step applied to the feature vector, so that each(transformed) feature has zero mean and unit variance (over the examples) (Ryan, 1997). The mean
feature vector is µ = (1/m) ∑m
i=1 xi , and the vector of feature standard deviations σ is defined byσj =m(1/m) ∑ (xi j − µ j )2
i=1

!1/2

,

j = 1, . . . , n,where xi j is the jth component of xi . The standardized feature vector is defined as
xstd = diag(σ)−1 (x − µ).
When the examples are standardized, we obtain the standardized data matrix
Astd = diag(b)(X − 1µT ) diag(σ)−1 = A diag(σ)−1 − bµT diag(σ)−1 ,(20)

where X = [x1 · · · xm ]T . We carry out logistic regression (possibly regularized) using the data matrix
Astd in place of A, to obtain (standardized) logistic model parameters w std , vstd . In terms of theoriginal feature vector, our logistic model isProb(b|x) =1
1=
1 + exp (−b(wstd T xstd + vstd )) 1 + exp (−b(wT x + v))

where
w = diag(σ)−1 wstd ,v = vstd − wstd T diag(σ)−1 µ.

We point out one subtlety here related to sparsity of the data matrix. For small or medium sized
problems, or when the original data matrix A is dense, forming the standardized data matrix A std
does no harm. But when the original data matrix A is sparse, which is the key to efficient solutionof large-scale `1 -regularized LRPs, forming Astd is disastrous, since Astd is in general dense, evenwhen A is sparse.
But we can get around this problem, when working with very large problems, by never actuallyforming the matrix Astd , as explained in Section 5.
1549

KOH , K IM AND B OYD

References
E. Allgower and K. Georg. Continuation and path following. Acta Numerica, 2:1–64, 1993.
U. Alon, N. Barkai, D. Notterman, K. Gish, S. Ybarra, D. Mack, and A. Levine. Broad patternsof gene expression revealed by clustering of tumor and normal colon tissues probed by oligonucleotide arrays. Proceedings of the National Academy of Sciences, 96:6745–6750, 1999.
S. Balakrishnan and D. Madigan. Algorithms for sparse linear classifiers in the massive data setting,2006. Manuscript. Available from www.stat.rutgers.edu/˜madigan/papers/.
O. Banerjee, L. El Ghaoui, A. d’Aspremont, and G. Natsoulis. Convex optimization techniques forfitting sparse Gaussian graphical models. In Proceedings of the 23rd International Conference
on Machine Learning, 2006.
D. Bertsekas. Nonlinear Programming. Athena Scientific, second edition, 1999.
A. Bhusnurmath and C. Taylor. Solving the graph cut problem via ` 1 norm minimization, 2007.
University of Pennsylvania CIS Tech Report number MS-CIS-07-10.
J. Borwein and A. Lewis. Convex Analysis and Nonlinear Optimization. Springer, 2000.
S. Boyd, S.-J. Kim, L. Vandenberghe, and A. Hassibi. A tutorial on geometric programming, 2006.
To appear in Optimization and Engineering.
S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
S. Boyd, L. Vandenberghe, A. El Gamal, and S. Yun. Design of robust global power and groundnetworks. In Proceedings of ACM/SIGDA International Symposium on Physical Design (ISPD),pages 60–65, 2001.
E. Candès, J. Romberg, and T. Tao. Stable signal recovery from incomplete and inaccurate measurements. Communications on Pure and Applied Mathematics, 59(8):1207–1223, 2005.
E. Candès, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal reconstruction
from highly incomplete frequency information. IEEE Transactions on Information Theory, 52(2):489–509, 2006.
E. Candès and T. Tao. Decoding by linear programming. IEEE Transactions on Information Theory,51(12):4203–4215, 2005.
K. Chaloner and K. Larntz. Optimal Bayesian design applied to logistic regression experiments.
Journal of Statistical Planning and Inferenc, 21:191–208, 1989.
S. Chen and D. Donoho. Basis pursuit. In Proceedings of the Twenty-Eighth Asilomar Conferenceon Signals, Systems and Computers, volume 1, pages 41–44, 1994.
S. Chen, D. Donoho, and M. Saunders. Atomic decomposition by basis pursuit. SIAM Review, 43(1):129–159, 2001.
J. Claerbout and F. Muir. Robust modeling of erratic data. Geophysics, 38(5):826–844, 1973.
1550

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

A. Conn, N. Gould, and Ph. Toint. LANCELOT: A Fortran package for large-scale nonlinear
optimization (Release A), volume 17 of Springer Series in Computational Mathematics. SpringerVerlag, 1992.
D. Cox. Regression models and life-tables. Journal of the Royal Statistical Society. Series B, 34(2):
187–220, 1972.
J. Dahl, V. Roychowdhury, and L. Vandenberghe. Maximum likelihood estimation of Gaussian
graphical models: Numerical implementation and topology selection. Submitted. Available fromwww.ee.ucla.edu/˜vandenbe/covsel.html, 2005.
A. d’Aspremont, L. El Ghaoui, M. Jordan, and G. Lanckriet. A direct formulation for sparse PCA
using semidefinite programming, 2005. In L. Saul, Y. Weiss and L. Bottou, editors, Advances in
Neural Information Processing Systems, 17, pp. 41-48, MIT Press.
R. Dembo and T. Steihaug. Truncated-Newton algorithms for large-scale unconstrained optimization. Math. Program., 26:190–212, 1983.
J. Demmel. Applied Numerical Linear Algebra. Society for Industrial and Applied Mathematics,
1997.
D. Donoho. Compressed sensing. IEEE Transactions on Information Theory, 52(4):1289–1306,
2006.
D. Donoho and M. Elad. Optimally sparse representation in general (non-orthogonal) dictionaries
via `1 minimization. Proc. Nat. Aca. Sci., 100(5):2197–2202, March 2003.
D. Donoho, I. Johnstone, G. Kerkyacharian, and D. Picard. Wavelet shrinkage: Asymptopia? J. R.
Statist. Soc. B., 57(2):301–337, 1995.
B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression. Annals of Statistics, 32(2):407–499, 2004.
J. Friedman, T. Hastie T, and R. Tibshirani. Pathwise coordinate optimization, 2007. Manuscript
available from www-stat.stanford.edu/˜hastie/pub.htm.
H. Fu, M. Ng, M. Nikolova, and J. Barlow. Efficient minimization methods of mixed ` 1 -`1 and `2 -`1
norms for image restoration. SIAM Journal on Scientific computing, 27(6):1881–1902, 2006.
A. Genkin, D. Lewis, and D. Madigan. Large-scale Bayesian logistic regression for text categorization, 2006. To appear in Technometrics. Available from www.stat.rutgers.edu/˜madigan/papers/.
A. George and J. Liu. Computer Solution of Large Sparse Positive Definite Systems. Prentice-Hall,
1981.
P. Gill, W. Murray, M. Saunders, and M. Wright. User’s guide for NPSOL (Version 4.0): A FORTRAN package for nonlinear programming. Technical Report SOL 86-2, Operations Research
Dept., Stanford University, Stanford, California 94305, January 1986.
1551

KOH , K IM AND B OYD

G. Golub and C. Van Loan. Matrix Computations, volume 13 of Studies in Applied Mathematics.
John Hopkins University Press, third edition, 1996.
T. Golub, D. Slonim, P. Tamayo, C. Gaasenbeek, J. Mesirov, H. Coller, M. Loh, J. Downing,
M. Caligiuri, C. Bloomfield, and E. Lander. Molecular classification of cancer: Class discovery and class prediction by gene expression monitoring. Science, 286:531–537, 1999.
J. Goodman. Exponential priors for maximum entropy models. In Proceedings of the Annual
Meetings of the Association for Computational Linguistics, 2004.
A. Hassibi, J. How, and S. Boyd. Low-authority controller design via convex optimization. In
Proceedings of the IEEE Conference on Decision and Control, pages 140–145, 1999.
T. Hastie, S. Rosset, R. Tibshirani, and J. Zhu. The entire regularization path for the support vector
machine. Journal of Machine Learning Research, 5:1391–1415, 2004.
T. Hastie, J. Taylor, R. Tibshirani, and G. Walther. Forward stagewise regression and the monotone
lasso. Electronic Journal of Statistics, 1:1–29, 2007.
T. Hastie, R. Tibshirani, and J. Friedman. The Elements of Statistical Learning. Springer Series in
Statistics. Springer-Verlag, New York, 2001. ISBN 0-387-95284-5.
T. Jaakkola and M. Jordan. Bayesian parameter estimation via variational methods. Statistics andComputing, 10:25–37, 2000.
S. Keerthi and D. DeCoste. A modified finite Newton method for fast solution of large scale linear
SVMs. Journal of Machine Learning Research, 6:341–361, 2005.
Y. Kim, J. Kim, and Y. Kim. Blockwise sparse regression. Statistica Sinica, 16:375–390, 2006.
P. Komarek. Logistic Regression for Data Mining and High-Dimensional Classification. PhD thesis,
Carnegie Mellon University, 2004.
B. Krishnapuram, L. Carin, M. Figueiredo, and A. Hartemink. Sparse multinomial logistic regression: Fast algorithms and generalization bounds. IEEE Transactions on Pattern Analysis andMachine Intelligence, 27(6):957–968, 2005.
B. Krishnapuram and A. Hartemink. Sparse multinomial logistic regression: Fast algorithms andgeneralization bounds. IEEE Transactions on Pattern Analysis and Mach. Intelligence, 27(6):957–968, 2005. ISSN 0162-8828.
K. Lang. Newsweeder: Learning to filter netnews. In Proceedings of the Twenty-First International
Conference on Machine learning (ICML), pages 331–339, 1995.
S. Lee, H. Lee, P. Abeel, and A. Ng. Efficient l1 -regularized logistic regression. In Proceedings ofthe 21st National Conference on Artificial Intelligence (AAAI-06), 2006.
S. Levy and P. Fullagar. Reconstruction of a sparse spike train from a portion of its spectrum andapplication to high-resolution deconvolution. Geophysics, 46(9):1235–1243, 1981.
1552

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

C.-J. Lin, R. Weng, and S. Keerthi. Trust region Newton methods for large-scale logistic regression, 2007. To appear in Proceedings of the 24th International Conference on Machine Learning(ICML).
M. Lobo, M. Fazel, and S. Boyd. Portfolio optimization with linear and fixed transaction costs.
Annals of Operations Research, 2005.
J. Lokhorst. The LASSO and generalised linear models, 1999. Honors Project, Department of
Statistics, The University of Adelaide, South Australia, Australia.
D. Luenberger. Linear and Nonlinear Programming. Addison-Wesley, second edition, 1984.
A. McCallum. Bow: A toolkit for statistical language modeling, text retrieval, classification andclustering. Available from www.cs.cmu.edu/˜mccallum/bow, 1996.
P. McCullagh and J. Nelder. Generalized Linear Models. Chapmand & Hall/CRC, second edition,
1989.
T. Minka. A comparison of numerical optimizers for logistic regression, 2003. Technical report.
Available from research.microsoft.com/˜minka/papers/logreg/.
MOSEK ApS. The MOSEK Optimization Tools Version 2.5. User’s Manual and Reference, 2002.
Available from www.mosek.com.
S. Nash. A survey of truncated-Newton methods. Journal of Computational and Applied Mathematics, 124:45–59, 2000.
Y. Nesterov and A. Nemirovsky. Interior-Point Polynomial Methods in Convex Programming, volume 13 of Studies in Applied Mathematics. SIAM, Philadelphia, PA, 1994.
D. Newman, S. Hettich, C. Blake, and C. Merz. UCI repository of machine learning databases,
1998. Available from www.ics.uci.edu/˜mlearn/MLRepository.html.
A. Ng. Feature selection, `1 vs. `2 regularization, and rotational invariance. In Proceedings of theTwenty-First International Conference on Machine learning (ICML), pages 78–85, New York,NY, USA, 2004. ACM Press. ISBN 1-58113-828-5.
J. Nocedal and S. Wright. Numerical Optimization. Springer Series in Operations Research.
Springer, 1999.
D. Oldenburg, T. Scheuer, and S. Levy. Recovery of the acoustic impedance from reflection seismograms. Geophysics, 48(10):1318–1337, 1983.
M. Osborne, B. Presnell, and B. Turlach. A new approach to variable selection in least squares
problems. IMA Journal of Numerical Analysis, 20(3):389–403, 2000.
M.-Y. Park and T. Hastie. An `1 regularization-path algorithm for generalized linear models, 2006a. To appear in Journal of the Royal Statistical Society, Series B. Available fromwww-stat.stanford.edu/˜hastie/pub.htm.
1553

KOH , K IM AND B OYD

M.-Y. Park and T. Hastie. Regularization path algorithms for detecting gene interactions, 2006b.
Manuscript. Available from www-stat.stanford.edu/˜hastie/Papers/glasso.pdf.
S. Perkins and J. Theiler. Online feature selection using grafting. In Proceedings of the Twenty-FirstInternational Conference on Machine learning (ICML), pages 592–599. ACM Press, 2003.
B. Polyak. Introduction to Optimization. Optimization Software, 1987. Translated from Russian.
L. Portugal, M. Resende, G. Veiga, and J. Júdice. A truncated primal-infeasible dual-feasible network interior point method. Networks, 35:91–108, 2000.
S. Rosset. Tracking curved regularized optimization solution paths. In L. Saul, Y. Weiss, and L. Bottou, editors, Advances in Neural Information Processing Systems 17. MIT Press, Cambridge, MA,
2005.
S. Rosset and J. Zhu. Piecewise linear regularized solution paths, 2007. To appear in Annals of
Statistics.
S. Rosset, J. Zhu, and T. Hastie. Boosting as a regularized path to a maximum margin classifier.
Journal of Machine Learning Research, 5:941–973, 2004.
V. Roth. The generalized LASSO. IEEE Transactions on Neural Networks, 15(1):16–28, 2004.
A. Ruszczynski. Nonlinear Optimization. Princeton university press, 2006.
T. Ryan. Modern Regression Methods. Wiley, 1997.
S. Shevade and S. Keerthi. A simple and efficient algorithm for gene selection using sparse logistic
regression. Bioinformatics, 19(17):2246–2253, 2003.
N. Z. Shor. Minimization Methods for Non-differentiable Functions. Springer Series in Computational Mathematics. Springer, 1985.
H. Taylor, S. Banks, and J. McCoy. Deconvolution with the l 1 norm. Geophysics, 44(1):39–52,
1979.
R. Tibshirani. Regression shrinkage and selection via the Lasso. Journal of the Royal StatisticalSociety, Series B, 58(1):267–288, 1996.
R. Tibshirani. The Lasso for variable selection in the Cox model. Statistics in Medicine, 16:385–
395, 1997.
R. Tibshirani, M. Saunders, S. Rosset, and J. Zhu. Sparsity and smoothness via the fused Lasso.
Journal of the Royal Statistical Society Series B, 67(1):91–108, 2005.
J. Tropp. Just relax: Convex programming methods for identifying sparse signals in noise. IEEE
Transactions on Information Theory, 52(3):1030–1051, 2006.
L. Vandenberghe and S. Boyd. A primal-dual potential reduction method for problems involving
matrix inequalities. Math. Program., 69:205–236, 1995.
1554

A N I NTERIOR -P OINT M ETHOD FOR L ARGE -S CALE `1 -R EGULARIZED L OGISTIC R EGRESSION

R. Vanderbei. LOQO User’s Manual — Version 3.10, 1997. Available from www.orfe.princeton.
edu/loqo.
M. Wainwright, P. Ravikumar, and J. Lafferty. High-dimensional graphical model selection using`1 -regularized logistic regression., 2007. To appear in Advances in Neural Information Processing Systems (NIPS) 19. Available from http://www.eecs.berkeley.edu/˜wainwrig/Pubs/
publist.html#High-dimension%al.
S. Wright. Primal-dual Interior-point Methods. Society for Industrial and Applied Mathematics,
Philadelphia, PA, USA, 1997. ISBN 0-89871-382-X.
Y. Ye. Interior Point Algorithms: Theory and Analysis. John Wiley & Sons, 1997.
M. Yuan and L. Lin. Model selection and estimation in regression with grouped variables. Journalof the Royal Statistical Society, Series B, 68(1):49–67, 2006.
Z. Zhang, J. Kwok, and D. Yeung. Surrogate maximization/minimization algorithms for AdaBoostand the logistic regression model. In Proceedings of the Twenty-First International Conference
on Machine learning (ICML), pages 927–934, New York, NY, USA, 2004. ACM Press.
P. Zhao, G. Rocha, and B. Yu. Grouped and hierarchical model selection through composite absolute
penalties, 2007. Tech Report 703. Stat Dept. UCB. Available from www.stat.berkeley.edu/
˜binyu/ps/703.pdf.
P. Zhao and B. Yu. On model selection consistency of Lasso. Journal of Machine Learning Research, 7:2541–2563, 2006.
J. Zhu, S. Rosset, T. Hastie, and R. Tibshirani. 1-norm support vector machines. In S. Thrun, L Saul,
and B. Schölkopf, editors, Advances in Neural Information Processing Systems 16, pages 49–56,
Cambridge, MA, 2004. MIT Press.
H. Zou. The adaptive Lasso and its oracle properties. Journal of the American Statistical Association, 101(476):1418–1429, 2006.
H. Zou, T. Hastie, and R. Tibshirani. Sparse principal component analysis. Journal of Computational and Graphical Statistics, 15(2):262–286, 2006.
H. Zou, T. Hastie, and R. Tibshirani. On the degrees of freedom of the Lasso, 2007. To appear in
Annals of Statistics.

1555