Journal of Machine Learning Research 12 (2011) 2121-2159

Submitted 3/10; Revised 3/11; Published 7/11

Adaptive Subgradient Methods forOnline Learning and Stochastic Optimization∗
John Duchi

JDUCHI @ CS . BERKELEY. EDU

Computer Science Division
University of California, Berkeley
Berkeley, CA 94720 USA

Elad Hazan

EHAZAN @ IE . TECHNION . AC . IL

Technion - Israel Institute of Technology
Technion City
Haifa, 32000, Israel

Yoram Singer

SINGER @ GOOGLE . COM

Google1600 Amphitheatre Parkway
Mountain View, CA 94043 USA

Editor: Tong Zhang

Abstract
We present a new family of subgradient methods that dynamically incorporate knowledge of thegeometry of the data observed in earlier iterations to perform more informative gradient-basedlearning. Metaphorically, the adaptation allows us to find needles in haystacks in the form of very
predictive but rarely seen features. Our paradigm stems from recent advances in stochastic optimization and online learning which employ proximal functions to control the gradient steps ofthe algorithm. We describe and analyze an apparatus for adaptively modifying the proximal function, which significantly simplifies setting a learning rate and results in regret guarantees that areprovably as good as the best proximal function that can be chosen in hindsight. We give several
efficient algorithms for empirical risk minimization problems with common and important regularization functions and domain constraints. We experimentally study our theoretical analysis andshow that adaptive subgradient methods outperform state-of-the-art, yet non-adaptive, subgradientalgorithms.
Keywords: subgradient methods, adaptivity, online learning, stochastic convex optimization

1. Introduction
In many applications of online and stochastic learning, the input instances are of very high dimension, yet within any particular instance only a few features are non-zero. It is often the case,however, that infrequently occurring features are highly informative and discriminative. The informativeness of rare features has led practitioners to craft domain-specific feature weightings, such as
TF-IDF (Salton and Buckley, 1988), which pre-emphasize infrequently occurring features. We use
this old idea as a motivation for applying modern learning-theoretic techniques to the problem of
online and stochastic learning, focusing concretely on (sub)gradient methods.
∗. A preliminary version of this work was published in COLT 2010.
c 2011 John Duchi, Elad Hazan and Yoram Singer.

D UCHI , H AZAN AND S INGER

Standard stochastic subgradient methods largely follow a predetermined procedural scheme thatis oblivious to the characteristics of the data being observed. In contrast, our algorithms dynamicallyincorporate knowledge of the geometry of the data observed in earlier iterations to perform more
informative gradient-based learning. Informally, our procedures give frequently occurring features
very low learning rates and infrequent features high learning rates, where the intuition is that eachtime an infrequent feature is seen, the learner should “take notice.” Thus, the adaptation facilitatesfinding and identifying very predictive but comparatively rare features.
1.1 The Adaptive Gradient Algorithm
Before introducing our adaptive gradient algorithm, which we term A DAG RAD, we establish notation. Vectors and scalars are lower case italic letters, such as x ∈ X . We denote a sequence of
vectors by subscripts, that is, xt , xt+1 , . . ., and entries of each vector by an additional subscript, forexample, xt, j . The subdifferential set of a function f evaluated at x is denoted ∂ f (x), and a particular vector in the subdifferential set is denoted by f ′ (x) ∈ ∂ f (x) or gt ∈ ∂ ft (xt ). When a function
is differentiable, we write ∇ f (x). We use hx, yi to denote the inner product between x and y. The
Bregman divergence associated with a strongly convex and differentiable function ψ isBψ (x, y) = ψ(x) − ψ(y) − h∇ψ(y), x − yi .
We also make frequent use of the following two matrices. Let g1:t = [g1 · · · gt ] denote the matrix
obtained by concatenating the subgradient sequence. We denote the ith row of this matrix, whichamounts to the concatenation of the ith component of each subgradient we observe, by g1:t,i . We
also define the outer product matrix Gt = ∑tτ=1 gτ gτ⊤ .
Online learning and stochastic optimization are closely related and basically interchangeable(Cesa-Bianchi et al., 2004). In order to keep our presentation simple, we confine our discussion andalgorithmic descriptions to the online setting with the regret bound model. In online learning, thelearner repeatedly predicts a point xt ∈ X ⊆ Rd , which often represents a weight vector assigning
importance values to various features. The learner’s goal is to achieve low regret with respect to astatic predictor x∗ in the (closed) convex set X ⊆ Rd (possibly X = Rd ) on a sequence of functions
ft (x), measured as
T

T

R(T ) = ∑ ft (xt ) − inf ∑ ft (x) .
x∈X

t=1

t=1

At every timestep t, the learner receives the (sub)gradient information gt ∈ ∂ ft (xt ). Standard subgradient algorithms then move the predictor xt in the opposite direction of gt while maintaining
xt+1 ∈ X via the projected gradient update (e.g., Zinkevich, 2003)
xt+1 = ΠX (xt − ηgt ) = argmin kx − (xt − ηgt )k22 .
x∈Xp
In contrast, let the Mahalanobis norm k·kA = h·, A·i and denote the projection of a point y onto X
according to A by ΠAX (y) = argminx∈X kx − ykA = argminx∈X hx − y, A(x − y)i. Using this notation,
our generalization of standard gradient descent employs the update
1/2

G

xt+1 = ΠX t



−1/2
xt − ηGt
gt .
2122

A DAPTIVE S UBGRADIENT M ETHODS

The above algorithm is computationally impractical in high dimensions since it requires computation of the root of the matrix Gt , the outer product matrix. Thus we specialize the update todiag(Gt )1/2
xt+1 = ΠX(1)
xt − η diag(Gt )−1/2 gt .

Both the inverse and root of diag(Gt ) can be computed in linear time. Moreover, as we discuss later,when the gradient vectors are sparse the update above can often be performed in time proportional
to the support of the gradient. We now elaborate and give a more formal discussion of our setting.
In this paper we consider several different online learning algorithms and their stochastic convex
optimization counterparts. Formally, we consider online learning with a sequence of composite
functions φt . Each function is of the form φt (x) = ft (x) + ϕ(x) where ft and ϕ are (closed) convex
functions. In the learning settings we study, ft is either an instantaneous loss or a stochastic estimateof the objective function in an optimization task. The function ϕ serves as a fixed regularization
function and is typically used to control the complexity of x. At each round the algorithm makes aprediction xt ∈ X and then receives the function ft . We define the regret with respect to the fixed(optimal) predictor x∗ as
T

T

t=1

t=1

Rφ (T ) , ∑ [φt (xt ) − φt (x∗ )] = ∑ [ ft (xt ) + ϕ(xt ) − ft (x∗ ) − ϕ(x∗ )] .

(2)

Our goal is to devise algorithms which are guaranteed to suffer asymptotically sub-linear regret,namely, Rφ (T ) = o(T ).
Our analysis applies to related, yet different, methods for minimizing the regret (2). The firstis Nesterov’s primal-dual subgradient method (2009), and in particular Xiao’s (2010) extension,
regularized dual averaging, and the follow-the-regularized-leader (FTRL) family of algorithms (seefor instance Kalai and Vempala, 2003; Hazan et al., 2006). In the primal-dual subgradient method
the algorithm makes a prediction xt on round t using the average gradient ḡt = 1t ∑tτ=1 gτ . The update
encompasses a trade-off between a gradient-dependent linear term, the regularizer ϕ, and a stronglyconvex term ψt for well-conditioned predictions. Here ψt is the proximal term. The update amounts
to solving1
xt+1 = argmin η hḡt , xi + ηϕ(x) + ψt (x) ,
(3)
t
x∈X
where η is a fixed step-size and x1 = argminx∈X ϕ(x). The second method similarly has numerous names, including proximal gradient, forward-backward splitting, and composite mirror descent(Tseng, 2008; Duchi et al., 2010). We use the term composite mirror descent. The composite mirror
descent method employs a more immediate trade-off between the current gradient gt , ϕ, and staying
close to xt using the proximal function ψ,
xt+1 = argmin η hgt , xi + ηϕ(x) + Bψt (x, xt ) .
(4)x∈X

Our work focuses on temporal adaptation √of the proximal function in a data driven way, while
previous work simply sets ψt ≡ ψ, ψt (·) = tψ(·), or ψt (·) = tψ(·) for some fixed ψ.
We provide formal analyses equally applicable to the above two updates and show how to automatically choose the function ψt so as to achieve asymptotically small regret. We describe andanalyze two algorithms. Both algorithms use squared Mahalanobis norms as their proximal functions, setting ψt (x) = hx, Ht xi for a symmetric matrix Ht  0. The first uses diagonal matrices while2123

D UCHI , H AZAN AND S INGER

the second constructs full dimensional matrices. Concretely, for some small fixed δ ≥ 0 (specified
later, though in practice δ can be set to 0) we set
1/2

Ht = δI + diag(Gt )1/2 (Diagonal) and Ht = δI + Gt

(Full) .

(5)

Plugging the appropriate matrix from the above equation into ψt in (3) or (4) gives rise to our
A DAG RAD family of algorithms. Informally, we obtain algorithms which are similar to secondorder gradient descent by constructing approximations to the Hessian of the functions ft , though we
use roots of the matrices.
1.2 Outline of Results
We now outline our results, deferring formal statements of the theorems to later sections. Recall thedefinitions of g1:t as the matrix of concatenated subgradients and Gt as the outer product matrix inthe prequel. The A DAG RAD algorithm with full matrix divergences entertains bounds of the form




1/2
1/2
∗
∗
Rφ (T ) = O kx k2 tr(GT ) and Rφ (T ) = O max kxt − x k2 tr(GT ) .
t≤T

We further show that



1/2

tr GT



v (
)
u
T
u
1/2 t
hgt , S−1 gt i : S  0, tr(S) ≤ d .
=d
inf
S

∑

t=1

These results are formally given in Theorem 7 and its corollaries. When our proximal function
ψt (x) = x, diag(Gt )1/2 x we have bounds attainable in time at most linear in the dimension d ofour problems of the form!
!
dRφ (T ) = O kx∗ k∞ ∑ kg1:T,i k2
i=1d

and Rφ (T ) = O max kxt − x∗ k∞ ∑ kg1:T,i k2 .
t≤T

i=1

Similar to the above, we will show that
v (
)
u
d
T
u
1/2 t
kg1:T,i k = d
hgt , diag(s)−1 gt i : s  0, h1, si ≤ d .
inf

∑

2

s

i=1

∑

t=1

We formally state the above two regret bounds in Theorem 5 and its corollaries.
Following are a simple example and corollary to Theorem 5 to illustrate one regime in which
we expect substantial improvements (see also the next subsection). Let ϕ ≡ 0 and consider Zinkevich’s online gradient descent algorithm. Given a compact convex set X ⊆ Rd and sequenceof convex functions√ft , Zinkevich’s algorithm makes the sequence of predictions x1 , . . . , xT with
xt+1 = ΠX (xt − (η/ t)gt ). If the diameter of X is bounded, thus supx,y∈X kx − yk2 ≤ D2 , then online gradient descent, with the optimal choice in hindsight for the stepsize η (see the bound (7) in
Section 1.4), achieves a regret bound of
s
T
T
T
√
(6)
∑ ft (xt ) − inf ∑ ft (x) ≤ 2D2 ∑ kgt k22 .
t=1

x∈X t=1t=1When X is bounded via supx,y∈X kx − yk∞ ≤ D∞ , the following corollary is a simple consequence of
our Theorem 5.
2124

A DAPTIVE S UBGRADIENT M ETHODSCorollary 1 Let the sequence {xt } ⊂ Rd be√generated by the update (4) and assume thatmaxt kx∗ − xt k∞ ≤ D∞ . Using stepsize η = D∞ / 2, for any x∗ , the following bound holds.
s
T
d
√
√
Rφ (T ) ≤ 2dD∞
inf
∑ kgt k2diag(s)−1 = 2D∞ ∑ kg1:T,i k2 .
s0,h1,si≤d t=1i=1

The important feature of the bound above is the infimum under the square root, which allows us toperform better than simply using the identity matrix, and the√fact that the stepsize is easy to set a
priori. For example, if the set X = {x : kxk∞ ≤ 1}, then D2 = 2 d while D∞ = 2, which suggests that
if we are learning a dense predictor over a box, the adaptive method should perform well. Indeed,in this case we are guaranteed that the bound in Corollary 1 is better than (6) as the identity matrix
belongs to the set over which we take the infimum.
To conclude the outline of results, we would like to point to two relevant research papers. First,
Zinkevich’s regret bound is tight and cannot be improved in a minimax sense (Abernethy et al.,2008). Therefore, improving the regret bound requires further reasonable assumptions on the input
space. Second, in a independent work, performed concurrently to the research presented in thispaper, McMahan and Streeter (2010) study competitive ratios, showing guaranteed improvementsof the above bounds relative to families of online algorithms.
1.3 Improvements and Motivating Example
As mentioned in the prequel, we expect our adaptive methods to outperform standard online learning
methods when the gradient vectors are sparse. We give empirical evidence supporting the improved
performance of the adaptive methods in Section 6. Here we give a few abstract examples that showthat for sparse data (input sequences where gt has many zeros) the adaptive methods herein havebetter performance than non-adaptive methods. In our examples we use the hinge loss, that is,
ft (x) = [1 − yt hzt , xi]+ ,where yt is the label of example t and zt ∈ Rd is the data vector.
For our first example, which was also given by McMahan and Streeter (2010), consider thefollowing sparse random data scenario, where the vectors zt ∈ {−1, 0, 1}d . Assume that at in eachround t, feature i appears with probability pi = min{1, ci−α } for some α ∈ (1, ∞) and a dimensionindependent constant c. Then taking the expectation of the gradient terms in the bound in Corollary 1, we have
q

d
d q
d p
d
|{t : |gt,i | = 1}| ≤ ∑ E|{t : |gt,i | = 1}| = ∑ pi T
E ∑ kg1:T,i k2 = ∑ E
i=1

i=1

i=1

i=1

by Jensen’s inequality. In the rightmost sum, we have c ∑di=1 i−α/2 = O(log d) for α ≥ 2, and= {x : kxk∞ ≤ 1}, then∑di=1 i−α/2 = O(d 1−α/2 ) for α ∈ (1, 2). If the domain X is a hypercube, say X √in Corollary 1 D∞ = 2, and the regret of A DAG RAD is O(max{log√d, d 1−α/2 } T ). For contrast, thestandard regret√bound (6) for online gradient descent has D2 = 2 d and kgt k22 ≥ 1, yielding bestcase regret O( dT ). So we see that in this sparse yet heavy tailed feature setting, A DAG RAD’s regret guarantee can be exponentially smaller in the dimension d than the non-adaptive regret bound.
Our remaining examples construct a sparse sequence for which there is a perfect predictor thatthe adaptive methods learn after d iterations, while standard online gradient descent (Zinkevich,
2125

D UCHI , H AZAN AND S INGER

2003) suffers significantly higher√ X is compact, so that for online√ loss. We assume the domain
gradient descent we set ηt = η/ t, which gives the optimal O( T ) regret (the setting of η does notmatter to the adversary we construct).
1.3.1 D IAGONAL A DAPTATION
Consider the diagonal version of our proposed update (4) with X = {x : kxk∞√≤ 1}. Evidently,
we can take D∞ = 2, and this choice simply results in the update xt+1 = xt − 2 diag(Gt )−1/2 gt
followed by projection (1) onto X for A DAG RAD (we use a pseudo-inverse if the inverse does notexist). Let ei denote the ith unit basis vector, and assume that for each t, zt = ±ei for some i. Also
let yt = sign(h1, zt i) so that there exists a perfect classifier x∗ = 1 ∈ X ⊂ Rd . We initialize x1 to be
the zero vector. Fix some ε > 0, and on rounds rounds t = 1, . . . , η2 /ε2 , set zt = e1 . After theserounds, simply choose zt = ±ei for index i ∈ {2, . . . , d} chosen at random. It is clear that the update
to parameter xi at these iterations is different, and amounts toη
xt+1 = xt + ei A DAG RAD
xt+1 = xt + √(Gradient Descent) .
t [−1,1]d(Here [·][−1,1]d denotes the truncation of the vector to [−1, 1]d ). In particular, after suffering d − 1more losses, A√
DAG RAD has a perfect classifier. However, on the remaining iterations gradientdescent has η/ t ≤ ε and thus evidently suffers loss at least d/(2ε). Of course, for small ε, wehave d/(2ε) ≫ d. In short, A DAG RAD achieves constant regret per dimension while online gradientdescent can suffer arbitrary loss (for unbounded t). It seems quite silly, then, to use a global learning
rate rather than one for each feature.
Full Matrix Adaptation. We use a similar construction to the diagonal case to show a situationin which the full matrix update from (5) gives substantially
gradient√ lower regret than stochasticd×d
kxkd}. Let V = [v1 . . . vd ] ∈ R
be andescent. For full divergences we set X = {x :2 ≤orthonormal matrix. Instead of having zt cycle through the unit vectors, we make zt cycle throughthe vi so that zt = ±vi . We let the label yt = sign( 1,V ⊤ zt ) = sign ∑di=1 hvi , zt i . We provide anelaborated explanation in Appendix A. Intuitively, with ψt (x) = hx, Ht xi and Ht set to be the fullmatrix from (5), A DAG RAD again needs to observe each orthonormal vector vi only once while
stochastic gradient descent’s loss can be made Ω(d/ε) for any ε > 0.
1.4 Related Work
Many successful algorithms have been developed over the past few years to minimize regret in
the online learning setting. A modern view of these algorithms casts the problem as the task offollowing the (regularized) leader (see Rakhlin, 2009, and the references therein) or FTRL in short.
Informally, FTRL methods choose the best decision in hindsight at every iteration. Verbatim usageof the FTRL approach fails to achieve low regret, however, adding a proximal1 term to the past
predictions leads to numerous low regret algorithms (Kalai and Vempala, 2003; Hazan and Kale,
2008; Rakhlin, 2009). The proximal term strongly affects the performance of the learning algorithm.
Therefore, adapting the proximal function to the characteristics of the problem at hand is desirable.
Our approach is thus motivated by two goals. The first is to generalize the agnostic online learning paradigm to the meta-task of specializing an algorithm to fit a particular data set. Specifically,
1. The proximal term is also referred to as regularization in the online learning literature. We use the phrase proximal
term in order to avoid confusion with the statistical regularization function ϕ.

2126

A DAPTIVE S UBGRADIENT M ETHODS

we change the proximal function to achieve performance guarantees which are competitive with thebest proximal term found in hindsight. The second, as alluded to earlier, is to automatically adjust
the learning rates for online learning and stochastic gradient descent on a per-feature basis. The
latter can be very useful when our gradient vectors gt are sparse, for example, in a classification
setting where examples may have only a small number of non-zero features. As we demonstrated
in the examples above, it is rather deficient to employ exactly the same learning rate for a feature
seen hundreds of times and for a feature seen only once or twice.
Our techniques stem from a variety of research directions, and as a byproduct we also extend a
few well-known algorithms. In particular, we consider variants of the follow-the-regularized leader(FTRL) algorithms mentioned above, which are kin to Zinkevich’s lazy projection algorithm. We
use Xiao’s recently analyzed regularized dual averaging (RDA) algorithm (2010), which builds upon
Nesterov’s (2009) primal-dual subgradient method. We also consider forward-backward splitting(F OBOS) (Duchi and Singer, 2009) and its composite mirror-descent (proximal gradient) generalizations (Tseng, 2008; Duchi et al., 2010), which in turn include as special cases projected gradients(Zinkevich, 2003) and mirror descent (Nemirovski and Yudin, 1983; Beck and Teboulle, 2003). Recent work by several authors (Nemirovski et al., 2009; Juditsky et al., 2008; Lan, 2010; Xiao, 2010)
considered efficient and robust methods for stochastic optimization, especially in the case when the
expected objective f is smooth. It may be interesting to investigate adaptive metric approaches insmooth stochastic optimization.
The idea of adapting first order optimization methods is by no means new and can be traced
back at least to the 1970s with the work on space dilation methods of Shor (1972) and variable
metric methods, such as the BFGS family of algorithms (e.g., Fletcher, 1970). This prior workoften assumed that the function to be minimized was differentiable and, to our knowledge, did notconsider stochastic, online, or composite optimization. In her thesis, Nedić (2002) studied variable
metric subgradient methods, though it seems difficult to derive explicit rates of convergence from theresults there, and the algorithms apply only when the constraint set X = Rd . More recently, Bordes
et al. (2009) proposed a Quasi-Newton stochastic gradient-descent procedure, which is similar in
spirit to our methods. However, their convergence results assume a smooth objective with positive
definite Hessian bounded away from 0. Our results apply more generally.
Prior to the analysis presented in this paper for online and stochastic optimization, the strongly
convex function ψ in the update equations (3) and (4) either remained intact or was simply multiplied
by a time-dependent scalar throughout the run of the algorithm. Zinkevich’s√ projected gradient,
for example, uses ψt (x) = kxk22 , while RDA (Xiao, 2010) employs ψt (x) = tψ(x) where ψ is astrongly convex function. The bounds for both types of algorithms are similar, and both rely on thenorm k·k (and its associated dual k·k∗ ) with respect to which ψ is strongly convex. Mirror-descenttype first order algorithms, such as projected gradient methods, attain regret bounds of the form
(Zinkevich, 2003; Bartlett et al., 2007; Duchi et al., 2010)

Rφ (T ) ≤η T
12Bψ (x∗ , x1 ) + ∑ ft′ (xt ) ∗ .
η
2 t=1

(7)

√
√
Choosing η ∝ 1/ √T gives Rφ (T ) = O( T ). When Bψ (x, x∗√) is bounded for all x ∈ X , we choose
step sizes ηt ∝ 1/ t which is equivalent to setting ψt (x) = tψ(x). Therefore, no assumption on
the time horizon is necessary. For RDA and follow-the-leader algorithms, the bounds are similar
2127

D UCHI , H AZAN AND S INGER(Xiao, 2010, Theorem 3):
Rφ (T ) ≤√1 T
2T ψ(x∗ ) + √ ∑ ft′ (xt ) ∗ .
2 T t=1(8)

The problem of adapting to data and obtaining tighter data-dependent bounds for algorithms
such as those above is a natural one and has been studied in the mistake-bound setting for online
learning in the past. A framework that is somewhat related to ours is the confidence weighted
learning scheme by Crammer et al. (2008) and the adaptive regularization of weights algorithm(AROW) of Crammer et al. (2009). These papers provide mistake-bound analyses for secondorder algorithms, which in turn are similar in spirit to the second-order Perceptron algorithm (CesaBianchi et al., 2005). The analyses by Crammer and colleagues, however, yield mistake bounds
dependent on the runs of the individual algorithms and are thus difficult to compare with our regretbounds.
AROW maintains a mean prediction vector µt ∈ Rd and a covariance matrix Σt ∈ Rd×d over µt
as well. At every step of the algorithm, the learner receives a pair (zt , yt ) where zt ∈ Rd is the tth
example and yt ∈ {−1, +1} is the label. Whenever the predictor µt attains a margin value smaller
than 1, AROW performs the update
βt =1, αt = [1 − yt hzt , µt i]+ ,hzt , Σt zt i + λµt+1 = µt + αt Σt yt zt , Σt+1 = Σt − βt Σt xt xt⊤ Σt .

(9)In the above scheme, one can force Σt to be diagonal, which reduces the run-time and storage
requirements of the algorithm but still gives good performance (Crammer et al., 2009). In contrast
to AROW, the A DAG RAD algorithm uses the root of the inverse covariance matrix, a consequence of
our formal analysis. Crammer et al.’s algorithm and our algorithms have similar run times, generally
linear in the dimension d, when using diagonal matrices. However, when using full matrices theruntime of AROW algorithm is O(d 2 ), which is faster than ours as it requires computing the root ofa matrix.
In concurrent work, McMahan and Streeter (2010) propose and analyze an algorithm which
is very similar to some of the algorithms presented in this paper. Our analysis builds on recent
advances in online learning and stochastic optimization (Duchi et al., 2010; Xiao, 2010), whereas
McMahan and Streeter use first-principles to derive their regret bounds. As a consequence of ourapproach, we are able to apply our analysis to algorithms for composite minimization with a known
additional objective term ϕ. We are also able to generalize and analyze both the mirror descent anddual-averaging family of algorithms. McMahan and Streeter focus on what they term the competitive ratio, which is the ratio of the worst case regret of the adaptive algorithm to the worst caseregret of a non-adaptive algorithm with the best proximal term ψ chosen in hindsight. We touch on
this issue briefly in the sequel, but refer the interested reader to McMahan and Streeter (2010) for
this alternative elegant perspective. We believe that both analyses shed insights into the problems
studied in this paper and complement each other.
There are also other lines of work on adaptive gradient methods that are not directly related to
our work but nonetheless relevant. Tighter regret bounds using the variation of the cost functions ft
were proposed by Cesa-Bianchi et al. (2007) and derived by Hazan and Kale (2008). Bartlett et al.
(2007) explore another adaptation technique for ηt where they adapt the step size to accommodate
2128

A DAPTIVE S UBGRADIENT M ETHODS

both strongly and weakly convex functions. Our approach differs from previous approaches as it
does not focus on a particular loss function or mistake bound. Instead, we view the problem ofadapting the proximal function as a meta-learning problem. We then obtain a bound comparable tothe bound obtained using the best proximal function chosen in hindsight.

2. Adaptive Proximal Functions
Examining the bounds (7) and (8), we see that most of the regret depends on dual norms of ft′ (xt ),and the dual norms in turn depend on the choice of ψ. This naturally leads to the question of whether
we can modify the proximal term ψ along the run of the algorithm in order to lower the contributionof the aforementioned norms. We achieve this goal by keeping second order information about thesequence ft and allow ψ to vary on each round of the algorithms.
We begin by providing two corollaries based on previous work that give the regret of our basealgorithms when the proximal function ψt is allowed to change. These corollaries are used in
the sequel in our regret analysis. We assume that ψt is monotonically non-decreasing, that is,ψt+1 (x) ≥ ψt (x). We also assume that ψt is 1-strongly convex with respect to a time-dependentsemi-norm k·kψt . Formally, ψ is 1-strongly convex with respect to k·kψ if
1
ψ(y) ≥ ψ(x) + h∇ψ(x), y − xi + kx − yk2ψ .
2Strong convexity is guaranteed if and only if Bψt (x, y) ≥ 21 kx − yk2ψt . We also denote the dual normof k·kψt by k·kψt∗ . For completeness, we provide the proofs of following two results in Appendix F,
as they build straightforwardly on work by Duchi et al. (2010) and Xiao (2010). For the primal-dual
subgradient update, the following bound holds.
Proposition 2 Let the sequence {xt } be defined by the update (3). For any x∗ ∈ X ,T∑ ft (xt ) + ϕ(xt ) − ft (x∗ ) − ϕ(x∗ ) ≤

t=1

η T
1
2
ψT (x∗ ) + ∑ ft′ (xt ) ψ∗ .
t−1
η
2 t=1

(10)

For composite mirror descent algorithms a similar result holds.
Proposition 3 Let the sequence {xt } be defined by the update (4). Assume w.l.o.g. that ϕ(x1 ) = 0.
For any x∗ ∈ X ,T∑ ft (xt ) + ϕ(xt ) − ft (x∗ ) − ϕ(x∗ )

t=1

≤

 η T
1 T −1 
1
2
Bψ1 (x∗ , x1 ) + ∑ Bψt+1 (x∗ , xt+1 ) − Bψt (x∗ , xt+1 ) + ∑ ft′ (xt ) ψ∗ .
t
η
η t=1
2 t=1

(11)

The above corollaries allow us to prove regret bounds for a family of algorithms that iteratively
modify the proximal functions ψt in attempt to lower the regret bounds.
2129

D UCHI , H AZAN AND S INGERI NPUT: η > 0, δ ≥ 0
VARIABLES: s ∈ Rd , H ∈ Rd×d , g1:t,i ∈ Rt for i ∈ {1, . . . , d}I NITIALIZE x1 = 0, g1:0 = []F OR t = 1 to T
Suffer loss ft (xt )Receive subgradient gt ∈ ∂ ft (xt ) of ft at xt
U PDATE g1:t = [g1:t−1 gt ], st,i = kg1:t,i k2
S ET Ht = δI + diag(st ), ψt (x) = 12 hx, Ht xi
Primal-Dual SubgradientUpdate( *+ (3):)
t1
1
xt+1 = argmin η∑ gτ , x + ηϕ(x) + t ψt (x) .
t τ=1
x∈X
Composite Mirror
 Descent Update (4):xt+1 = argmin ηhgt , xi + ηϕ(x) + Bψt (x, xt ) .
x∈XFigure 1: A DAG RAD with diagonal matrices

3. Diagonal Matrix Proximal Functions
We begin by restricting ourselves to using diagonal matrices to define matrix proximal functionsand (semi)norms. This restriction serves a two-fold purpose. First, the analysis for the general caseis somewhat complicated and thus the analysis of the diagonal restriction serves as a proxy for better
understanding. Second, in problems with high dimension where we expect this type of modification
to help, maintaining more complicated proximal functions is likely to be prohibitively expensive.
Whereas earlier analysis requires a learning rate to slow changes between predictors xt and xt+1 , wewill instead automatically grow the proximal function we use to achieve asymptotically low regret.
To remind the reader, g1:t,i is the ith row of the matrix obtained by concatenating the subgradientsfrom iteration 1 through t in the online algorithm.
To provide some intuition for the algorithm we show in Algorithm 1, let us examine the problem
2gt,is.t. s  0, h1, si ≤ c .
t=1 i=1 siTdmin ∑ ∑s

This problem is solved by setting si = kg1:T,i k2 and scaling s so that hs, 1i = c. To see this, we canwrite the Lagrangian of the minimization problem by introducing multipliers λ  0 and θ ≥ 0 to get
kg1:T,i k22− hλ, si + θ(h1, si − c).
si
i=1
d

L (s, λ, θ) = ∑Taking partial derivatives to find the infimum of L , we see that − kg1:T,i k22 /s2i − λi + θ = 0, and complementarity conditions on λi si (Boyd and Vandenberghe, 2004) imply that λi = 0. Thus we have1
si = θ− 2 kg1:T,i k2 , and normalizing appropriately using θ gives that si = c kg1:T,i k2 / ∑dj=1 g1:T, j 2 .
2130

A DAPTIVE S UBGRADIENT M ETHODSAs a final note, we can plug si into the objective above to see()
T d g2
1t,i
inf ∑ ∑: s  0, h1, si ≤ c =sc
t=1 i=1 sid∑ kg1:T,i k2

i=1

!2

.

(12)

Let diag(v) denote the diagonal matrix with diagonal v. It is natural to suspect that for s achieving
the infimum in Equation (12), if we use a proximal function similar to ψ(x) = hx, diag(s)xi withassociated squared dual norm kxk2ψ∗ = x, diag(s)−1 x , we should do well lowering the gradientterms in the regret bounds (10) and (11).
To prove a regret bound for our Algorithm 1, we note that both types of updates suffer losses that
include a term depending solely on the gradients obtained along their run. The following lemma
is applicable to both updates, and was originally proved by Auer and Gentile (2000), though we
provide a proof in Appendix C. McMahan and Streeter (2010) also give an identical lemma.
Lemma 4 Let gt = ft′ (xt ) and g1:t and st be defined as in Algorithm 1. Then
T

d

t=1

i=1

∑ gt , diag(st )−1 gt ≤ 2 ∑ kg1:T,i k2 .

To obtain a regret bound, we need to consider the terms consisting of the dual-norm of the subgradient in the regret bounds (10) and (11), which is k ft′ (xt )k2ψt∗ . When ψt (x) = hx, (δI + diag(st ))xi,
it is easy to see that the associated dual-norm is
kgk2ψt∗ = g, (δI + diag(st ))−1 g .
From the definition of st in Algorithm 1, we clearly have k ft′ (xt )k2ψt∗ ≤ gt , diag(st )−1 gt . Note that
if st,i = 0 then gt,i = 0 by definition of st,i . Thus, for any δ ≥ 0, Lemma 4 implies
T

∑

t=1

d

ft′ (xt ) ψ∗ ≤ 2 ∑ kg1:T,i k2 .
2

t

(13)

i=1

To obtain a bound for a primal-dual subgradient method, we set δ ≥ maxt kgt k∞ , in which casekgt k2ψ∗ ≤ gt , diag(st )−1 gt , and we follow the same lines of reasoning to achieve the inequalt−1
ity (13).
It remains to bound the various Bregman divergence terms for Corollary 3 and the term ψT (x∗ )for Corollary 2. We focus first on the composite mirror-descent update. Examining the bound (11)and Algorithm 1, we notice that1 ∗hx − xt+1 , diag(st+1 − st )(x∗ − xt+1 )i2
1≤ max(xi∗ − xt+1,i )2 kst+1 − st k1 .
2 i

Bψt+1 (x∗ , xt+1 ) − Bψt (x∗ , xt+1 ) =Since kst+1 − st k1 = hst+1 − st , 1i and hsT , 1i = ∑di=1 kg1:T,i k2 , we have
1 T −1T −1∑ Bψ (x∗ , xt+1 ) − Bψ (x∗ , xt+1 ) ≤ 2 ∑ kx∗ − xt+1 k2∞ hst+1 − st , 1i
t+1

t

t=1

t=1

d

≤

1
1
max kx∗ − xt k2∞ ∑ kg1:T,i k2 − kx∗ − x1 k2∞ hs1 , 1i .
2 t≤T
2
i=1
2131

(14)

D UCHI , H AZAN AND S INGER

We also haved

ψT (x∗ ) = δ kx∗ k22 + hx∗ , diag(sT )x∗ i ≤ δ kx∗ k22 + kx∗ k2∞ ∑ kg1:T,i k2 .
i=1

Combining the above arguments with Corollaries 2 and 3, and using (14) with the fact that Bψ1 (x∗ , x1 ) ≤2
1∗
2 kx − x1 k∞ h1, s1 i, we have proved the following theorem.
Theorem 5 Let the sequence {xt } be defined by Algorithm 1. For xt generated using the primaldual subgradient update (3) with δ ≥ maxt kgt k∞ , for any x∗ ∈ X ,Rφ (T ) ≤d
δ ∗ 2 1 ∗ 2 d
kx k2 + kx k∞ ∑ kg1:T,i k2 + η ∑ kg1:T,i k2 .
η
η
i=1
i=1

For xt generated using the composite mirror-descent update (4), for any x∗ ∈ X
Rφ (T ) ≤dd
1
max kx∗ − xt k2∞ ∑ kg1:T,i k2 + η ∑ kg1:T,i k2 .
2η t≤T
i=1
i=1

The above theorem is a bit unwieldy. We thus perform a few algebraic simplifications to get thenext corollary, which has a more intuitive form. Let us assume that X is compact and set D∞ =supx∈X kx − x∗ k∞ . Furthermore, define)(dγT , ∑ kg1:T,i k2 = inf
i=1

s

T

d

t=1

i=1

∑ gt , diag(s)−1 gt : h1, si ≤ ∑ kg1:T,i k2 , s  0.

Also w.l.o.g. let 0 ∈ X . The√ following corollary is immediate (this is equivalent to Corollary 1,
though we have moved the d term in the earlier bound).
Corollary 6 Assume that D∞ and γT are defined as above. For {xt } generated by Algorithm 1 using
the primal-dual subgradient update (3) with η = kx∗ k∞ , for any x∗ ∈ X we haveRφ (T ) ≤ 2 kx∗ k∞ γT + δkx∗ k22≤ 2 kx∗ k∞ γT + δ kx∗ k1 .
kx∗ k∞√
Using the composite mirror descent update (4) to generate {xt } and setting η = D∞ / 2, we haveRφ (T ) ≤d√√
2D∞ ∑ kg1:T,i k2 = 2D∞ γT .
i=1

We now give a short derivation of Corollary 1 from the introduction: use Theorem 5, Corollary 6,
and the fact that!2)(
T d g2
1 dt,i: s  0, h1, si ≤ d =inf ∑ ∑
∑ kg1:T,i k2 .
s
d i=1
t=1 i=1 si
as in (12)√ in the beginning of Section 3. Plugging the γT term in from Corollary 6 and multiplying
D∞ by d completes the proof of the corollary.
2132

A DAPTIVE S UBGRADIENT M ETHODSAs discussed in the introduction, Algorithm 1 should have lower regret than non-adaptive algorithms on sparse data, though this depends on the geometry of the underlying optimization space
X . For example, suppose that our learning problem is a logistic regression with 0/1-valued features.
Then the gradient terms are likewise based on 0/1-valued featuresand sparse, so the gradient terms√
dkg
kin the bound ∑i=1 1:T,i 2 should all be much smaller than T . If some features appear much more
frequently than others, then the infimal representation of γT and the infimal equality in Corollary 1
show that we have significantly lower regret by using higher learning rates for infrequent featuresand lower learning rates on commonly appearing features. Further, if the optimal predictor is relatively dense, as is often the case in predictions problems with sparse inputs, then kx∗ k∞ is the best
p-norm we can have in the regret.
More precisely, McMahan and Streeter (2010) show that if X is contained within an ℓ∞ ballof radius √R and contains an ℓ∞ ball of radius r, then the bound in the above corollary is within afactor of 2R/r of the regret of the best diagonal proximal matrix, chosen in hindsight. So, forexample, if X = {x ∈ Rd : kxk p ≤ C}, then R/r = d 1/p , which shows that the domain X does effect
the guarantees we can give on optimality of A DAG RAD.

4. Full Matrix Proximal Functions
In this section we derive and analyze new updates when we estimate a full matrix for the divergenceψt instead of a diagonal one. In this generalized case, we use the root of the matrix of outer productsof the gradients that we have observed to update our parameters. As in the diagonal case, we build
on intuition garnered from an optimization problem, and in particular, we seek a matrix S which isthe solution to the following minimization problem:Tmin ∑ gt , S−1 gt
S

t=1s.t. S  0, tr(S) ≤ c .

(15)

The solution is obtained by defining Gt = ∑tτ=1 gτ gτ⊤ and setting S to be a normalized version of
1/2
1/2
the root of GT , that is, S = c GT / tr(GT ). For a proof, see Lemma 15 in Appendix E, which also
shows that when GT is not full rank
D we canE instead use its pseudo-inverse. If we iteratively use
1/2
divergences of the form ψt (x) = x, Gt x , we might expect as in the diagonal case to attain low
regret by collecting gradient information. We achieve our low regret goal by employing a similar
doubling lemma to Lemma 4 and bounding the gradient norm terms. The resulting algorithm isgiven in Algorithm 2, and the next theorem provides a quantitative analysis of the brief motivation
above.
Theorem 7 Let Gt be the outer product matrix defined above and the sequence {xt } be defined by
Algorithm 2. For xt generated using the primal-dual subgradient update of (3) and δ ≥ maxt kgt k2 ,for any x∗ ∈ X
δ
1
1/2
1/2
Rφ (T ) ≤ kx∗ k22 + kx∗ k22 tr(GT ) + η tr(GT ).
ηη
For xt generated with the composite mirror-descent update of (4), if x∗ ∈ X and δ ≥ 0
Rφ (T ) ≤

1
δ ∗ 2
1/2
1/2
kx k2 +
max kx∗ − xt k22 tr(GT ) + η tr(GT ).
η2η t≤T
2133

D UCHI , H AZAN AND S INGERI NPUT: η > 0, δ ≥ 0
VARIABLES: St ∈ Rd×d , Ht ∈ Rd×d , Gt ∈ Rd×d
I NITIALIZE x1 = 0, S0 = 0, H0 = 0, G0 = 0
F OR t = 1 to T
Suffer loss ft (xt )Receive subgradient gt ∈ ∂ ft (xt ) of ft at xt
1

U PDATE Gt = Gt−1 + gt gt⊤ , St = Gt2
S ET Ht = δI + St , ψt (x) = 21 hx, Ht xi

Primal-Dual SubgradientUpdate( *+ ((3)):)
t1
1
xt+1 = argmin η∑ gτ , x + ηϕ(x) + t ψt (x) .
t τ=1
x∈X
Composite Mirror
 Descent Update ((4)):
xt+1 = argmin ηhgt , xi + ηϕ(x) + Bψt (x, xt ) .
x∈XFigure 2: A DAG RAD with full matrices
Proof To begin, we consider the difference between the divergence terms at time t + 1 and time t
from the regret (11) in Corollary 3. Let λmax (M) denote the largest eigenvalue of a matrix M. We
have
E
1D ∗
∗
∗
1/2
1/2
∗
Bψt+1 (x , xt+1 ) − Bψt (x , xt+1 ) =x − xt+1 , (Gt+1 − Gt )(x − xt+1 )
2
1 ∗
1
1/2
1/2
1/2
1/2
kx − xt+1 k22 λmax (Gt+1 − Gt ) ≤ kx∗ − xt+1 k22 tr(Gt+1 − Gt ) .
≤2
2
For the last inequality we used the fact that the trace of a matrix is equal to the sum of its eigenvalues
1/2
along with the property Gt+1 1/2 − Gt 1/2  0 (see Lemma 13 in Appendix B) and therefore tr(Gt+1 −
1/2

1/2

1/2

Gt ) ≥ λmax (Gt+1 − Gt ). Thus, we get
T −1∑ Bψt+1 (x∗ , xt+1 ) − Bψt (x∗ , xt+1 ) ≤

t=1



1 T −1 ∗
1/2
1/2
2
kx
k
tr(G
−
x
)
−
tr(G
)
.
t+1 2
t
∑
t+1
2 t=1

Now we use the fact that G1 is a rank 1 PSD matrix with non-negative trace to see that
T −1



1/2
1/2
2
∗
kx
k
−
x
tr(G
)
−
tr(G
)
t+1 2
t
∑
t+1

t=1

1/2

≤ max kx∗ − xt k22 tr(GT 1/2 ) − kx∗ − x1 k22 tr(G1 ) .
t≤T

(16)

It remains to bound the gradient terms common to all our bounds. We use the following threelemmas, which essentially directly applicable. We prove the first two in Appendix D.
Lemma 8 Let B  0 and B−1/2 denote the root of the inverse of B when B ≻ 0 and the root of thepseudo-inverse of B otherwise. For any ν such that B − νgg⊤  0 the following inequality holds.
2 tr((B − νgg⊤ )1/2 ) ≤ 2 tr(B1/2 ) − ν tr(B−1/2 gg⊤ ) .
2134

A DAPTIVE S UBGRADIENT M ETHODSD
1/2 E
Lemma 9 Let δ ≥ kgk2 and A  0, then g, (δI + A1/2 )−1 g ≤ g, (A + gg⊤ )†g .

Lemma 10 Let St = Gt 1/2 be as defined in Algorithm 2 and A† denote the pseudo-inverse of A.
Then
E
E
T D
T D
∑ gt , St† gt ≤ 2 ∑ gt , ST† gt = 2 tr(GT 1/2 ) .
t=1

t=1

Proof We prove the lemma by induction. The base case is immediate, since we haveD
E hg , g i
1 1= kg1 k2 ≤ 2 kg1 k2 .
g1 , (G†1 )1/2 g1 =kg1 k2
Now, assume the lemma is true for T − 1, so from the inductive assumption we get
E D
E
E
D
T −1 D
†
†
†
g
,
S
+
g
,
S
.
g
g
≤
2
g
,
S
g
T T T
∑ t T −1 t
∑ t t t
T

t=1

t=1

E
D
T −1
gt , ST† −1 gt as
Since ST −1 does not depend on t we can rewrite ∑t=1
ST† −1 ,

tr

T −1

∑

gt gt⊤

t=1!

= tr((G†T −1 )1/2 GT −1 ) ,where the right-most equality follows from the definitions of St and Gt . Therefore, we getT∑ gt , St† gt

t=1

D

E

E
D
≤ 2 tr((G†T −1 )1/2 GT −1 ) + gT , (G†T )1/2 gTED
1/2= 2 tr(GT −1 ) + gT , (G†T )1/2 gT .

Using Lemma 8 with the substitution B = GT , ν = 1, and g = gt lets us exploit the concavity of the
1/2
function tr(A1/2 ) to bound the above sum by 2 tr(GT ).
N
We can now finalize our proof of the theorem. As in the diagonal case, we have that the squared
dual norm (seminorm when δ = 0) associated with ψt is
kxk2ψt∗ = x, (δI + St )−1 x .
E
DThus it is clear that kgt k2ψt∗ ≤ gt , St† gt . For the dual-averaging algorithms, we use Lemma 9 aboveEDshow that kgt k2ψ∗ ≤ gt , St† gt so long as δ ≥ kgt k2 . Lemma 10’s doubling inequality then implies
t−1
that
T

T

∑ ft′ (xt ) ψ ≤ 2 tr(GT ) and ∑ ft′ (xt ) ψ

t=1

2

1/2

∗
t

t=1

2

∗
t−1

1/2

≤ 2 tr(GT )(17)

for the mirror-descent and primal-dual subgradient algorithm, respectively.
1/2
To finish the proof, Note that Bψ1 (x∗ , x1 ) ≤ 21 kx∗ − x1 k22 tr(G1 ) when δ = 0. By combining thisT −1with the first of the bounds (17) and the bound (16) on ∑t=1
Bψt+1 (x∗ , xt+1 ) − Bψt (x∗ , xt+1 ), Corollary 3 gives the theorem’s statement for the mirror-descent family of algorithms. Combining the2135

D UCHI , H AZAN AND S INGER

1/2

T
k ft′ (xt )k2ψ∗ ≤ 2 tr(GT ) and the bound (16) with Corollary 2 gives the desired boundfact that ∑t=1
t−1
on Rφ (T ) for the primal-dual subgradient algorithms, which completes the proof of the theorem.

As before, we can give a corollary that simplifies the bound implied by Theorem 7. The infimalequality in the corollary uses Lemma 15 in Appendix B. The corollary underscores that for learning problems in which there is a rotation U of the space for which the gradient vectors gt havesmall inner products hgt ,Ugt i (essentially a sparse basis for the gt ) then using full-matrix proximal
functions can attain significantly lower regret.
Corollary 11 Assume that ϕ(x1 ) = 0. Then the regret of the sequence {xt } generated by Algorithm 2
when using the primal-dual subgradient update with η = kx∗ k2 is
1/2

Rφ (T ) ≤ 2 kx∗ k2 tr(GT ) + δ kx∗ k2 .

√Let X be compact set so that supx∈X kx − x∗ k2 ≤ D. Taking η = D/ 2 and using the composite
mirror descent update with δ = 0, we haveRφ (T ) ≤√1/2

2D tr(GT ) =

√

v (
)
u
T
u
2dDtinf ∑ gt⊤ S−1 gt : S  0, tr(S) ≤ d .
S

t=1

5. Derived Algorithms
In this section, we derive updates using concrete regularization functions ϕ and settings of thedomain X for the A DAG RAD framework. We focus on showing how to solve Equations (3) and (4)with the diagonal matrix version of the algorithms we have presented. We focus on the diagonalcase for two reasons. First, the updates often take closed-form in this case and carry some intuition.
Second, the diagonal case is feasible to implement in very high dimensions, whereas the full matrix
version is likely to be confined to a few thousand dimensions. We also discuss how to efficiently
compute the updates when the gradient vectors are sparse.
We begin by noting a simple but useful fact. Let Gt denote either the outer product matrix of
1/2
gradients or its diagonal counterpart and let Ht = δI + Gt , as usual. Simple algebraic manipulations yield that each of the updates (3) and (4) in the prequel can be written in the following form
(omitting the stepsize η):1
xt+1 = argmin hu, xi + ϕ(x) + hx, Ht xi .
(18)
2
x∈X
In particular, at time t for the RDA update, we have u = ηt ḡt . For the composite gradient update (4),111
η hgt , xi + hx − xt , Ht (x − xt )i = hηgt − Ht xt , xi + hx, Ht xi + hxt , Ht xt i222
so that u = ηgt − Ht xt . We now derive algorithms for solving the general update (18). Since mostof the derivations are known, we generally provide only the closed-form solutions or algorithms forthe solutions in the remainder of the subsection, deferring detailed derivations to Appendix G for
the interested reader.
2136

A DAPTIVE S UBGRADIENT M ETHODS

5.1 ℓ1 -regularization
We begin by considering how to solve the minimization problems necessary for Algorithm 1 withdiagonal matrix divergences and ϕ(x) = λ kxk1 . We consider the two updates we proposed anddenote the ith diagonal element of the matrix Ht = δI + diag(st ) from Algorithm 1 by Ht,ii = δ +kg1:t,i k2 . For the primal-dual subgradient update, the solution to (3) amounts to the following simple
update for xt+1,i :ηtxt+1,i = sign (−ḡt,i )[|ḡt,i | − λ]+ .
(19)
Ht,ii
Comparing the update (19) to the standard dual averaging update (Xiao, 2010), which is√
xt+1,i = sign (−ḡt,i ) η t [|ḡt,i | − λ]+ ,it is clear that the difference distills to the step size employed for each coordinate. Our generalizationof RDA yields a dedicated step size for each coordinate inversely proportional to the time-based
norm of the coordinate in the sequence of gradients. Due to the normalization by this term the step
size scales linearly with t, so when Ht,ii is small, gradient information on coordinate i is quickly
incorporated.
The composite mirror-descent update (4) has a similar form that essentially amounts to iterative
shrinkage and thresholding, where the shrinkage differs per coordinate:



η
λη
η
xt+1,i = sign xt,i −
gt,igt,i −
xt,i −
.
Ht,ii
Ht,ii
Ht,ii +
We compare the actual performance of the newly derived algorithms to previously studied versions
in the next section.
For both updates it is clear that we can perform “lazy” computation when the gradient vectors
are sparse, a frequently occurring setting when learning for instance from text corpora. Supposethat from time step t0 through t, the ith component of the gradient is 0. Then we can evaluate theabove updates on demand since Ht,ii remains intact. For composite mirror-descent, at time t whenxt,i is needed, we update
ληxt,i = sign(xt0 ,i ) |xt0 ,i | −(t − t0 ).
Ht0 ,ii+
Even simpler just in time evaluation can be performed for the the primal-dual subgradient update.
Here we need to keep an unnormalized version of the average ḡt . Concretely, we keep track of
ut = t ḡt = ∑tτ=1 gτ = ut−1 + gt , then use the update (19):ηt |ut,i |,−λxt,i = sign(−ut,i )
Ht,ii
t
+
where Ht can clearly be updated lazily in a similar fashion.
5.2 ℓ1 -ball Projections
We next consider the setting in which ϕ ≡ 0 and X = {x : kxk1 ≤ c}, for which it is straightforward to adapt efficient solutions to continuous quadratic knapsack problems (Brucker, 1984). We
2137

D UCHI , H AZAN AND S INGERI NPUT: v  0, a  0, c ≥ 0.
I F ∑i vi ≤ c RETURNz∗ = v 
S ORT vi /ai inton µ = vi j /ai j s.t. vi j /ai j ≥ vi j+1 /aoi j+1ρvi

ρS ET ρ := max ρ : ∑ j=1 ai j vi j − aiρ ∑ j=1 a2i j < cρS ET θ =ρ∑ j=1 ai j vi j −cρ

∑ j=1 a2i j

R ETURN z∗ where z∗i = [vi − θai ]+ .
Figure 3: Project v  0 to {z : ha, zi ≤ c, z  0}.
use the matrix Ht = δI + diag(Gt )1/2 from Algorithm 1. We provide a brief derivation sketch andan O(d log d) algorithm in this section. First, we convert the problem (18) into a projection problem onto a scaled ℓ1 -ball. By making the substitutions z = H 1/2 x and A = H −1/2 , it is clear thatproblem (18) is equivalent tominzz + H −1/2 u

2
2

s.t. kAzk1 ≤ c .
−1/2

Now, by appropriate choice of v = −H −1/2 u = −ηtHt
ḡt for the primal-dual update (3) and
1/2−1/2
v = Ht xt − ηHtgt for the mirror-descent update (4), we arrive at the problem
min
z

d
1
kz − vk22 s.t. ∑ ai |zi | ≤ c .
2
i=1

(20)
−1/2

We can clearly recover xt+1 from the solution z∗ to the projection (20) via xt+1 = Htz∗ .
By the symmetry of the objective (20), we can assume without loss of generality that v  0 andconstrain z  0, and a bit of manipulation with the Lagrangian (see Appendix G) for the problem
shows that the solution z∗ has the form
vi − θ∗ ai if vi ≥ θ∗ ai∗
zi =0otherwisefor some θ∗ ≥ 0. The algorithm in Figure 3 constructs the optimal θ and returns z∗ .
5.3 ℓ2 Regularization
We now turn to the case where ϕ(x) = λ kxk2 while X = Rd . This type of regularization is usefulfor zeroing multiple weights in a group, for example in multi-task or multiclass learning (Obozinski
et al., 2007). Recalling the general proximal step (18), we must solve1
min hu, xi + hx, Hxi + λ kxk2 .
x
2

(21)

There is no closed form solution for this problem, but we give an efficient bisection-based procedure
for solving (21). We start by deriving the dual. Introducing a variable z = x, we get the equivalent
problem of minimizing hu, xi + 21 hx, Hxi + λ kzk2 subject to x = z. With Lagrange multipliers α forthe equality constraint, we obtain the Lagrangian12

L (x, z, α) = hu, xi + hx, Hxi + λ kzk2 + hα, x − zi .
2138

A DAPTIVE S UBGRADIENT M ETHODSI NPUT: u ∈ Rd , H  0, λ > 0.
I F kuk2 ≤ λ
R ETURN x = 0
S ET v = H −1 u, θmax = kvk2 /λ − 1/σmin (H)θmin = kvk2 /λ − 1/σmax (H)W HILE θmax − θmin > ε
S ET θ = (θmax + θmin )/2, α(θ) = −(H −1 + θI)−1 v
I F kα(θ)k2 > λS ET θmin = θ
E LSES ET θmax = θ
R ETURN x = −H −1 (u + α(θ))
Figure 4: Minimize hu, xi + 12 hx, Hxi + λ kxk2
Taking the infimum of L with respect to the primal variables x and z, we see that the infimum isattained at x = −H −1 (u + α). Coupled with the fact that infz λ kzk2 − hα, zi = −∞ unless kαk2 ≤ λ,in which case the infimum is 0, we arrive at the dual form
inf L (x, z, α) =x,z

− 21 u + α, H −1 (u + α)−∞

if kαk2 ≤ λotherwise.

Setting v = H −1 u, we further distill the dual tomin hv, αi +α1α, H −1 α2s.t. kαk2 ≤ λ .

(22)

We can solve problem (22) efficiently using a bisection search of its equivalent representation in
Lagrange form,θ
1min hv, αi + α, H −1 α + kαk22 ,α22where θ > 0 is an unknown scalar. The solution to the latter as a function of θ is clearly α(θ) =−(H −1 + θI)−1 v = −(H −1 + θI)−1 H −1 u. Since kα(θ)k2 is monotonically decreasing in θ (consider
the the eigen-decomposition of the positive definite H −1 ), we can simply perform a bisection searchover θ, checking at each point whether kα(θ)k2 ≷ λ.
To find initial upper and lower bounds on θ, we note that(1/σmax (H) + θ)−1 kvk2 ≤ kα(θ)k2 ≤ (1/σmin (H) + θ)−1 kvk2
where σmax (H) denotes the maximum singular value of H and σmin (H) the minimum. To guaranteekα(θmax )k2 ≤ λ, we thus set θmax = kvk2 /λ − 1/σmax (H). Similarly, for θmin we see that so long asθ ≥ kvk2 /λ − 1/σmin (H) we have kα(θ)k2 ≥ λ. The fact that ∂ kxk2 = {z : kzk2 ≤ 1} when x = 0
implies that the solution for the original problem (21) is x = 0 if and only if kuk2 ≤ λ. We provide
pseudocode for solving (21) in Algorithm 4.
2139

D UCHI , H AZAN AND S INGER

5.4 ℓ∞ Regularization
We again let X = Rd but now choose ϕ(x) = λ kxk∞ . This type of update, similarly to ℓ2 , zeroes
groups of variables, which is handy in finding structurally sparse solutions for multitask or multiclass problems. Solving the ℓ∞ regularized problem amounts to1
min hu, xi + hx, Hxi + λ kxk∞ .
x
2

(23)

The dual of this problem is a modified ℓ1 -projection problem. As in the case of ℓ2 regularization,
we introduce an equality constrained variable z = x with associated Lagrange multipliers α ∈ Rd toobtain1
L (x, z, α) = hu, xi + hx, Hxi + λ kzk∞ + hα, x − zi .
2
Performing identical manipulations to the ℓ2 case, we take derivatives and get that x = −H −1 (u + α)and, similarly, unless kαk1 ≤ λ, infz L (x, z, α) = −∞. Thus the dual problem for (23) is1
max − (u + α)H −1 (u + α) s.t. kαk1 ≤ λ .
α2When H is diagonal we can find the optimal α∗ using the generalized ℓ1 -projection in Algorithm 3,then reconstruct the optimal x via x = −H −1 (u + α∗ ).
5.5 Mixed-norm Regularization
Finally, we combine the above results to show how to solve problems with matrix-valued inputs
X ∈ Rd×k , where X = [x1 · · · xd ]⊤ . We consider mixed-norm regularization, which is very usefulfor encouraging sparsity across several tasks (Obozinski et al., 2007). Now ϕ is an ℓ1 /ℓ p norm, thatis, ϕ(X) = λ ∑di=1 kxi k p . By imposing an ℓ1 -norm over p-norms of the rows of X, entire rows arenulled at once.
When p ∈ {2, ∞} and the proximal H in (18) is diagonal, the previous algorithms can be readily
used to solve the mixed norm problems. We simply maintain diagonal matrix information for eachof the rows x̄i of X separately, then solve one of the previous updates for each row independently.
We use this form of regularization in our experiments with multiclass prediction problems in thenext section.

6. Experiments
We performed experiments with several real world data sets with different characteristics: the ImageNet image database (Deng et al., 2009), the Reuters RCV1 text classification data set (Lewis
et al., 2004), the MNIST multiclass digit recognition problem, and the census income data set fromthe UCI repository (Asuncion and Newman, 2007). For uniformity across experiments, we focus on
the completely online (fully stochastic) optimization setting, in which at each iteration the learningalgorithm receives a single example. We measure performance using two metrics: the online loss
or error and the test set performance of the predictor the learning algorithm outputs at the end of asingle pass through the training data. We also give some results that show how imposing sparsity
constraints (in the form of ℓ1 and mixed-norm regularization) affects the learning algorithm’s performance. One benefit of the A DAG RAD framework is its ability to straightforwardly generalize to2140

A DAPTIVE S UBGRADIENT M ETHODS

ECAT
CCAT
GCAT
MCAT

RDA
.051 (.099).064 (.123).046 (.092).037 (.074)

FB.058 (.194).111 (.226).056 (.183).056 (.146)

A DAG RAD-RDA.044 (.086).053 (.105).040 (.080).035 (.063)

A DAG RAD-FB.044 (.238).053 (.276).040 (.225).034 (.176)

PA
.059
.107
.066
.053

AROW
.049
.061
.044
.039

Table 1: Test set error rates and proportion non-zero (in parenthesis) on Reuters RCV1.
domain constraints X 6= Rd and arbitrary regularization functions ϕ, in contrast to previous adaptive
online algorithms.
We experiment with RDA (Xiao, 2010), F OBOS(Duchi and Singer, 2009), adaptive RDA, adaptive F OBOS, the Passive-Aggressive (PA) algorithm (Crammer et al., 2006), and AROW (Crammer
et al., 2009). To remind the reader, PA is an online learning procedure with the updatext+1 = argmin [1 − yt hzt , xi]+ +x

λkx − xt k22 ,2

where λ is a regularization parameter. PA’s update is similar to the update employed by AROW(see (9)), but the latter maintains second order information on x. By using a representer theorem
it is also possible to derive efficient updates for PA and AROW when the loss is the logistic loss,
log(1 + exp(−yt hzt , xt i)). We thus we compare the above six algorithms using both hinge andlogistic loss.
6.1 Text Classification
The Reuters RCV1 data set consists of a collection of approximately 800,000 text articles, eachof which is assigned multiple labels. There are 4 high-level categories, Economics, Commerce,
Medical, and Government (ECAT, CCAT, MCAT, GCAT), and multiple more specific categories.
We focus on training binary classifiers for each of the four major categories. The input features
we use are 0/1 bigram features, which, post word stemming, give data of approximately 2 million
dimensions. The feature vectors are very sparse, however, and most examples have fewer than 5000
non-zero features.
We compare the twelve different algorithms mentioned in the prequel as well as variants of
F OBOS and RDA with ℓ1 -regularization. We summarize the results of the ℓ1 -regularized runs aswell as AROW and PA in Table 1. The results for both hinge and logistic losses are qualitatively
and quantitatively very similar, so we report results only for training with the hinge loss in Table 1.
Each row in the table represents the average of four different experiments in which we hold out 25%of the data for a test set and perform an online pass on the remaining 75% of the data. For RDA
and F OBOS, we cross-validate the stepsize parameter η by simply running multiple passes and thenchoosing the output of the learner that had the fewest mistakes during training. For PA and AROW
we choose λ using the same approach. We use the same regularization multiplier on the ℓ1 term for
RDA and F OBOS, selected so that RDA achieved approximately 10% non-zero predictors.
It is evident from the results presented in Table 1 that the adaptive algorithms (AROW and A DA G RAD) are far superior to non-adaptive algorithms in terms of error rate on test data. The A DA G RAD algorithms naturally incorporate sparsity as well since they are run with ℓ1 -regularization,
though RDA has significantly higher sparsity levels (PA and AROW do not have any sparsity). Furthermore, although omitted from the table to avoid clutter, in every test with the RCV1 corpus, the2141

D UCHI , H AZAN AND S INGERAlg.
A DAG RAD RDA
AROW
PA
RDA

Avg. Prec.
0.60220.58130.55810.5042

P@1
0.85020.85970.84550.7496

P@30.83070.83690.81840.7185P@50.81300.81650.79570.6950P@10
0.78110.78160.75760.6545

Prop. nonzero
0.72671.00001.00000.8996

Table 2: Test set precision for ImageNet
adaptive algorithms outperformed the non-adaptive algorithms. Moreover, both A DAG RAD-RDA
and A DAG RAD-Fobos outperform AROW on all the classification tasks. Unregularized RDA andF OBOS attained similar results as did the ℓ1 -regularized variants (of course without sparsity), but
we omit the results to avoid clutter and because they do not give much more understanding.
6.2 Image Ranking
ImageNet (Deng et al., 2009) consists of images organized according to the nouns in the WordNet
hierarchy, where each noun is associated on average with more than 500 images collected fromthe web. We selected 15,000 important nouns from the hierarchy and conducted a large scale image ranking task for each noun. This approach is identical to the task tackled by Grangier andBengio (2008) using the Passive-Aggressive algorithm. To solve this problem, we train 15,000
ranking machines using Grangier and Bengio’s visterms features, which represent patches in an image with 79-dimensional sparse vectors. There are approximately 120 patches per image, resulting
in a 10,000-dimensional feature space.
Based on the results in the previous section, we focus on four algorithms for solving this task:
AROW, A DAG RAD with RDA updates and ℓ1 -regularization, vanilla RDA with ℓ1 , and PassiveAggressive. We use the ranking hinge loss, which is [1 − hx, z1 − z2 i]+ when z1 is ranked above
z2 . We train a ranker xc for each of the image classes individually, cross-validating the choice ofinitial stepsize for each algorithm on a small held-out set. To train an individual ranker for classc, at each step of the algorithm we randomly sample a positive image z1 for the category c andan image z2 from the training set (which with high probability is a negative example for class c)
and perform an update on the example z1 − z2 . We let each algorithm take 100,000 such steps foreach image category, we train four sets of rankers with each algorithm, and the training set includes
approximately 2 million images.
For evaluation, we use a distinct test set of approximately 1 million images. To evaluate a set of
rankers, we iterate through all 15,000 classes in the data set. For each class we take all the positive
image examples in the test set and sample 10 times as many negative image examples. Following
Grangier and Bengio, we then rank the set of positive and negative images and compute precisionat-k for k = {1, . . . , 10} and the average precision for each category. The precision-at-k is definedas the proportion of examples ranked in the top k for a category c that actually belong to c, andthe average precision is the average of the precisions at each position in which a relevant pictureappears. Letting Pos(c) denote the positive examples for category c and p(i) denote the position of
the ith returned picture in list of images sorted by inner product with xc , the average precision is|Pos(c)|1i∑ p(i) .
|Pos(c)| i=12142

A DAPTIVE S UBGRADIENT M ETHODS

10000
9000
8000

Mistakes

7000

PA
Ada RDARDA
Ada RDA L1/L2
RDA L1/L2

6000
5000
4000
3000
2000
1000
0
0

1

2

3
4
Examples seen

5

6
4

x 10

Figure 5: Learning curves on MNIST
We compute the mean of each measurement across all classes, performing this twelve times foreach of the sets of rankers trained. Table 2 summarizes our results. We do not report variance as thevariance was on the order of 10−5 for each algorithm. One apparent characteristic to note from thetable is that A DAG RAD RDA achieves higher levels of sparsity than the other algorithms—using
only 73% of the input features it achieves very high performance. Moreover, it outperforms all thealgorithms in average precision. AROW has better results than the other algorithms in terms ofprecision-at-k for k ≤ 10, though A DAG RAD’s performance catches up to and eventually surpasses
AROW’s as k grows.
6.3 Multiclass Optical Character Recognition
In the well-known MNIST multiclass classification data set, we are given 28 × 28 pixel images ai ,and the learner’s task is to classify each image as a digit in {0, . . . , 9}. Linear classifiers do notwork well on a simple pixel-based representation. Thus we learn classifiers built on top of a kernel
machine with Gaussian kernels, as do Duchi and Singer (2009), which gives a different (and nonsparse) structure to the feature space in contrast to our previous experiments.
In particular,
for the21ith example and jth feature, the feature value is zi j = K(ai , a j ) , exp − 2σ2 ai − a j 2 . We use a
support set of approximately 3000 images to compute the kernels and trained multiclass predictors,
which consist of one vector xc ∈ R3000 for each class c, giving a 30,000 dimensional problem.
There is no known multiclass AROW algorithm. We therefore compare adaptive RDA with andwithout mixed-norm ℓ1 /ℓ2 and ℓ1 /ℓ∞ regularization (see Section 5.5), RDA, and multiclass Passive
Aggressive to one another using the multiclass hinge loss (Crammer et al., 2006). For each algorithm
we used the first 5000 of 60,000 training examples to choose the stepsize η (for RDA) and λ (forPA).
In Figure 5, we plot the learning curves (cumulative mistakes made) of multiclass PA, RDA,
RDA with ℓ1 /ℓ2 regularization, adaptive RDA, and adaptive RDA with ℓ1 /ℓ2 regularization (ℓ1 /ℓ∞
2143

D UCHI , H AZAN AND S INGER

PA
Ada-RDARDA
Ada-RDA λ = 5 · 10−4
RDA λ = 5 · 10−4
Ada-RDA λ = 10−3
RDA λ = 10−3

Test error rate
0.062
0.066
0.108
0.100
0.138
0.137
0.192

Prop. nonzero
1.000
1.000
1.000
0.569
0.878
0.144
0.532

Table 3: Test set error rates and sparsity proportions on MNIST. The scalar λ is the multiplier on
the ℓ1 /ℓ2 regularization term.
is similar). From the curves, we see that Adaptive RDA seems to have similar performance to PA,
and the adaptive versions of RDA are vastly superior to their non-adaptive counterparts. Table 3
further supports this, where we see that the adaptive RDA algorithms outperform their non-adaptive
counterparts both in terms of sparsity (the proportion of non-zero rows) and test set error rates.
6.4 Income Prediction
The KDD census income data set from the UCI repository (Asuncion and Newman, 2007) contains
census data extracted from 1994 and 1995 population surveys conducted by the U.S. Census Bureau.
The data consists of 40 demographic and employment related variables which are used to predictwhether a respondent has income above or below $50,000. We quantize each feature into bins (5per feature for continuous features) and take products of features to give a 4001 dimensional feature
space with 0/1 features. The data is divided into a training set of 199,523 instances and test set of
99,762 test instances.
As in the prequel, we compare AROW, PA, RDA, and adaptive RDA with and without ℓ1 regularization on this data set. We use the first 10,000 examples of the training set to select thestep size parameters λ for AROW and PA and η for RDA. We perform ten experiments on random
shuffles of the training data. Each experiment consists of a training pass through some proportionof the data (.05, .1, .25, .5, or the entire training set) and computing the test set error rate of thelearned predictor. Table 4 and Figure 6 summarize the results of these experiments. The varianceof the test error rates is on the order of 10−6 so we do not report it. As earlier, the table and figure
make it clear that the adaptive methods (AROW and A DAG RAD-RDA) give better performance
than non-adaptive methods. Further, as detailed in the table, the A DAG RAD methods can give
extremely sparse predictors that still give excellent test set performance. This is consistent withthe experiments we have seen to this point, where A DAG RAD gives sparse but highly accurate
predictors.
6.5 Experiments with Sparsity-Accuracy Tradeoffs
In our final set of experiments, we investigate the tradeoff between the level of sparsity and theclassification accuracy for the A DAG RAD-RDA algorithms. Using the same experimental setup
as for the initial text classification experiments described in Section 6.1, we record the averagetest-set performance of A DAG RAD-RDA versus the proportion of features that are non-zero in thepredictor A DAG RAD outputs after a single pass through the training data. To achieve this, we run
2144

A DAPTIVE S UBGRADIENT M ETHODS

0.056

AROW
PA
RDA
Ada RDA

Test error rate

0.054

0.052

0.05

0.048

0.046

0.044
0

0.2

0.4
0.6
Proportion train

0.8

1

Figure 6: Test set error rates as function of proportion of training data seen on Census Income data
set.
Prop. Train
AROW
PA
RDA
Ada-RDAℓ1 RDAℓ1 Ada-RDA

0.05
0.049
0.055
0.055
0.053
0.056 (0.075)
0.052 (0.062)

0.10
0.048
0.052
0.054
0.051
0.054 (0.066)
0.051 (0.053)

0.25
0.046
0.050
0.052
0.049
0.053 (0.058)
0.050 (0.044)

0.50
0.045
0.049
0.051
0.048
0.052 (0.053)
0.050 (0.040)

1.00
0.044
0.048
0.050
0.047
0.051 (0.050)
0.049 (0.037)

Table 4: Test set error rates as function of proportion of training data seen (proportion of non-zeros
in parenthesis where appropriate) on Census Income data set.
A DAG RAD with ℓ1 -regularization, and we sweep the regularization multiplier λ from 10−8 to 10−1 .
These values result in predictors ranging from a completely dense predictor to an all-zeros predictor,
respectively.
We summarize our results in Figure 7, which shows the test set performance of A DAG RAD
for each of the four categories ECAT, CCAT, GCAT, and MCAT. Within each plot, the horizontal
black line labeled AROW designates the baseline performance of AROW on the text classificationtask, though we would like to note that AROW generates fully dense predictors. The plots all
portray a similar story. With high regularization values, A DAG RAD exhibits, as expected, poor
performance as it retains no predictive information from the learning task. Put another way, when
the regularization value is high A DAG RAD is confined to an overly sparse predictor which exhibits
poor generalization. However, as the regularization multiplier λ decreases, the learned predictor
becomes less sparse and eventually the accuracy of A DAG RAD exceeds AROW’s accuracy. It isinteresting to note that for these experiments, as soon as the predictor resulting from a single pass2145

D UCHI , H AZAN AND S INGER

ECAT

CCAT

0.16

0.16

AdaGrad
AROW

0.12

0.1

0.08

0.06

0.04

0.02 −510

AdaGrad
AROW

0.14

Test−set error rateTest−set error rate

0.14

0.12

0.1

0.08

0.06

0.04

−4

10

−3

10

−2

10

−1

10

0.02 −5
10

0

10

−4

10

Proportion non−zero

−3

10

GCAT

0

10

0.16

AdaGrad
AROW

0.14

AdaGrad
AROW

0.14

0.12

Test−set error rateTest−set error rate

−1

10

MCAT

0.16

0.1

0.08

0.06

0.04

0.02 −510−2

10Proportion non−zero

0.12

0.1

0.08

0.06

0.04

−4

10

−3

10

−2

10

−1

10

0

10

0.02 −510

Proportion non−zero

−4

10

−3

10

−2

10

−1

10

0

10

Proportion non−zeroFigure 7: Test set error rates as a function of proportion of non-zeros in predictor x output by A DA G RAD (AROW plotted for reference).
through the data has more than 1% non-zero coefficients, A DAG RAD’s performance matches that of
AROW. We also would like to note that the variance in the test-set error rates for these experiments
is on the order of 10−6 , and we thus do not draw error bars in the graphs. The performance of
A DAG RAD as a function of regularization for other sparse data sets, especially in relation to that of
AROW, was qualitatively similar to this experiment.

7. Conclusions
We presented a paradigm that adapts subgradient methods to the geometry of the problem at hand.
The adaptation allows us to derive strong regret guarantees, which for some natural data distributions
achieve better performance guarantees than previous algorithms. Our online regret bounds can benaturally converted into rate of convergence and generalization bounds (Cesa-Bianchi et al., 2004).
Our experiments show that adaptive methods, specifically A DAG RAD-F OBOS, A DAG RAD-RDA,
and AROW clearly outperform their non-adaptive counterparts. Furthermore, the A DAG RAD fam2146

A DAPTIVE S UBGRADIENT M ETHODS

ily of algorithms naturally incorporates regularization and gives very sparse solutions with similar
performance to dense solutions. Our experiments with adaptive methods use a diagonal approximation to the matrix obtained by taking outer products of subgradients computed along the run of thealgorithm. It remains to be tested whether using the full outer product matrix can further improve
performance.
To conclude we would like to underscore a possible elegant generalization that interpolates
between full-matrix proximal functions and diagonal approximations using block diagonal matrices.
k⊤ ⊤diSpecifically, for v ∈ Rd let v = [v⊤[1] · · · v[k] ] where v[i] ∈ R are subvectors of v with ∑i=1 di = d.
We can define the associated block-diagonal approximation to the outer product matrix ∑tτ=1 gτ g⊤τbygτ,[1] g⊤
0
···
0
τ,[1]


..

t 
.
0
gτ,[2] g⊤
0


τ,[2]
Gt = ∑ 
.
..
.
.


.
.
τ=1 
.
.
0
.

0
···
0 gτ,[k] g⊤τ,[k]In this case, a combination of Theorems 5 and 7 gives the next corollary.

Corollary 12 Let Gt be the block-diagonal outer productD matrix
E defined above and the sequence
1/2{xt } be defined by the RDA update of (3) with ψt (x) = x, Gt x . Then, for any x∗ ∈ X ,Rφ (T ) ≤

2
1
1/2
1/2
∗
max x[i]tr(GT ) + η tr(GT ).
η i
2

A similar bound holds for composite mirror-descent updates, and it is straightforward to get infimalequalities similar to those in Corollary 11 with the infimum taken over block-diagonal matrices.
Such an algorithm can interpolate between the computational simplicity of the diagonal proximal
functions and the ability of full matrices to capture correlation in the gradient vectors.
A few open questions stem from this line of research. The first is whether we can efficiently
use full matrices in the proximal functions, as in Section 4. A second open issue is whether nonEuclidean proximal functions, such as the relative entropy, can be used. We also think that thestrongly convex case—when ft or ϕ is strongly convex—presents interesting challenges that wehave not completely resolved. We hope to investigate both empirical and formal extensions of thiswork in the near future.

Acknowledgments
There are many people to whom we owe our sincere thanks for this research. Fernando Pereira
helped push us in the direction of working on adaptive online methods and has been a constantsource of discussion and helpful feedback. Samy Bengio provided us with a processed version of
the ImageNet data set and was instrumental in helping to get our experiments running, and Adam
Sadovsky gave many indispensable coding suggestions. The anonymous reviewers also gave several
suggestions that improved the quality of the paper. Lastly, Sam Roweis was a sounding board forsome of our earlier ideas on the subject, and we will miss him dearly.
2147

D UCHI , H AZAN AND S INGER

Appendix A. Full Matrix Motivating Example
As in the diagonal case, as the adversary we choose ε > 0 and on rounds t = 1, . . . , η2 /ε2 play thevector ±v1 . After the first η2 /ε2 rounds, the adversary simply cycles through the vectors v2 , . . . , vd .
Thus, for Zinkevich’s projected gradient, we have xt = αt,1 v1 for some multiplier αt,1 > 0 whent ≤ η2 /ε2 . After the first η2 /ε2 rounds, we perform the updates


η
√
√
vi
xt+1 = Πkxk ≤ d xt +2t√for some index i, but as in the diagonal case, η/ t ≤ ε, and by orthogonality of vi , v j , we have
xt = V αt for some αt  0, and the projection step can only shrink the multiplier αt,i for index
i. Thus, each coordinate incurs loss at least 1/(2ε), and projected gradient descent suffers lossesΩ(d/ε).
On the other hand, A DAG RAD suffers loss at most d. Indeed, since g1 = v1 and kv1 k2 = 1, we
1

†
⊤
⊤
2
have G21 = v1 v⊤
1 v1 v1 = v1 v1 = G1 , so G1 = G1 = G1 , andx2 = x1 + G†1 = x1 + v1 v⊤
1 v1 = x1 + v1 .
Since hx2 , v1 i = 1, we see that A DAG RAD suffers no loss (and Gt = G1 ) until a vector zt = ±vi for
i 6= 1 is played by the adversary. However, an identical argument shows that Gt is simply updated⊤
to v1 v⊤
1 + vi . Indeed, an inductive argument shows that until all the
1 + vi vi , in which case xt = v√
vectors vi are seen, we have kxt k2 < d by orthogonality, and eventually we havesdd√xt = ∑ vi and kxt k2 = ∑ kvi k22 = di=1so that xt ∈ X = {x : kxk2 ≤1 and suffer no loss.

√

i=1d} for A DAG RAD for all t. All future predictions thus achieve margin

Appendix B. Technical Lemmas
Lemma 13 Let A  B  0 be symmetric d × d PSD matrices. Then A1/2  B1/2 .
Proof This is Example 3 of Davis (1963). We include a proof for convenience of the reader.
Let λ be any eigenvalue (with corresponding eigenvector x) of A1/2 − B1/2 ; we show that λ ≥ 0.
2
Clearly A1/2 x − λx = B1/2 x. Taking the inner product of both sides with A1/2 x, we have A1/2 x 2 −λ A1/2 x, x = A1/2 x, B1/2 x . We use the Cauchy-Schwarz inequality:
A1/2 x2
2DE− λ A1/2 x, x ≤ A1/2 x

2

B1/2 x2=p2hAx, xi hBx, xi ≤ hAx, xi = A1/2 x
2

where the last inequality follows from the assumption that A  B. Thus we must have λ A1/2 x, x ≥0, which implies λ ≥ 0.

The gradient of the function tr(X p ) is easy to compute for integer values of p. However, when p isreal we need the following lemma. The lemma tacitly uses the fact that there is a unique positive
semidefinite X p when X  0 (Horn and Johnson, 1985, Theorem 7.2.6).
2148

A DAPTIVE S UBGRADIENT M ETHODS

Lemma 14 Let p ∈ R and X ≻ 0. Then ∇X tr(X p ) = pX p−1 .
Proof We do a first order expansion of (X + A) p when X ≻ 0 and A is symmetric. Let X = UΛU ⊤be the symmetric eigen-decomposition of X and V DV ⊤ be the decomposition of Λ−1/2U ⊤ AUΛ−1/2 .
Then(X + A) p = (UΛU ⊤ + A) p = U(Λ +U ⊤ AU) pU ⊤ = UΛ p/2 (I + Λ−1/2U ⊤ AUΛ−1/2 ) p Λ p/2U ⊤= UΛ p/2V ⊤ (I + D) pV Λ p/2U ⊤ = UΛ p/2V ⊤ (I + pD + o(D))V Λ p/2U ⊤= UΛ pU ⊤ + pUΛ p/2V ⊤ DV Λ p/2U ⊤ + o(UΛ−/2V ⊤ DV Λ p/2U ⊤ )= X p +UΛ(p−1)/2U ⊤ AUΛ(p−1)/2U ⊤ + o(A) = X p + pX (p−1)/2 AX (p−1)/2 + o(A).
In the above, o(A) is a matrix that goes to zero faster than A → 0, and the second line follows via afirst-order Taylor expansion of (1 + di ) p . From the above, we immediately havetr((X + A) p ) = tr X p + p tr(X p−1 A) + o(tr A),
which completes the proof.

Appendix C. Proof of Lemma 4
We prove the lemma by considering an arbitrary real-valued sequence {ai } and its vector representation a1:i = [a1 · · · ai ]. We are next going to show thatT

at2
≤ 2 ka1:T k2 ,∑t=1 ka1:t k2(24)

where we define 00 = 0. We use induction on T to prove inequality (24). For T = 1, the inequality
trivially holds. Assume the bound (24) holds true for T − 1, in which case
T

T −1

a2

a2

a2

a2

∑ ka1:tt k = ∑ ka1:tt k + ka1:TT k ≤ 2 ka1:T −1 k2 + ka1:TT k ,
2

t=1

2

t=1

2

2

T
where the inequalityq
follows from the inductive hypothesis. We define bT = ∑t=1
at2 and use con√cavity to obtain that bT − a2T ≤ bT − a2T 2√1b so long as bT − a2T ≥ 0.2 Thus,T

2 ka1:T −1 k2 +

q
p
a2T
a2
= 2 bT − a2T + √ T ≤ 2 bT = 2 ka1:T k2 .
ka1:T k2bT

Having proved the bound (24), we note that by construction that st,i = kg1:t,i k2 , so
2dgt,i
∑ gt , diag(st ) gt = ∑ ∑ kg1:t,i k ≤ 2 ∑ kg1:T,i k2 .
2
t=1 i=1
t=1
i=1
T

−1

T

d

2. We note that we use an identical technique in the full-matrix case. See Lemma 8.

2149

D UCHI , H AZAN AND S INGER

Appendix D. Proof of Lemmas 8 and 9
We begin with the more difficult proof of Lemma 8.
Proof of Lemma 8 The core of the proof is based on the concavity of the function tr(A1/2 ). However, careful analysis is required as A might not be strictly positive definite. We also use the previous
lemma which implies that the gradient of tr(A1/2 ) is 12 A−1/2 when A ≻ 0.
First, A p is matrix-concave for A ≻ 0 and 0 ≤ p ≤ 1 (see, for example, Corollary 4.1 in Ando,
1979 or Theorem 16.1 in Bondar, 1994). That is, for A, B ≻ 0 and α ∈ [0, 1] we have
(αA + (1 − α)B) p  αA p + (1 − α)B p .

(25)

Now suppose simply A, B  0 (but neither is necessarily strict). Then for any δ > 0, we haveA + δI ≻ 0 and B + δI ≻ 0 and therefore(α(A + δI) + (1 − α)(B + δI)) p  α(A + δI) p + (1 − α)(B + δI) p  αA p + (1 − α)B p ,where we used Lemma 13 for the second matrix inequality. Moreover, αA + (1 − α)B + δI →αA + (1 − α)B as δ → 0. Since A p is continuous (when we use the unique PSD root), this line ofreasoning proves that (25) holds for A, B  0. Thus, we proved thattr((αA + (1 − α)B) p ) ≥ α tr(A p ) + (1 − α) tr(B p ) for 0 ≤ p ≤ 1 .
Recall now that Lemma 14 implies that the gradient of tr(A1/2 ) is 12 A−1/2 when A ≻ 0. Therefore, from the concavity of A1/2 and the form of its gradient, we can use the standard first-orderinequality for concave functions so that for any A, B ≻ 0,
1
tr(A1/2 ) ≤ tr(B1/2 ) + tr(B−1/2 (A − B)) .
2(26)

Let A = B − νgg⊤  0 and suppose only that B  0. We must take some care since B−1/2 may
not necessarily exist, and the above inequality does not hold true in the pseudo-inverse sense when
B 6≻ 0. However, for any δ > 0 we know that 2∇B tr((B+δI)1/2 ) = (B+δI)−1/2 , and A−B = −νgg⊤ .
From (26) and Lemma 13, we have
2 tr(B − tgg⊤ )1/2 = 2 tr(A1/2 ) ≤ 2 tr((A + δI)1/2 )≤ 2 tr(B + δI)1/2 − ν tr((B + δI)−1/2 gg⊤ ) .

(27)

Note that g ∈ Range(B), because if it were not, we could choose some u with Bu = 0 and hg, ui 6= 0,
which would give u, (B − cgg⊤ )u = −c hg, ui2 < 0, a contradiction. Now let B = V diag(λ)V ⊤ bethe eigen-decomposition of B. Since g ∈ Range(B), pg⊤ (B + δI)−1/2 g = g⊤V diag 1/ λi + δ V ⊤ g
=

1

∑ λi
∑ √λi + δ (g⊤ vi )2 −→
δ↓0

−1/2

(g⊤ vi )2 = g⊤ (B† )1/2 g .

i:λi >0i:λi >0Thus, by taking δ ↓ 0 in (27), and since both tr(B + δI)1/2 and tr((B + δI)−1/2 gg⊤ ) are evidently
continuous in δ, we complete the proof.

2150

A DAPTIVE S UBGRADIENT M ETHODS

Proof of Lemma 9 We begin by noting that δ2 I  gg⊤ , so from Lemma 13 we get (A + gg⊤ )1/2 
2 1/2(A√ A and I are simultaneously diagonalizable, we can generalize the inequality√ + δ I) √. Sincea + b ≤ a + b, which holds for a, b ≥ 0, to positive semi-definite matrices, thus,(A + δ2 I)1/2  A1/2 + δI .
Therefore, if A + gg⊤ is of full rank, we have (A + gg⊤ )−1/2  (A1/2 + δI)−1 (Horn and Johnson,
1985, Corollary 7.7.4(a)). Since g ∈ Range((A + gg⊤ )1/2 ), we can apply an analogous limiting argument to the one used in the proof of Lemma 8 and discard all zero eigenvalues of A + gg⊤ , which
completes the lemma.

Appendix E. Solution to Problem (15)We prove here a technical lemma that is useful in characterizing the solution of the optimization
problem below. Note that the second part of the lemma implies that we can treat the inverse of thesolution matrix S−1 as S† . We consider solving
min tr(S−1 A) subject to S  0, tr(S) ≤ c where A  0 .
S

1

(28)

1

Lemma 15 If A is of full rank, then the minimizer of (28) is S = cA 2 / tr(A 2 ). If A is not of full rank,
11then setting S = cA 2 / tr(A 2 ) givestr(S† A) = inf tr(S−1 A) : S  0, tr(S) ≤ c .
S

1

In either case, tr(S† A) = tr(A 2 )2 /c.

Proof Both proofs rely on constructing the Lagrangian for (28). We introduce θ ∈ R+ for the trace
constraint and Z  0 for the positive semidefinite constraint on S. In this case, the Lagrangian is

L (S, θ, Z) = tr(S−1 A) + θ(tr(S) − c) − tr(SZ).
The derivative of L with respect to S is−S−1 AS−1 + θI − Z.

(29)

If S is full rank, then to satisfy the generalized complementarity conditions for the problem (Boydand Vandenberghe, 2004), we must have Z = 0. Therefore, we get S−1 AS−1 = θI. We now can
1multiply by S on the right and the left to get that A = θS2 , which implies that S ∝ A 2 . If A is of fullrank, the optimal solution for S ≻ 0 forces θ to be positive so that tr(S) = c. This yields the solution
1
1
S = cA 2 / tr(A 2 ). In order to verify optimality of this solution, we set Z = 0 and θ = c−2 tr(A1/2 )2which gives ∇S L (S, θ, Z) = 0, as is indeed required.
Suppose now that A is not full rank and that


Λ 0 ⊤
A=Q
Q
0 0
2151

D UCHI , H AZAN AND S INGER

is the eigen-decomposition of A. Let n be the dimension of the null-space of A (so the rank of A isd − n). Define the variables


 1


 1
c
1
0 0
Λ2 0
Λ2 0
⊤
Q , S(δ) =Q⊤ .
Z(θ) =, S(θ, δ) = √ Q
Q
1
0 θI
0 δI
0 δI
θ
tr(A 2 ) + δn
It is easy to see that tr S(δ) = c, and1

11

lim tr(S(δ)−1 A) = tr(S(0)† A) = tr(A 2 ) tr(Λ 2 )/c = tr(A 2 )2 /c.

δ→0

Further, let g(θ) = infS L (S, θ, Z(θ)) be the dual of (28). From the above analysis and (29), it isevident that



 −1
1
0 0
0
Λ 2 ΛΛ− 2⊤−1−1
Q + θI −= 0.
−S(θ, δ) AS(θ, δ) + θI − Z(θ) = −θQ0 θI0δ−2 I · 0
So S(θ, δ) achieves the infimum in the dual for any δ > 0, tr(S(0)Z(θ)) = 0, and
√
√
√
1
1
g(θ) = θ tr(Λ 2 ) + θ tr(Λ 2 ) + θδn − θc.
1

1

1

1

Setting θ = tr(Λ 2 )2 /c2 gives g(θ) = tr(Λ 2 )2 /c − δn tr(Λ 2 )/c. Taking δ → 0 gives g(θ) = tr(A 2 )2 /c,
1which means that limδ→0 tr(S(δ)−1 A) = tr(A 2 )2 /c = g(θ). Thus the duality gap for the original
problem is 0 so S(0) is the limiting solution.
1
1
The last statement of the lemma is simply plugging S† = (A† ) 2 tr(A 2 )/c in to the objective being
minimized.

Appendix F. Proofs of Propositions 2 and 3
We begin with the proof of Proposition 2. The proof essentially builds upon Xiao (2010) and
Nesterov (2009), with some modification to deal with the indexing of ψt . We include the proof forcompleteness.
Proof of Proposition 2 Define ψt∗ to be the conjugate dual of tϕ(x) + ψt (x)/η:no1ψt∗ (g) = sup hg, xi − tϕ(x) − ψt (x) .
ηx∈XSince ψt /η is 1/η-strongly convex with respect to the norm k·kψt , the function ψt∗ has η-Lipschitz
continuous gradients with respect to k·kψt∗ :k∇ψt∗ (g1 ) − ∇ψt∗ (g2 )kψt ≤ η kg1 − g2 kψt∗

(30)

for any g1 , g2 (see, e.g., Nesterov, 2005, Theorem 1 or Hiriart-Urruty and Lemaréchal, 1996, Chapter X). Further, a simple argument with the fundamental theorem of calculus gives that if f has
L-Lipschitz gradients, f (y) ≤ f (x) + h∇ f (x), y − xi + (L/2) ky − xk2 , and


1
∗
(31)
∇ψt (g) = argmin − hg, xi + tϕ(x) + ψt (x) .
η
x∈X
2152

A DAPTIVE S UBGRADIENT M ETHODS

Using the bound (30) and identity (31), we can give the proof of the corollary. Indeed, letting
gt ∈ ∂ ft (xt ) and defining zt = ∑tτ=1 gτ , we haveT∑ ft (xt ) + ϕ(xt ) − ft (x∗ ) − ϕ(x∗ )t=1

T≤ ∑ hgt , xt − x∗ i − ϕ(x∗ ) + ϕ(xt )
t=1

(

)
1
≤ ∑ hgt , xt i + ϕ(xt ) + sup − ∑ hgt , xi − T ϕ(x) − ψT (x) + ψT (x∗ )
η
x∈X
t=1
t=1
T

=

T

T
1
ψT (x∗ ) + ∑ hgt , xt i + ϕ(xt ) + ψ∗T (−zT ) .
η
t=1

Since ψt+1 ≥ ψt , it is clear thatT1ψ∗T (−zT ) = − ∑ hgt , xT +1 i − T ϕ(xT +1 ) − ψT (xT +1 )
η
t=1
T
1
≤ − ∑ hgt , xT +1 i − (T − 1)ϕ(xT +1 ) − ϕ(xT +1 ) − ψT −1 (xT +1 )
η
t=1


1
≤ sup − hzT , xi − (T − 1)ϕ(x) − ψT −1 (x) − ϕ(xT +1 ) = ψ∗T −1 (−zT ) − ϕ(xT +1 ).
ηx∈X

The Lipschitz continuity of ∇ψt∗ , the identity (31), and the fact that zT − zT −1 = −gT giveT∑ ft (xt ) + ϕ(xt+1) − ft (x∗ ) − ϕ(x∗ )

t=1

≤

T
1
ψT (x∗ ) + ∑ hgt , xt i + ϕ(xt+1 ) + ψ∗T −1 (−zT ) − ϕ(xT +1 )
η
t=1

≤

T
1
ψT (x∗ ) + ∑ hgt , xt i + ϕ(xt+1 ) − ϕ(xT +1 )η
t=1

η
kgT k2ψ∗
T −1
2
T −1
η
1
= ψT (x∗ ) + ∑ hgt , xt i + ϕ(xt+1 ) + ψ∗T −1 (−zT −1 ) + kgT k2ψ∗ .
T −1
η
2
t=1
+ ψ∗T −1 (−zT −1 ) − ∇ψ∗T −1 (zT −1 ), gT +We can repeat the same sequence of steps that gave the last equality to see thatT∑ ft (xt ) + ϕ(xt+1 ) − ft (x∗ ) − ϕ(x∗ ) ≤

t=1

1
η T
ψT (x∗ ) + ∑ kgt k2ψ∗ + ψ∗0 (−z0 ).
t−1
η
2 t=1

Recalling that x1 = argminx∈X {ϕ(x)} and that ψ∗0 (0) = 0 completes the proof.
We now turn to the proof of Proposition 3. We begin by stating and fully proving an (essentially)
immediate corollary to Lemma 2.3 of Duchi et al. (2010).
2153

D UCHI , H AZAN AND S INGER

Lemma 16 Let {xt } be the sequence defined by the update (4) and assume that Bψt (·, ·) is strongly
convex with respect to a norm k·kψt . Let k·kψt∗ be the associated dual norm. Then for any x∗ ,
η ( ft (xt ) − ft (x∗ )) + η (ϕ(xt+1 ) − ϕ(x∗ )) ≤ Bψt (x∗ , xt ) − Bψt (x∗ , xt+1 ) +η2 ′2f (xt ) ψ∗t
2 t

Proof The optimality of xt+1 for (4) implies for all x ∈ X and ϕ′ (xt+1 ) ∈ ∂ϕ(xt+1 )x − xt+1 , η f ′ (xt ) + ∇ψt (xt+1 ) − ∇ψt (xt ) + ηϕ′ (xt+1 ) ≥ 0.

(32)

In particular, this obtains for x = x∗ . From the subgradient inequality for convex functions, we haveft (x∗ ) ≥ ft (xt ) + h ft′ (xt ), x∗ − xt i, or ft (xt ) − ft (x∗ ) ≤ h ft′ (xt ), xt − x∗ i, and likewise for ϕ(xt+1 ). We
thus haveη [ ft (xt ) + ϕ(xt+1 ) − ft (x∗ ) − ϕ(x∗ )]≤ η xt − x∗ , ft′ (xt ) + η xt+1 − x∗ , ϕ′ (xt+1 )= η xt+1 − x∗ , ft′ (xt ) + η xt+1 − x∗ , ϕ′ (xt+1 ) + η xt − xt+1 , ft′ (xt )=x∗ − xt+1 , ∇ψt (xt ) − ∇ψt (xt+1 ) − η ft′ (xt ) − ηϕ′ (xt+1 )+ hx∗ − xt+1 , ∇ψt (xt+1 ) − ∇ψt (xt )i + η xt − xt+1 , ft′ (xt ) .
Now, by (32), the first term in the last equation is non-positive. Thus we have that
η [ ft (xt ) + ϕ(xt+1 ) − ft (x∗ ) − ϕ(x∗ )]≤ hx∗ − xt+1 , ∇ψt (xt+1 ) − ∇ψt (xt )i + η xt − xt+1 , ft′ (xt )= Bψt (x∗ , xt ) − Bψt (xt+1 , xt ) − Bψt (x∗ , xt+1 ) + η xt − xt+1 , ft′ (xt )D 1
E√= Bψt (x∗ , xt ) − Bψt (xt+1 , xt ) − Bψt (x∗ , xt+1 ) + η η− 2 (xt − xt+1 ), η ft′ (xt )

1
η2 ′
2
≤ Bψt (x∗ , xt ) − Bψt (xt+1 , xt ) − Bψt (x∗ , xt+1 ) + kxt − xt+1 k2ψt +ft (xt ) ψ∗
t
2
2
2
η
2
≤ Bψt (x∗ , xt ) − Bψt (x∗ , xt+1 ) +ft′ (xt ) ψ∗ .
t
2

In the above, the first equality follows from simple algebra with Bregman divergences, the second
to last inequality follows from Fenchel’s inequality applied to the conjugate functions 21 k·k2ψt and
2
1
2 k·kψt∗ (Boyd and Vandenberghe, 2004, Example 3.27), and the last inequality follows from theassumed strong convexity of Bψt with respect to the norm k·kψt .
Proof of Proposition 3 Sum the equation in the conclusion of Lemma 16.

Appendix G. Derivations of Algorithms
In this appendix, we give the formal derivations of the solution to the A DAG RAD update for ℓ1 regularization and projection to an ℓ1 -ball, as described originally in Section 5.
2154

A DAPTIVE S UBGRADIENT M ETHODS

G.1 ℓ1 -regularization
We give the derivation for the primal-dual subgradient update, as composite mirror-descent is entirely similar. We need to solve update (3), which amounts tomin η hḡt , xi +x1
1
δ kxk22 + hx, diag(st )xi + ηλ kxk1 .
2t
2t

Let x̂ denote the optimal solution of the above optimization problem. Standard subgradient calculus
implies that when |ḡt,i | ≤ λ the solution is x̂i = 0. Similarly, when ḡt,i < −λ, then x̂i > 0, theobjective is differentiable, and the solution is obtained by setting the gradient to zero:ηḡt,i +Ht,iiηtx̂i + ηλ = 0 , so that x̂i =(−ḡt,i − λ) .
t
Ht,ii

Likewise, when ḡt,i > λ then x̂i < 0, and the solution is x̂i = Hηtt,ii (−ḡt,i + λ). Combining the threecases, we obtain the simple update (19) for xt+1,i .
G.2 ℓ1 -ball projections
The derivation we give is somewhat terse, and we refer the interested reader to Brucker (1984) orPardalos and Rosen (1990) for more depth. Recall that our original problem (20) is symmetric in its
objective and constraints, so we assume without loss of generality that v  0 (otherwise, we reverse
the sign of each negative component in v, then flip the sign of the corresponding component in thesolution vector). This givesminz1kz − vk22 s.t. ha, zi ≤ c, z  0 .
2

Clearly, if ha, vi ≤ c the optimal z∗ = v, hence we assume that ha, vi > c. We also assume without lossof generality that vi /ai ≥ vi+1 /ai+1 for simplicity of our derivation. (We revisit this assumption at
the end of the derivation.) Introducing Lagrange multipliers θ ∈ R+ for the constraint that ha, zi ≤ cand α ∈ Rd+ for the positivity constraint on z, we get
12

L (z, α, θ) = kz − vk22 + θ(ha, zi − c) − hα, zi .
Computing the gradient of L , we have ∇z L (z, α, θ) = z − v + θa − α. Suppose that we knew the
optimal θ∗ ≥ 0. Using the complementarity conditions on z and α for optimality of z (Boyd andVandenberghe, 2004), we see that the solution z∗i satisfiesvi − θ∗ ai if vi ≥ θ∗ ai∗
zi =0otherwise .
Analogously, the complimentary conditions on ha, zi ≤ c show that given θ∗ , we have


d
d
2 vi
∗
∗
∑ ai [vi − θ ai ]+ = c or ∑ ai ai − θ = c .
+
i=1
i=1
Conversely, had we obtained a value θ ≥ 0 satisfying the above equation, then θ would evidently
induce the optimal z∗ through the equation zi = [vi − θai ]+ .
2155

D UCHI , H AZAN AND S INGER

Now, let ρ be the largest index in {1, . . . , d} such that vi − θ∗ ai > 0 for i ≤ ρ and vi − θ∗ ai ≤ 0
for i > ρ. From the assumption that vi /ai ≤ vi+1 /ai+1 , we have vρ+1 /aρ+1 ≤ θ∗ < vρ /aρ . Thus, had
we known the last non-zero index ρ, we would have obtained


ρ
ρ
vρ
vρ ρ 22 vi∑ ai vi − aρ ∑ ai = ∑ ai ai − aρ < c ,
i=1
i=1
i=1


ρ
ρ
ρ+1
vρ+1
vρ+1
2
2 vi
−
a
v
−
a
=
a
∑ i i aρ+1 ∑ i ∑ i ai aρ+1 ≥ c .
i=1
i=1
i=1
Given ρ satisfying the above inequalities, we can reconstruct the optimal θ∗ by noting that the latter
inequality should equal c exactly when we replace vρ /aρ with θ, that is,ρ∑ ai vi − cθ = i=1ρ 2 .
∑i=1 ai∗(33)

The above derivation results in the following procedure (when ha, vi > c). We sort v in descendρρ−1ing order of vi /ai and find the largest index ρ such that ∑i=1 ai vi − (vρ /aρ ) ∑i=1 a2i < c. We thenreconstruct θ∗ using equality (33) and return the soft-thresholded values of vi (see Algorithm 3). It
is easy to verify that the algorithm can be implemented in O(d log d) time. A randomized searchwith bookkeeping (Pardalos and Rosen, 1990) can be straightforwardly used to derive a linear timealgorithm.

References
J. Abernethy, P. L. Bartlett, A. Rakhlin, and A. Tewari. Optimal strategies and minimax lower
bounds for online convex games. In Proceedings of the Twenty First Annual Conference on
Computational Learning Theory, 2008.
T. Ando. Concavity of certain maps on positive definite matrices and applications to Hadamard
products. Linear Algebra and its Applications, 26:203–241, 1979.
A. Asuncion and D. J. Newman. UCI machine learning repository, 2007. URL http://www.ics.
uci.edu/˜mlearn/MLRepository.html.
P. Auer and C. Gentile. Adaptive and self-confident online learning algorithms. In Proceedings of
the Thirteenth Annual Conference on Computational Learning Theory, 2000.
P. L. Bartlett, E. Hazan, and A. Rakhlin. Adaptive online gradient descent. In Advances in Neural
Information Processing Systems 20, 2007.
A. Beck and M. Teboulle. Mirror descent and nonlinear projected subgradient methods for convex
optimization. Operations Research Letters, 31:167–175, 2003.
J. V. Bondar. Comments on and complements to Inequalities: Theory of Majorization and Its
Applications. Linear Algebra and its Applications, 199:115–129, 1994.
A. Bordes, L. Bottou, and P. Gallinari. Sgd-qn: Careful quasi-newton stochastic gradient descent.
Journal of Machine Learning Research, 10:1737–1754, 2009.
2156

A DAPTIVE S UBGRADIENT M ETHODS

S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
P. Brucker. An O(n) algorithm for quadratic knapsack problems. Operations Research Letters, 3(3):163–166, 1984.
N. Cesa-Bianchi, A. Conconi, and C. Gentile. On the generalization ability of on-line learningalgorithms. IEEE Transactions on Information Theory, 50(9):2050–2057, September 2004.
N. Cesa-Bianchi, A. Conconi, , and C. Gentile. A second-order perceptron algorithm. SIAM Journal
on Computing, 34(3):640–668, 2005.
N. Cesa-Bianchi, Y. Mansour, and G. Stoltz. Improved second-order bounds for prediction withexpert advice. Machine Learning, 66:321–352, 2007.
K. Crammer, O. Dekel, J. Keshet, S. Shalev-Shwartz, and Y. Singer. Online passive aggressivealgorithms. Journal of Machine Learning Research, 7:551–585, 2006.
K. Crammer, M. Dredze, and F. Pereira. Exact convex confidence-weighted learning. In Advancesin Neural Information Processing Systems 22, 2008.
K. Crammer, M. Dredze, and A. Kulesza. Adaptive regularization of weight vectors. In Advancesin Neural Information Processing Systems 23, 2009.
C. Davis. Notions generalizing convexity for functions defined on spaces of matrices. In Proceedings of the Symposia in Pure Mathematics, volume 7, pages 187–201. American MathematicalSociety, 1963.
J. Deng, W. Dong, R. Socher, L. Li, K. Li, and L. Fei-Fei. ImageNet: a large-scale hierarchical image database. In Proceedings of the IEEE Conference on Computer Vision and Pattern
Recognition, 2009.
J. Duchi and Y. Singer. Efficient online and batch learning using forward backward splitting. Journalof Machine Learning Research, 10:2873–2908, 2009.
J. Duchi, S. Shalev-Shwartz, Y. Singer, and A. Tewari. Composite objective mirror descent. In
Proceedings of the Twenty Third Annual Conference on Computational Learning Theory, 2010.
R. Fletcher. A new approach to variable metric algorithms. Computer Journal, 13:317–322, 1970.
D. Grangier and S. Bengio. A discriminative kernel-based model to rank images from text queries.
IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(8):1371–1384, 2008.
E. Hazan and S. Kale. Extracting certainty from uncertainty: regret bounded by variation in costs.
In Proceedings of the Twenty First Annual Conference on Computational Learning Theory, 2008.
E. Hazan, A. Kalai, S. Kale, and A. Agarwal. Logarithmic regret algorithms for online convex
optimization. In Proceedings of the Nineteenth Annual Conference on Computational Learning
Theory, 2006.
J. B. Hiriart-Urruty and C. Lemaréchal. Convex Analysis and Minimization Algorithms II. SpringerVerlag, 1996.
2157

D UCHI , H AZAN AND S INGER

R. A. Horn and C. R. Johnson. Matrix Analysis. Cambridge University Press, 1985.
A. Juditsky, A. Nemirovski, and C. Tauvel. Solving variational inequalities with the stochastic
mirror-prox algorithm. http://arxiv.org/abs/0809.0815, 2008.
A. Kalai and S. Vempala. Efficient algorithms for online decision problems. Journal of Computerand System Sciences, 71(3):291–307, 2003.
G. Lan. An optimal method for stochastic composite optimization. Mathematical Programming
Series A, 2010. Online first; to appear.
D. Lewis, Y. Yang, T. Rose, and F. Li. RCV1: A new benchmark collection for text categorization
research. Journal of Machine Learning Research, 5:361–397, 2004.
H. B. McMahan and M. Streeter. Adaptive bound optimization for online convex optimization. In
Proceedings of the Twenty Third Annual Conference on Computational Learning Theory, 2010.
A. Nedić. Subgradient Methods for Convex Minimization. PhD thesis, Massachusetts Institute ofTechnology, 2002.
A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro. Robust stochastic approximation approach tostochastic programming. SIAM Journal on Optimization, 19(4):1574–1609, 2009.
A. S. Nemirovski and D. B. Yudin. Problem Complexity and Efficiency in Optimization. John Wileyand Sons, 1983.
Y. Nesterov. Smooth minimization of nonsmooth functions. Mathematical Programming, 103:
127–152, 2005.
Y. Nesterov. Primal-dual subgradient methods for convex problems. Mathematical Programming,
120(1):221–259, 2009.
G. Obozinski, B. Taskar, and M. Jordan. Joint covariate selection for grouped classification. Technical Report 743, Dept. of Statistics, University of California Berkeley, 2007.
P. M. Pardalos and J. B. Rosen. An algorithm for a singly constrained class of quadratic programs
subject to upper and lower bounds. Mathematical Programming, 46:321–328, 1990.
A. Rakhlin. Lecture notes on online learning. For the Statistical Machine Learning Course at
University of California, Berkeley, 2009.
G. Salton and C. Buckley. Term weighting approaches in automatic text retrieval. Information
Processing and Management, 24(5), 1988.
N. Z. Shor. Utilization of the operation of space dilation in the minimization of convex functions.
Cybernetics and Systems Analysis, 6(1):7–15, 1972. Translated from Kibernetika.
P. Tseng. On accelerated proximal gradient methods for convex-concave optimization. Technicalreport, Department of Mathematics, University of Washington, 2008.
L. Xiao. Dual averaging methods for regularized stochastic learning and online optimization. Technical Report MSR-TR-2010-23, Microsoft Research, 2010.
2158

A DAPTIVE S UBGRADIENT M ETHODS

M. Zinkevich. Online convex programming and generalized infinitesimal gradient ascent. In Proceedings of the Twentieth International Conference on Machine Learning, 2003.

2159