summaries: Inspired by the analysis in [10], this paper shows bounds on multitask online convex optimization w.r.t. an arbitrary strongly convex regularizer. The task relatedness is expressed I.t.o. task variance, defined by means of a task relatedness matrix. Experiments on real data are presented at the end.

summaries: This paper studies the online multi-task learning setting and proposes a multi-task variant of online mirror descent that works with convex losses, a general class of Bregman divergences, and any positively weighted task-interaction matrix. They provide regret guarantees that, on two different instantiations, improve over the baseline of running independent mirror descent algorithms. The algorithms have straightforward implementations that are used on experiments on four datasets.

summaries: The paper proposes MT-OMD, which allows update sharing among tasks. It is shown that it achieves a regret of O(sqrt{sigma^2(N-1)}sqrt{T}). The paper also extends the OGD and EG and obtains closed-form updates. Numerical results are presented to validate the theoretical findings.

summaries: The paper presents MT-OMD and related variants on multi-task online learning with strongly convex loss functions.  The proposals can improve naive optimization when the task similarity is small measured by the distance of reference vectors. Specifically, it proves a regret bound of the form $sqrt{1+\sigma^2(N-1)}\sqrt{T}$, as opposed to the standard $\sqrt{nT}$. This $\sigma$ measures the task similarity (The smaller $sigma$ is, the closer the tasks are).

summaries: This paper considered multitask online convex optimization (OCO), in which each of $N$ agents learns a possibly different task on a common convex decision set.  The goal of this paper is to minimize multitask regret, defined as the sum of the regret of all agents. A simple way to independently run $N$ instances of OMD  for $N$ agents can achieve a multitask regret bound of $O(\sqrt{NT})$, where $T$ is the time horizon. The main contribution of this paper is to develop multitask OMD (MT-OMD) and establish the regret bound of $O(\sqrt{1+\sigma^2(N-1)}\sqrt{T})$, where $\sigma^2$ is the task variance. If $\sigma<1$, their regret bound is better than that established by the simple way.

abstractText: We introduce and analyze MT-OMD, a multitask generalization of Online Mirror 1 Descent (OMD) which operates by sharing updates between tasks. We prove that 2 the regret of MT-OMD is of order √ 1 + σ2(N − 1) √ T , where σ is the task 3 variance according to the geometry induced by the regularizer, N is the number 4 of tasks, and T is the time horizon. Whenever tasks are similar, that is, σ ≤ 1, 5 this improves upon the √ NT bound obtained by running independent OMDs on 6 each task. Our multitask extensions of Online Gradient Descent and Exponentiated 7 Gradient, two important instances of OMD, are shown to enjoy closed-form updates, 8 making them easy to use in practice. Finally, we provide numerical experiments on 9 several real-world datasets which support our theoretical findings. 10

authors:
Jacob Abernethy
Peter Bartlett
Alexander Rakhlin

title: Multitask learning with expert advice

authors:
Gotz Alefeld
Norbert Schneider

title: On square roots of m-matrices

venue: Linear algebra and its 345 applications,

authors:
Pierre Alquier
Massimiliano Pontil

title: Regret bounds for lifelong learning

venue: In Artificial 347 Intelligence and Statistics,

authors:
Andreas Argyriou
Charles A Micchelli
Massimiliano Pontil
Yiming Ying

title: A spectral 349 regularization framework for multi-task structure learning

venue: Neural Information Processing 350 Systems,

authors:
Bart Bakker
Tom Heskes

venue: Journal of Machine Learning Research,

authors:
Maria-Florina Balcan
Mikhail Khodak
Ameet Talwalkar

title: Provable guarantees for gradient- 354 based meta-learning

authors:
Heinz H Bauschke
Patrick L Combettes

title: Convex analysis and monotone operator theory 357 in Hilbert spaces, volume 408

title: Multitask learning

venue: Machine learning,

authors:
Giovanni Cavallanti
Nicolò Cesa-Bianchi
Claudio Gentile

title: Linear algorithms for online 362 multitask classification

venue: The Journal of Machine Learning Research,

authors:
Nicolò Cesa-Bianchi
Claudio Gentile
Giovanni Zappella

title: A gang of bandits

venue: In Proceed- 364 ings of the 26th International Conference on Neural Information Processing Systems,

authors:
Ofer Dekel
Philip M Long
Yoram Singer

title: Online learning of multiple tasks with a shared 367 loss

authors:
Canh T Dinh
Tung T Vu
Nguyen H Tran
Minh N Dao
Hongyu Zhang

title: Fedu: A unified 372 framework for federated multi-task learning with Laplacian regularization

venue: arXiv preprint 373

authors:
Theodoros Evgeniou
Charles A Micchelli
Massimiliano Pontil
John Shawe-Taylor

title: Learning multiple tasks with kernel methods

venue: Journal of machine learning research,

authors:
Theodoros Evgeniou
Massimiliano Pontil

title: Regularized multi–task learning

venue: In Proceedings 377 of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining,

authors:
Chelsea Finn
Aravind Rajeswaran
Sham Kakade
Sergey Levine

title: Online meta-learning

title: The robustness of the p-norm algorithms

venue: Machine Learning, 53(3):265–299,

authors:
A.S. Georghiades
P.N. Belhumeur
D.J. Kriegman

title: From few to many: illumination cone 384 models for face recognition under variable lighting and pose

venue: IEEE Transactions on Pattern 385 Analysis and Machine Intelligence,

authors:
Adam J Grove
Nick Littlestone
Dale Schuurmans

title: General convergence results for linear 387 discriminant updates

venue: In Proceedings of the tenth annual conference on Computational learning 388 theory,

authors:
Quanquan Gu
Zhenhui Li
Jiawei Han

title: Learning a kernel for multi-task clustering

title: Introduction to online convex optimization

venue: Foundations and Trends in Optimiza-

authors:
David P Helmbold
Robert E Schapire
Yoram Singer
Manfred K Warmuth

title: On-line 394 portfolio selection using multiplicative updates

venue: Mathematical Finance,

authors:
Mark Herbster
Guy Lever

title: Predicting the labelling of a graph via minimum p-seminorm 396 interpolation

authors:
Sham M Kakade
Shai Shalev-Shwartz
Ambuj Tewari

title: Regularization techniques for 399 learning with matrices

authors:
Guangxia Li
Steven CH Hoi
Kuiyu Chang
Wenting Liu
Ramesh Jain

title: Collaborative online 401 multitask learning

authors:
Rui Li
Fenglong Ma
Wenjun Jiang
Jing Gao

title: Online federated multitask learning

authors:
Keerthiram Murugesan
Hanxiao Liu
Jaime Carbonell
Yiming Yang

title: Adaptive smoothed 406 online multi-task learning

title: A modern introduction to online learning

venue: arXiv preprint arXiv:1912.13213,

authors:
Anastasia Pentina
Christoph H Lampert

title: Multi-task learning with labeled and unlabeled 411 tasks

authors:
Gianluigi Pillonetto
Francesco Dinuzzo
Giuseppe De Nicolao

title: Bayesian online multitask 413 learning of gaussian processes

authors:
Avishek Saha
Piyush Rai
Hal Daumé
Suresh Venkatasubramanian

title: Online learning of 416 multiple tasks and their relationships

venue: In Proceedings of the Fourteenth International Conference 417 on Artificial Intelligence and Statistics,

venue: PhD thesis, The 419 Hebrew University of Jerusalem,

title: Graphical multi-task learning

venue: Technical report,

authors:
Changjian Shui
Mahdieh Abbasi
Louis-Émile Robitaille
Boyu Wang
Christian Gagné

title: A 423 principled approach for learning task similarity in multitask learning

venue: In Proceedings of the 28th 424 International Joint Conference on Artificial Intelligence

venue: Data Mining and 426 Knowledge Discovery,

authors:
Jonathan Tuck
Shane Barratt
Stephen Boyd

title: A distributed method for fitting Laplacian 428 regularized stratified models

authors:
Yang Yang
Zhigang Ma
Yi Yang
Feiping Nie
Heng Tao Shen

title: Multitask spectral clustering 430 by exploring intertask correlation

venue: IEEE transactions on cybernetics,

authors:
Chi Zhang
Peilin Zhao
Shuji Hao
Yeng Chai Soh
Bu Sung Lee
Chunyan Miao
Steven CH 432 Hoi

venue: Machine 433 Learning,

authors:
Yu Zhang
Qiang Yang

title: A survey on multi-task learning

venue: IEEE Transactions on Knowledge 435 and Data Engineering,

text:
We introduce and analyze MT-OMD, a multitask generalization of Online Mirror1 Descent (OMD) which operates by sharing updates between tasks. We prove that2 the regret of MT-OMD is of order √ 1 + σ2(N − 1) √ T , where σ2 is the task3 variance according to the geometry induced by the regularizer, N is the number4 of tasks, and T is the time horizon. Whenever tasks are similar, that is, σ2 ≤ 1,5 this improves upon the √ NT bound obtained by running independent OMDs on6
each task. Our multitask extensions of Online Gradient Descent and Exponentiated7 Gradient, two important instances of OMD, are shown to enjoy closed-form updates,8 making them easy to use in practice. Finally, we provide numerical experiments on9 several real-world datasets which support our theoretical findings.10
1 Introduction11
In multitask learning [9], one faces a set of tasks to solve, and tries to leverage their similarities to learn12 faster. Task similarity is often formalized in terms of Euclidean distances among the best performing13 models for each task (see [16] for an example). However, in online convex optimization, and Online14 Mirror Descent (OMD) in particular, it is well known that using different geometries to measure15 distances in the model space can bring substantial advantages — see, e.g., [22, 29]. For instance,16 when the model space is the probability simplex in Rd, running OMD with the KL divergence17 (corresponding to an entropic regularizer) allows to learn at a rate depending only logarithmically18 on d. It is thus natural to investigate to what extent measuring task similarities using geometries that19 are possibly non-Euclidean could improve the analysis of online multitask learning.20
In this work, we introduce MT-OMD, a multitask generalization of OMD which applies to any21 strongly convex regularizer. We present a regret analysis establishing that MT-OMD outperforms22 OMD (run independently on each task) whenever tasks are similar according to the geometry induced23 by the regularizer. Our work is inspired by the multitask extension of the Perceptron algorithm24 developed in [10], where prior knowledge about task similarities is expressed through a symmetric25 positive definite interaction matrix A. Typically, A = I + L, where L is the Laplacian of a task26 relatedness graph with adjacency matrix W . The authors then show that the number of mistakes27 depends on ∑ i ‖ui‖22 + ∑ i,jWij‖ui − uj‖22, where each ui denotes the best model for task i. This28 expression can be seen as a measure of the task dispersion with respect to matrix W and norm ‖ · ‖2.29 The Euclidean norm appears because the Perceptron is an instance of OMD for the hinge loss with the30 Euclidean regularizer, so that distances in the model space are measured through the corresponding31 Bregman divergence, which is the Euclidean squared norm.32
For an arbitrary strongly convex regularizer ψ, the regret of OMD is controlled by a Bregman33 divergence and a term inversely proportional to the curvature of the regularizer. The key challenge34 we face is how to extend the OMD regularizer to the multitask setting so that the dispersion term35 captures task similarities. A natural strategy would be to choose a regularizer whose Bregman36 divergence includes ∑ i,jWijBψ(ui, uj). Although this mimics the Euclidean dispersion term of the37
Submitted to 35th Conference on Neural Information Processing Systems (NeurIPS 2021). Do not distribute.
Perceptron, the associated regularizer has a small curvature, compromising the divergence-curvature38 balance which, as we said, controls the regret. Observing that the Perceptron’s dispersion term can be39 rewritten as ‖A1/2u‖22, where A is a block version (across tasks) of A and u is the concatenation40 of the task vectors ui, our solution consists in using the regularizer ψ ( A1/2 · ) , where ψ is the41 compound version of any base regularizer ψ defined on the model space. While exhibiting the right42 curvature, this regularizer has still the drawback thatA1/2u might be outside the domain of ψ. To43 get around this difficulty, we introduce a notion of variance aligned with the geometry induced by44 ψ, such that the corresponding Bregman divergence Bψ ( A1/2u,A1/2v ) is always defined for sets45 of tasks with small variance. We then show that the Bregman divergence can be upper bounded in46 terms of the task variance σ2, and by tuning appropriately the matrix A we obtain a regret bound for47 MT-OMD that scales as √ 1 + σ2(N − 1), which is approximately σ √ N for σ > 0. In contrast, the48
regret of independent OMD scales as √ N , highlighting the advantage brought by MT-OMD when49 tasks have a small variance. We stress that this improvement is independent of the chosen regularizer,50 thereby offering a substantial acceleration in a wide range of scenarios. The rest of the paper is51 organized as follows. In Section 2, we formally introduce the problem and describe MT-OMD, our52 multitask extension to solve it. In Section 3, we derive a regret analysis for MT-OMD, highlighting53 its advantage when tasks are similar. In Section 4, we detail some algorithmic implementations of54 MT-OMD, while numerical experiments are reported in Section 5.55
Related works. Starting from the seminal work [9], multitask learning has been intensively studied56 for more than two decades, see [40] for a recent survey. Similarly to [10], our work is inspired57 by the Laplacian multitask framework of [15]. This framework has been extended to kernel-based58 learning [34], kernel-based unsupervised learning [21], contextual bandits [11], spectral clustering59 [38], stratified model learning [37], and, more recently, federated learning [14]. See also [24] for60 different applications of Laplacians in online learning. A multitask version of OMD has been61 previously proposed in [25]. Their approach, unlike ours, is cast in terms of matrix learning, and uses62 group norms and Schatten p-norm regularizers. Their bounds scale with the diameter of the model63 space according to these norms (as opposed to scaling with the task variance, as in our analysis).64 Moreover, their learning bias is limited to the choice of the matrix norm regularizer and does not65 explicitly include a notion of task similarity matrix. In [1, 12], different multitask extensions of66 online learning are investigated, see also [3, 17, 6, 13] for related extensions to meta-learning. Some67 online multitask applications are studied in [31, 26, 27], but without providing any regret analyses.68 In [32, 39] authors extend the results of [10] to dynamically updated interaction matrices. However,69 no regret bounds are provided. In [28] they look at distributed online classification in a synchronous70 setting. Other approaches for learning task similarities have been investigated in [41, 30, 35].71
2 Multitask online learning72
In this section, we introduce the problem and describe our approach to solve it. We use a cooperative73 and asynchronous multiagent formalism, where an online algorithm is run in a distributed fashion by74 agents that communicate to solve their tasks, but make predictions at different time steps.75
Problem formulation and reminders on OMD. We consider an online convex optimization76 setting with a set of N ∈ N agents, each learning a possibly different task on a common convex77 decision set V ⊂ Rd. At each time step t = 1, 2, . . . some agent it ≤ N makes a prediction xt ∈ V78 for its task, incurs loss `t(xt), and observes a subgradient of `t at xt, where `t is a convex loss79 function. We say that it is the active agent at time t. Both the sequence i1, i2, . . . of active agents and80 the sequence `1, `2, . . . of convex losses are chosen adversarially and hidden from the agents. Our81 goal is to minimize the multitask regret, defined as the sum of the individual regrets:82
RT = N∑ i=1 ( ∑ t : it=i `t(xt)− inf u∈V ∑ t : it=i `t(u) ) = T∑ t=1 `t(xt)− N∑ i=1 inf u∈V ∑ t : it=i `t(u). (1)
A first natural idea to minimize (1) is to run N independent OMDs, one for each agent. Recall that83 OMD refers to a family of algorithms, typically used to minimize a regret of the form ∑ t `t(xt)−84
infu∈V ∑ t `t(u), for any sequence of proper convex loss functions `t. An instance of OMD is85 parameterized by a λ-strongly convex regularizer ψ : Rd → R, and has the update rule86 xt+1 = argmin
x∈V 〈ηtgt, x〉+Bψ(x, xt), (2)
where gt ∈ Rd is a subgradient of `t at point xt, Bψ(x, y) = ψ(x)−ψ(y)−〈∇ψ(y), x− y〉 denotes87 the Bregman divergence associated to ψ, and ηt > 0 is a tunable learning rate. Standard results allow88 to bound the regret achieved by the sequence of iterates produced by OMD. For example, for a fixed89 learning rate η and for any initial point x1 ∈ V , we have [29, Theorem 6.8]:90
∀u ∈ V, T∑ t=1 `t(xt)− `t(u) ≤ Bψ(u, x1) η + η 2λ T∑ t=1 ‖gt‖2?, (3)
with ‖ · ‖? the dual norm of the norm with respect to which ψ is strongly convex (see Definition 4).91 The choice of the regularizer ψ shapes the above bound through the quantities Bψ(u, x1) and ‖gt‖?.92 When ψ = 12‖ · ‖ 2 2, we have Bψ(x, y) = 1 2‖x − y‖ 2 2, ‖ · ‖? = ‖ · ‖2, λ = 1, and the algorithm93 is known as Online Gradient Descent (OGD). However, depending on the problem, a different94 choice of the regularizer might be better suited to the underlying geometry. A well-known example95 is Exponentiated Gradient (EG), an instance of OMD in which V is the probability simplex in96 Rd, such that V = ∆ := {x ∈ Rd+ : ∑ j xj = 1}. EG uses the negative entropy regularizer97
x 7→ ∑ j xj ln(xj), and assuming that ‖gt‖∞ ≤ L, one achieves bounds of order O(L √ T ln d),98
while OGD yields bounds of orderO(L √ dT ). We emphasize that our cooperative extension adapts to99 several types of regularizers, and can thus exploit these improvements with respect to the dependence100 on d, see Proposition 5. Let C be the generic constant such that C √ T upper bounds the regret101
incurred by the chosen OMD, e.g., C = L √ ln d, or C = L √ d in the above examples. Using Jensen’s102 inequality, one can show that the regret achieved by N independent OMDs satisfies:103
RT ≤ N∑ i=1 C √ Ti ≤ C √ NT, (4)
where Ti = ∑T t=1 I{it = i} denotes the number of times agent i was active. Our goal is to show that104 introducing communication between the agents may significantly improve on (4) with respect to the105 dependence on N . We now describe the multitask extension we propose.106
A multitask extension. In order to build some intuition about our approach, we first focus on OGD.107 For i ≤ N and t ≤ T , let xi,t ∈ Rd denote the prediction maintained by agent i at time step t. By108 completing the square in Equation (2) for ψ = ψEuc := 12‖ · ‖ 2 2, the independent OGDs updates can109 be rewritten for all i ≤ N and t such that it = i:110
xi,t+1 = ΠV,‖·‖2 ( xi,t − ηtgt ) , (5)
where ΠV,‖·‖ denotes the projection operator onto the convex set V according to the norm ‖ · ‖, such111 that ΠV,‖·‖(x) = argminy∈V ‖x− y‖. Our analysis next relies on compound representations, that112 we explain next. We use bold notation to refer to compound vectors, such that for u1, . . . , uN ∈ Rd,113 the compound vector is u = [u1, . . . , uN ] ∈ RNd. For i ≤ N , we use u(i) to refer to the ith block of114 u, such that u(i) = ui in the above example. So xt is the compound vector of the (xi,t)Ni=1, such that115 xt = x (it) t , and the multitask regret rewrites as RT (u) = ∑T t=1 `t ( x (it) t ) − `t ( u(it) ) . For any set116 V ⊂ Rd, let V = V ⊗N ⊂ RNd denote the compound set such that u ∈ V is equivalent to u(i) ∈ V117 for all i ≤ N . Equipped with this notation, the independent OGD updates (5) rewrite as:118
xt+1 = ΠV ,‖·‖2 ( xt − ηtḡt ) , (6)
with ḡt ∈ RNd such that ḡ(i)t = gt for i = it, and 0Rd otherwise. In other words, only the active agent119 has a non-zero gradient and therefore makes an update. Our goal is to incorporate communication120 into this independent update. To that end, we consider the general idea of sharing updates by121 considering gradients of the formA−1ḡt, whereA−1 ∈ RNd×Nd is a shortcut notation for A−1⊗ Id122 and A ∈ RN×N is any symmetric positive definite interaction matrix. Note that A is a parameter of123 the algorithm playing the role of a learning bias. While our main result (Theorem 1) holds for any124 choice of A, our more specialized bounds (Propositions 2 to 5) apply to a parameterized family of125 matrices A. A simple computation shows that (A−1ḡt)(i) = A−1iit gt. Thus, every agent i makes an126 update proportional to A−1iit at each time step t. In other words, the active agent (the only one to suffer127 a loss) shares its update with the other agents. The crux of the proofs consists in finding a matrix128
A−1 (or equivalently A) such that A−1iit captures the similarity between tasks i and it. Intuitively, an129 agent i is interested in sharing the update (gradient) of agent it to the extent their respective tasks are130 similar. Overall, the MT-OGD update (our cooperative extension of OGD) thus writes:131
xt+1 = ΠV ,‖·‖A ( xt − ηtA−1ḡt ) , (7)
where ‖ · ‖M denotes the Mahalanobis norm such that ‖u‖M = √ u>Mu. In comparison to (6), the132 need for changing the norm in the projection, although unclear at first sight, can be explained in133 multiple ways. First, it is key to the analysis, as we see by examining the proof of Theorem 1. Second,134 it can be interpreted as another way of exchanging information between agents, see Remark 1. Finally,135 observe that update (7) can be decomposed as:136
x̃t+1 = argmin x∈RNd
〈ηtḡt,x〉+ 1
2 ‖x− xt‖2A,
xt+1 = argmin x∈V
1 2 ‖x− x̃t+1‖2A,
(8)
showing that it is natural to keep the same norm in both updates. Most importantly, what update (8)137 tells us, is that the MT-OGD update rule is actually an OMD update — see, e.g., [29, Section 6.4]138 — with regularizer x 7→ 12‖x‖ 2 A = ψEuc ( A1/2x ) . This provides a natural path for extending our139 multitask approach to any regularizer. Given a base regularizer ψ : Rd → R, the compound regularizer140 ψ is given byψ : x ∈ RNd 7→ ∑N i=1 ψ ( x(i) ) . For base regularizers for which there exists φ : R→ R141
such that ψ(x) = ∑ j φ(xj), the compound regularizer is the natural extension of ψ to RNd. Note,142 however, that the relationship can be more complex, e.g., when the regularizer is the squared p-norm143 x 7→ 12‖x‖ 2 p. Now, using regularizer ψ ( A1/2 · ) , whose Bregman writes Bψ ( A1/2x,A1/2x′ ) , the144 MT-OMD update reads:145 xt+1 = argmin
x∈V 〈ηtḡt,x〉+Bψ
( A1/2x,A1/2xt ) . (9)
Clearly, when ψ = ψEuc, we recover the MT-OGD update. Observe also that whenever A = IN ,146 MT-OMD is equivalent to N independent OMDs. We conclude this exposition with an important147 remark shedding light on the way MT-OMD introduces communication between agents.148 Remark 1. Using the substitution yt = A1/2xt, update (9) can be rewritten as:149
xt+1 = A −1/2 argmin y∈A1/2(V ) 〈ηtA−1/2ḡt,y〉+Bψ(y,yt). (10)
The two occurrences ofA−1/2 reveal that agents communicate in two distinct ways: one through the150 shared update (the innermost occurrence ofA−1/2), and one through computing the final prediction151 xt+1 as a linear combination of the solution to the optimization problem. Multiplying both sides of152 Equation (10) byA1/2, MT-OMD can also be seen as standard OMD on the transformed iterate yt.153
3 Regret analysis154
We now provide a regret analysis for MT-OMD. We start with a general theorem presenting two155 bounds, for constant and time-varying learning rates. These results are then instantiated to different156 types of regularizer and variance in Propositions 1 to 5. The main difficulty in the proof is to157 characterize the strong convexity of ψ ( A1/2 · ) , see Lemmas 1 and 2. Note that the set V is always158 assumed to be convex. All technical proofs can be found in Appendix A.159 Theorem 1. Let V ⊂ Rd, ψ : Rd → R be λ-strongly convex with respect to norm ‖ · ‖ on V , and160 A ∈ RN×N be symmetric positive definite. Let (`t)Tt=1 be a sequence of proper convex loss functions,161 and set x1 ∈ V . Then, MT-OMD with ηt := η produces a sequence of iterates (xt)Tt=1 which satisfies162 for all u ∈ V :163
RT (u) ≤ Bψ ( A1/2u,A1/2x1 ) η + max i≤N A−1ii η 2λ T∑ t=1 ‖gt‖2?. (11)
Moreover, for any sequence of nonincreasing learning rates (ηt)Tt=1, MT-OMD produces a sequence164 of iterates (xt)Tt=1 which satisfies for all u ∈ V :165
RT (u) ≤ max t≤T
Bψ ( A1/2u,A1/2xt ) ηT + max i≤N A−1ii 1 2λ T∑ t=1 ηt‖gt‖2?. (12)
3.1 Multitask Online Gradient Descent (MT-OGD)166
For ψ = 12‖ · ‖ 2 2, A = IN (independent updates), unit-norm reference vectors (u (i))Ni=1, L-Lipschitz167 losses, and x1 = 0, bound (11) becomes: ND2/2η + ηTL2/2. Choosing η = D √ N/L √ T , we168
recover the DL √ NT bound of Equation (4). Our goal is to design interaction matrices A that make169 the right-hand side of (11) smaller. In absence of additional assumptions on the set of comparators,170 it is however impossible to get a systematic improvement: the bounds are the sums of two terms,171 and introducing interactions typically reduces one term but increases the other. To get around this172 difficulty, we introduce a simple condition on the task similarity, that allows to control the increase of173 Bψ ( A1/2u,A1/2x1 ) for a carefully designed class of interaction matrices.174
Definition 1. Let ‖ · ‖ : Rd → R be any norm. We define the variance of u ∈ RNd w.r.t. ‖ · ‖ as:175
Var‖·‖(u) = 1
N − 1 N∑ i=1 ∥∥u(i) − ū∥∥2, where ū = 1 N N∑ i=1 u(i).
Let V ⊂ Rd, D = supu∈V ‖u‖, and σ > 0. The set of comparators with variance smaller than176 σ2D2 is denoted by:177 V‖·‖,σ = { u ∈ V : Var‖·‖(u) ≤ σ2D2 } . (13)
For sets of comparators of the form (13), we show that MT-OGD achieves significant improvements178 over its independent counterpart. The rationale behind this gain is fairly natural: the tasks associated179 with comparators in (13) are similar due to the variance constraint, so that communication may indeed180 help. Note that this variance condition fundamentally involves information about the entire compound181 reference vector u, and is consequently irrelevant when OGD instances are run independently. The182 class of communication matrices we consider is parameterized as follows. Define L = IN − 11>/N .183 We consider matrices of the form A(b) = IN + b L, where b ≥ 0 quantifies the communication184 magnitude. For more intuition about this choice, we refer the reader to Remark 2. We can now state a185 first result highlighting the advantage of MT-OGD.186 Proposition 1. Let V ⊂ Rd, D = supx∈V ‖x‖2, ψ = 12‖ · ‖ 2 2, and σ ≤ 1. Let (`t)Tt=1 be a sequence187 of proper convex loss functions such that ‖gt‖2 ≤ L for all gt ∈ ∂`t(x), x ∈ V . Set b = N , x1 = 0,188 and η = D √ N(N + 1)(1 + (N − 1)σ2)/L √ 2T . Then, MT-OGD produces a sequence of iterates189
(xt) T t=1 which satisfies for all u ∈ V‖·‖2,σ:190 RT (u) ≤ DL √ 1 + σ2(N − 1) √ 2T . (14)
Proof sketch. For the particular choice ψ = 12‖ · ‖ 2 2, x1 = 0, and a generic b, one can show that191 2Bψ ( A(b)1/2u,A(b)1/20 ) = ‖u‖22 + b(N − 1)Var‖·‖2(u). Furthermore, it is easy to check that192 [
A(b)−1 ] ii = b+N(1+b)N for all i ≤ N . Substituting these values into (11), we get the upper bound:193
‖u‖22 + b(N − 1)Var‖·‖2(u) 2η + ηTL2 2 b+N (1 + b)N ≤ ND2(1 + bN−1N σ 2) 2η + ηTL2 2 b+N (1 + b)N .
Setting η = NDL
√( 1+bN−1N σ 2 ) (1+b)
(b+N)T , and then b = N , yields the desired bound.194 Thus, MT-OGD enjoys a √ 1 + σ2(N − 1) dependence, which is smaller than √ N whenever tasks195 have a variance smaller than 1. Considering limit cases is also informative. When σ = 0 (all tasks are196 equal), MT-OGD scales as if there were only one task. When σ ≥ 1, the analysis suggests to choose197 b = 0, i.e., A = IN , and one recovers the performance of independent OGDs. Note that the additional198 √
2 factor can be removed for limit cases through a better optimization in b: the bound obtained is199 DL √ C(σ)T , with C(0) = 1 and C(1) = N . However, the function C lacks of interpretability200 outside of the limit cases (see Appendix A.2) motivating our choice to present the looser but more201 interpretable bound (14). For large values of N and σ > 0, we have √ 1 + σ2(N − 1) ≈ σ √ N . The202 improvement brought by MT-OGD is thus roughly proportional to the square root of the task variance.203 From now on, we refer to this gain as the variance acceleration.204
Another way to gain intuition about this upper bound is to compare it to the independent OGD lower205 bound. Consider the following setting: let N = 2d, u0 ∈ Rd be such that ‖u0‖22 = 1 − σ2, and206
set for all i ≤ d: u(2i−1) = u0 − σ′ei, and u(2i) = u0 + σ′ei, where σ′ = σ √
(N − 1)/N , and207 (ei)i≤d is an orthonormal basis of Rd. It is easy to check that ‖u(i)‖22 ≥ 1 − 2σ2 for all i ≤ N ,208 and that Var‖·‖2(u) = σ
2. Then, a standard lower bound for OGD [29, Theorem 5.1] applied209 to the independent instances gives RT ≥ ∑ i ‖u(i)‖2 √ 2Ti/4 ≥ √ 2(1− 2σ2)NT/4, assuming210 unit-norm loss gradients and Ti = T/N . This lower bound is strictly greater than the MT-OGD upper211 bound (14) as soon as √ (1− 2σ2)N/4 ≥ √ 1 + σ2N , or equivalently σ2 ≤ (N − 16)/(18N).212
3.2 Extension to other norm regularizers213
Now, a natural question is: can the variance acceleration be achieved with other regularizers? Indeed,214 the proof of Proposition 1 crucially relies on the fact that the Bregman term can be exactly expressed215 in terms of ‖u‖22 and Var‖·‖2(u). In the following proposition, we show that such an improvement216 is also possible for all regularizers of the form 12‖ · ‖
2, for arbitrary norms ‖ · ‖, up to an additional217 multiplicative constant. A crucial application is the use of the p-norm on the probability simplex,218 which is known to exhibit a logarithmic dependence with respect to d when p = (2 ln d) / (2 ln d− 1).219
Proposition 2. Let ‖·‖ : Rd → R be any norm, V ⊂ Rd, D = supx∈D ‖x‖, ψ = 12‖·‖ 2, and σ ≤ 1.220 Let (`t)Tt=1 be a sequence of proper convex loss functions such that ‖gt‖? ≤ L for all gt ∈ ∂`t(x),221 x ∈ V . Set b = N , x1 = 0, and η = D √ N(N + 1)(1 + (N − 1)σ2)/L √ 2T . Then, MT-OMD222 produces a sequence of iterates (xt)Tt=1 which satisfies for all u ∈ V‖·‖,σ:223
RT (u) ≤ DL √ 1 + σ2(N − 1) √ 8T .
In particular, for d ≥ 3 and V = ∆, choosing ‖ ·‖ = ‖ ·‖p, for p = 2 ln d/(2 ln d−1), and assuming224 that ‖gt‖∞ ≤ L, it holds for all u ∈∆‖·‖p,σ:225
RT (u) ≤ L √ 1 + σ2(N − 1) √ 16e T ln d.
In comparison, under the same assumptions, bound (14) would write as: L √ 1 + σ2(N − 1) √ 2Td.226
Projecting onto Vσ . The previous results reveal that whenever tasks are similar (i.e., u ∈ Vσ),227 then using the regularizer ψ ( A1/2 · ) introduces communication and accelerates the convergence.228 However, introducing communication is not the only way to leverage the small variance condition.229 For instance, one may also use this information to directly project onto Vσ , by considering the update:230
231
xt+1 = argmin x∈Vσ
〈ηtḡt,x〉+Bψ ( A1/2x,A1/2xt ) . (15)
Although not necessary in general (Propositions 1 and 2 show that gradient communication only is232 sufficient to get an improvement), this refinement presents several advantages. First, it might actually233 be simpler to compute in practice, see Section 4. Second, it allows for adaptive learning rates, that234 preserve the guarantees while being independent from the horizon T (Proposition 3). Finally, it235 allows to derive L? bounds with the variance acceleration for smooth loss functions (Proposition 4).236 The results are stated for arbitrary norms, but bounds sharper by a factor 2 can be obtained for ‖ · ‖2.237
Proposition 3. Let ‖·‖ : Rd → R be any norm, V ⊂ Rd, D = supx∈D ‖x‖, ψ = 12‖·‖ 2, and σ ≤ 1.238 Let (`t)Tt=1 be a sequence of proper convex loss functions. Set ηt = D √ N(N + 1)(1 + (N − 1)σ2)239
/ √∑t
i=1 ‖gi‖2?, and b = N . Then, update (15) produces a sequence of iterates (xt)Tt=1 which240 satisfies for all u ∈ V‖·‖,σ:241
RT (u) ≤ 8D √ 1 + σ2(N − 1) √√√√ T∑ t=1 ‖gt‖2?.
Proposition 4. Let ‖·‖ : Rd → R be any norm, V ⊂ Rd, D = supx∈D ‖x‖, ψ = 12‖·‖ 2, and σ ≤ 1.242 Let (`t)Tt=1 be a sequence of convexM -smooth loss functions, i.e., ‖∇`t(x)−∇`t(y)‖? ≤M‖x−y‖243 for all t ≤ T , x, y ∈ V . Set b and ηt as in Proposition 3. Then, update (15) produces a sequence of244 iterates (xt)Tt=1 which satisfies for all u ∈ V‖·‖,σ:245
RT (u) ≤ 16D √ 1 + σ2(N − 1) 2MD√1 + σ2(N − 1) + √√√√M T∑
t=1
`t ( u(it) ) .
3.3 Regularizers on the simplex246
As can be seen from Propositions 1 to 4, MT-OMD induces a variance acceleration for several247 regularizers (Euclidean norm, p-norms), and under different assumptions on the losses (Lipschitz248 continuity, smoothness). This improvement is thus orthogonal to the one achievable by choosing the249 regularizer appropriately, and bounds benefiting from both improvements can be derived, such as the250 second claim of Proposition 2. However, all regularizers studied so far share a crucial feature: they251 are defined on the entire space Rd. As a consequence, the Bregman divergence Bψ(A1/2u,A1/2x)252 is always well defined, which might not be true in general. Of particular interest is the case where253 regularizers are defined on the probability simplex ∆. A workaround consists in assigning the value254 +∞ to the Bregman whenever either of the arguments is outside of the compound simplex ∆ = ∆⊗N .255 The choice of the interaction matrixA then becomes critical to prevent the bound from exploding, and256 calls for a new definition of the variance. Indeed, for i ≤ N we have (A(b)1/2)(i) = (1+b)u(i)−bū.257 If all u(i) are equal (say to u0 ∈ ∆), then all (A(b)1/2u)(i) are also equal to u0 andA(b)1/2u ∈∆.258 However, if they are different, by definition of ū, for all j ≤ d, there exists i ≤ N such that259 u
(i) j ≤ ūj . Then, for b large enough, (1 + b)u (i) j − būj becomes negative, and (A(b)1/2u)(i) is260 out of the simplex. Hopefully, the maximum acceptable value for b can be easily deduced from the261 following variance definition.262
Definition 2. Let u ∈ Rd. For all j ≤ d, let umaxj = maxi≤N u (i) j and u min j = mini≤N u (i) j . Then,263 given σ ≤ 1 and with the convention 0/0 = 0 we define:264
Var∆(u) = max j≤d
( umaxj − uminj
umaxj
)2 , and ∆σ = { u ∈∆ : Var∆(u) ≤ σ2 } .
Equipped with this new variance definition, we can now analyze regularizers defined on the simplex.265 Proposition 5. Let σ ≤ 1, and ψ : ∆→ R be λ-strongly convex w.r.t. norm ‖ · ‖. Assume that there266 exist x∗ ∈ ∆ and C<+∞ such that for all x ∈ ∆, Bψ(x, x∗) ≤ C. Let (`t)Tt=1 be a sequence of267 proper convex loss functions such that ‖gt‖?≤L for all gt ∈ ∂`t(x), x ∈ ∆. Set b = (1− σ2)/σ2,268 x1 = [x ∗, . . . , x∗], and η = N √ 2λ(1 + b)C/L √
(b+N)T . Then, MT-OMD produces a sequence269 of iterates (xt)Tt=1 which satisfies for all u ∈∆σ:270
RT (u) ≤ L √ 1 + σ2(N − 1) √ 2CT/λ.
In particular, for the negative entropy we have x∗=1/d and C= ln d. Hence, for subgradients271 satisfying ‖gt‖∞ ≤ L we obtain:272
RT (u) ≤ L √ 1 + σ2(N − 1) √ 2T ln d.
Proposition 5 shows that if the variance definition is adapted to the geometry of the problem, then the273 variance acceleration can be extended to many regularizers. We conclude the theoretical analysis by274 discussing more local variance definitions and different classes of interaction matrices.275 Remark 2. Observe that we can rewrite Var‖·‖2(u) as:276
(N − 1)Var‖·‖2(u) = ∑ i ‖u(i) − ū‖22 = 1 N ∑ i,j ‖u(i) − u(j)‖22 = u>Lu,
where L = IN − 11>/N is the Laplacian matrix of the weighted clique graph over {1, . . . , N},277 with edges of 1/N . Therefore, a natural extension consists in considering variances of the form:278 VarW‖·‖2(u) = ∑N i,j=1Wij‖u(i)−u(j)‖22 = u>LWu, for any adjacency matrixW and its Laplacian279 LW . For instance, if we expect tasks to be concentrated in clusters, it is natural to consider Wij = 1280 if u(i) and u(j) (are thought to) belong to the same cluster, and 0 otherwise. This local version is281 interesting as it allows to satisfy the variance condition with a smaller σ, which improves the MT-OMD282 regret bound. The proof of Theorem 1 can be readily adapted to this definition by considering the class283 of interaction matrices {A(b) = IN + bLW }. The complete analysis however requires computing284 maxi≤N [A(b)
−1]ii, which depends on W in a nontrivial way, preventing from stating a general285 result for an arbitrary W . Note that local variances can be similarly considered on the simplex.286 Instead of involving the global umaxj , the variance formula then features for each task/node i a local287
maximum (respectively minimum) taken over the neighbours N(i) of i: umax,ij = maxk∈N(i) u (k) j .288
4 Algorithmic implementation289
We now show that MT-OGD and MT-EG enjoy closed-form updates, making them easy to implement.290
Computing MT-OGD. Consider V = {u ∈ Rd : ‖u‖2 ≤ D} and σ ≤ 1. Recall that we have V = V ⊗N = {u ∈ RNd : ‖u‖2,∞ ≤ D}, and V‖·‖2,σ = {u ∈ V : Var‖·‖2(u) ≤ σ}. Solving the first equation in (8), the iterate xt+1 produced by MT-OGD is the solution to:
min x∈RNd {∥∥xt − ηtA−1ḡt − x∥∥2A : ‖x‖2,∞ ≤ D}. However, computing this update is made difficult by the discrepancy between the norms used in291 the objective and the constraint. A simple solution consists in projecting onto the Mahalanobis292 ball VA = {u ∈ RNd : ‖u‖2A ≤ (1 + bσ2)ND2}. It is easy to check that V‖·‖2,σ ⊂ VA, so that293 every result derived in Section 3 for MT-OGD remains valid (only the fact that comparators and294 iterates are in VA is actually used). Then, using the substitution yt = A1/2xt and the notation295 Proj(x, τ) = min { 1, τ‖x‖2 } x, the MT-OGD update rewrites (see Appendix B.1 for details):296
yt+1 = Proj ( yt − ηtA−1/2ḡt, √ (1 + bσ2)ND ) . (16)
Note that (16) can be easily turned back into an update on the xt’s by making the inverse substitution.297 However, (16) is closer to practice, xt being actually only computed to make predictions.298
Computing MT-EG. Using the two-step formulation (8) with yt = A1/2xt, MT-EG is given by:299
ỹt+1 = argmin y∈RNd
〈ηtA−1/2gt,y〉+Bψ(y,yt), (17)
yt+1 = argmin y∈A1/2(∆) Bψ(y, ỹt+1), (18)
where ψ is the compound negative entropy regularizer such that ψ(x) = ∑N i=1 ∑d j=1 x (i) j ln ( x (i) j ) .300 One can show (see Appendix B.2 for details) that the update can be rewritten as:301
∀i ≤ N, j ≤ d,  ỹ (i) t+1,j = y (i) t,j exp ( −ηtA−1/2iit gt,j − 1 ) , y (i) t+1,j = ỹ (i) t+1,j∑d
k=1 ỹ (i) t+1,k
.
Combining both equations, we finally obtain that302
∀i ≤ N, j ≤ d, y(i)t+1,j = y
(i) t,j e −ηtA−1/2iit gt,j∑d k=1 y (i) t,k e −ηtA−1/2iit gt,k . (19)
Update (19) enjoys a natural interpretation. Each block y(i) is operating an individual standard EG303 update, but with gradient A−1/2iit gt. When A = IN , only the active block is updated. Otherwise, the304 update of block i is proportional to A−1/2iit , which quantifies the similarity between tasks i and it.305 Note that the derivation of (19) works for any matrix A, as long as A−1/2 is stochastic. This is true306 in particular for matrices of the form A = IN − LW , see Appendix B.2 for technical details.307
5 Experiments308
In this section, we empirically compare the cumulative performance of MT-OGD and MT-EG against309 their task-independent versions, I-OGD and I-EG. All experiments show that MT-OMD outperforms310 I-OMD for all values of T in the range of the data, see Figure 1.311
OGD. We run OGD using the parameterized clique A(b) = (1 + b)IN − 11>/N and the square312 loss `t(w) = ( w>xt − yt )2 . The paramenters b (MT-OGD) and η (MT-OGD and I-OGD) were set313 according to their optimal theoretical values. The values D = 10 (diameter) and σ = 0.01 (variance),314 which are required to tune η in MT-OGD, were found through grid-search (oracle tuning) and then315
used in all experiments. We use the two following datasets. Lenk [4] consists of 2880 computer316 ratings made by 180 individuals (the tasks) on the basis of 14 binary features. Each computer is rated317 on a discrete scale from 0 to 10, expressing the likelihood of an individual buying that computer.318 School [5] consists of 15362 student performances based on 28 features and collected from 139319 schools representing the tasks. Features and labels are in the range [0, 1]. Both plots are averaged320 over 3 independent runs with respect to the random choice of the initial vector for OGD (note that321 Lenk is simpler than School, hence the lack of variance in the plot).322
EG. For EG we considered the two following datasets. Extended Yale B [19], a face classification323 dataset containing 64 images (1024 grey levels) for each one of 38 individuals. We took a subset324 consisting of the first 15 individuals and ran random forest to select 10 features. We generated325 (
15 2 ) = 105 binary classification tasks via the all-pairs reduction and used the logistic loss `t(w) =326
log ( 1 + e−ytw >xt ) . The interaction matrix A used by MT-EG is a regular graph, connecting pairs of327 tasks sharing the same person. We set η in MT-EG using oracle tuning, whereas η in I-EG was set328 according to its optimal theoretical value. NYSE [23] is a portfolio selection dataset consisting of329 5651 price relatives for 36 stocks. The goal at each round is to select a portfolio (i.e., a probability330 vector w). We created 100 tasks by assigning uniformly at random the data to the tasks. The reward331 incurred at each round is the log wealth log(w>xt). We ran MT-EG using A = IN + L, where L is332 the Laplacian of the clique. The value of η in MT-EG was set using oracle tuning, whereas for I-EG333 we chose the value η = 0.05 suggested in [23]. The plot is averaged over 3 independent runs with334 respect to the random assignment of data to the tasks.335
6 Conclusion336
We introduced and analyzed MT-OMD, a multitask extension of OMD whose regret is shown to337 improve as the task variance, expressed in terms of the geometry induced by the regularizer, decreases.338 Natural and interesting directions for future research are: (1) extending the variance acceleration to339 strongly convex losses; (2) proving similar regret bounds without resorting to the knowledge of σ2.340
References341 [1] Jacob Abernethy, Peter Bartlett, and Alexander Rakhlin. Multitask learning with expert advice.342 In International Conference on Computational Learning Theory, pages 484–498. Springer,343 2007.344
[2] Gotz Alefeld and Norbert Schneider. On square roots of m-matrices. Linear algebra and its345 applications, 42:119–132, 1982.346
[3] Pierre Alquier, Massimiliano Pontil, et al. Regret bounds for lifelong learning. In Artificial347 Intelligence and Statistics, pages 261–269. PMLR, 2017.348
[4] Andreas Argyriou, Charles A Micchelli, Massimiliano Pontil, and Yiming Ying. A spectral349 regularization framework for multi-task structure learning. Neural Information Processing350 Systems, pages 25–32, 2007.351
[5] Bart Bakker and Tom Heskes. Task clustering and gating for bayesian multitask learning.352 Journal of Machine Learning Research, 4:83–99, 2003.353
[6] Maria-Florina Balcan, Mikhail Khodak, and Ameet Talwalkar. Provable guarantees for gradient-354 based meta-learning. In International Conference on Machine Learning, pages 424–433. PMLR,355 2019.356
[7] Heinz H Bauschke, Patrick L Combettes, et al. Convex analysis and monotone operator theory357 in Hilbert spaces, volume 408. Springer, 2011.358
[8] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge359 university press, 2004.360
[9] Rich Caruana. Multitask learning. Machine learning, 28(1):41–75, 1997.361
[10] Giovanni Cavallanti, Nicolò Cesa-Bianchi, and Claudio Gentile. Linear algorithms for online362 multitask classification. The Journal of Machine Learning Research, 11:2901–2934, 2010.363
[11] Nicolò Cesa-Bianchi, Claudio Gentile, and Giovanni Zappella. A gang of bandits. In Proceed-364 ings of the 26th International Conference on Neural Information Processing Systems, pages365 737–745, 2013.366
[12] Ofer Dekel, Philip M Long, and Yoram Singer. Online learning of multiple tasks with a shared367 loss. Journal of Machine Learning Research, 8(10), 2007.368
[13] Giulia Denevi, Dimitris Stamos, Carlo Ciliberto, and Massimiliano Pontil. Online-within-369 online meta-learning. In Advances in Neural Information Processing Systems 32 (NIPS 2019),370 volume 32, pages 1–11, 2019.371
[14] Canh T Dinh, Tung T Vu, Nguyen H Tran, Minh N Dao, and Hongyu Zhang. Fedu: A unified372 framework for federated multi-task learning with Laplacian regularization. arXiv preprint373 arXiv:2102.07148, 2021.374
[15] Theodoros Evgeniou, Charles A Micchelli, Massimiliano Pontil, and John Shawe-Taylor.375 Learning multiple tasks with kernel methods. Journal of machine learning research, 6(4), 2005.376
[16] Theodoros Evgeniou and Massimiliano Pontil. Regularized multi–task learning. In Proceedings377 of the tenth ACM SIGKDD international conference on Knowledge discovery and data mining,378 pages 109–117, 2004.379
[17] Chelsea Finn, Aravind Rajeswaran, Sham Kakade, and Sergey Levine. Online meta-learning.380 In International Conference on Machine Learning, pages 1920–1930. PMLR, 2019.381
[18] Claudio Gentile. The robustness of the p-norm algorithms. Machine Learning, 53(3):265–299,382 2003.383
[19] A.S. Georghiades, P.N. Belhumeur, and D.J. Kriegman. From few to many: illumination cone384 models for face recognition under variable lighting and pose. IEEE Transactions on Pattern385 Analysis and Machine Intelligence, 23(6):643–660, 2001.386
[20] Adam J Grove, Nick Littlestone, and Dale Schuurmans. General convergence results for linear387 discriminant updates. In Proceedings of the tenth annual conference on Computational learning388 theory, pages 171–183, 1997.389
[21] Quanquan Gu, Zhenhui Li, and Jiawei Han. Learning a kernel for multi-task clustering. In390 Proceedings of the AAAI Conference on Artificial Intelligence, volume 25, 2011.391
[22] Elad Hazan. Introduction to online convex optimization. Foundations and Trends in Optimiza-392 tion, 2(3-4):157–325, 2016.393
[23] David P Helmbold, Robert E Schapire, Yoram Singer, and Manfred K Warmuth. On-line394 portfolio selection using multiplicative updates. Mathematical Finance, 8(4):325–347, 1998.395
[24] Mark Herbster and Guy Lever. Predicting the labelling of a graph via minimum p-seminorm396 interpolation. In NIPS Workshop 2010: Networks Across Disciplines: Theory and Applications,397 2009.398
[25] Sham M Kakade, Shai Shalev-Shwartz, and Ambuj Tewari. Regularization techniques for399 learning with matrices. The Journal of Machine Learning Research, 13(1):1865–1890, 2012.400
[26] Guangxia Li, Steven CH Hoi, Kuiyu Chang, Wenting Liu, and Ramesh Jain. Collaborative online401 multitask learning. IEEE Transactions on Knowledge and Data Engineering, 26(8):1866–1876,402 2014.403
[27] Rui Li, Fenglong Ma, Wenjun Jiang, and Jing Gao. Online federated multitask learning. In404 2019 IEEE International Conference on Big Data (Big Data), pages 215–220. IEEE, 2019.405
[28] Keerthiram Murugesan, Hanxiao Liu, Jaime Carbonell, and Yiming Yang. Adaptive smoothed406 online multi-task learning. In Advances in Neural Information Processing Systems, pages407 4296–4304, 2016.408
[29] Francesco Orabona. A modern introduction to online learning. arXiv preprint arXiv:1912.13213,409 2019.410
[30] Anastasia Pentina and Christoph H Lampert. Multi-task learning with labeled and unlabeled411 tasks. In International Conference on Machine Learning, pages 2807–2816. PMLR, 2017.412
[31] Gianluigi Pillonetto, Francesco Dinuzzo, and Giuseppe De Nicolao. Bayesian online multitask413 learning of gaussian processes. IEEE Transactions on Pattern Analysis and Machine Intelligence,414 32(2):193–205, 2008.415
[32] Avishek Saha, Piyush Rai, Hal Daumé, Suresh Venkatasubramanian, et al. Online learning of416 multiple tasks and their relationships. In Proceedings of the Fourteenth International Conference417 on Artificial Intelligence and Statistics, pages 643–651, 2011.418
[33] Shai Shalev-Shwartz. Online learning: Theory, algorithms, and applications. PhD thesis, The419 Hebrew University of Jerusalem, 2007., 2007.420
[34] Daniel Sheldon. Graphical multi-task learning. Technical report, Computing and Information421 Science Technical Reports, Cornell University, 2008.422
[35] Changjian Shui, Mahdieh Abbasi, Louis-Émile Robitaille, Boyu Wang, and Christian Gagné. A423 principled approach for learning task similarity in multitask learning. In Proceedings of the 28th424 International Joint Conference on Artificial Intelligence (IJCAI-19), pages 3446–3452, 2019.425
[36] Suvrit Sra. Fast projections onto mixed-norm balls with applications. Data Mining and426 Knowledge Discovery, 25(2):358–377, 2012.427
[37] Jonathan Tuck, Shane Barratt, and Stephen Boyd. A distributed method for fitting Laplacian428 regularized stratified models. Journal of Machine Learning Research, 22(60):1–37, 2021.429
[38] Yang Yang, Zhigang Ma, Yi Yang, Feiping Nie, and Heng Tao Shen. Multitask spectral clustering430 by exploring intertask correlation. IEEE transactions on cybernetics, 45(5):1083–1094, 2014.431
[39] Chi Zhang, Peilin Zhao, Shuji Hao, Yeng Chai Soh, Bu Sung Lee, Chunyan Miao, and Steven CH432 Hoi. Distributed multi-task classification: a decentralized online learning approach. Machine433 Learning, 107(4):727–747, 2018.434
[40] Yu Zhang and Qiang Yang. A survey on multi-task learning. IEEE Transactions on Knowledge435 and Data Engineering, 2021.436
[41] Yu Zhang and Dit Yan Yeung. A convex formulation for learning task relationships in multi-task437 learning. In Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence, UAI438 2010, page 733, 2010.439
Checklist440
1. For all authors...441 (a) Do the main claims made in the abstract and introduction accurately reflect the paper’s442 contributions and scope? [Yes]443 (b) Did you describe the limitations of your work? [Yes] For every theoretical claim444 the required assumptions are clearly stated. In general, the acceleration provided by445 MT-OMD is conditioned upon the fact that tasks do present some similarity (i.e., have446 a variance σ2 smaller than 1), what is reflected by the reduced sets of competitors Vσ .447 (c) Did you discuss any potential negative societal impacts of your work? [No] This is a448 theoretically oriented paper, without any foreseeable impact on society.449
(d) Have you read the ethics review guidelines and ensured that your paper conforms to450 them? [Yes]451
2. If you are including theoretical results...452 (a) Did you state the full set of assumptions of all theoretical results? [Yes]453 (b) Did you include complete proofs of all theoretical results? [Yes]454
3. If you ran experiments...455 (a) Did you include the code, data, and instructions needed to reproduce the main experi-456 mental results (either in the supplemental material or as a URL)? [Yes]457 (b) Did you specify all the training details (e.g., data splits, hyperparameters, how they458 were chosen)? [Yes] Some of the details are discussed in the Appendices459 (c) Did you report error bars (e.g., with respect to the random seed after running experi-460 ments multiple times)? [Yes]461 (d) Did you include the total amount of compute and the type of resources used (e.g., type462 of GPUs, internal cluster, or cloud provider)? [N/A]463 4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...464
(a) If your work uses existing assets, did you cite the creators? [N/A]465 (b) Did you mention the license of the assets? [N/A]466 (c) Did you include any new assets either in the supplemental material or as a URL?[N/A]467 (d) Did you discuss whether and how consent was obtained from people whose data you’re468 using/curating? [N/A]469 (e) Did you discuss whether the data you are using/curating contains personally identifiable470 information or offensive content? [N/A]471 5. If you used crowdsourcing or conducted research with human subjects...472
(a) Did you include the full text of instructions given to participants and screenshots, if473 applicable? [N/A]474
(b) Did you describe any potential participant risks, with links to Institutional Review475 Board (IRB) approvals, if applicable? [N/A]476
(c) Did you include the estimated hourly wage paid to participants and the total amount477 spent on participant compensation?[N/A]478

title: Multitask Online Mirror Descent


summaries: The paper considers methods for approximating underdamped Langevin diffusion processes with sum-decomposable strongly convex potential. The main contribution in the paper is the proposal of the Accelerated ULD-MCMC algorithm and its variance-reduced variants. This is followed by a detailed analysis of convergence for each of the algorithms, and the derivation of an information-based lower bound for gradient complexity in the task of estimating ULD that matches the upper bound (in dimensionality, component number, and target accuracy). The authors conclude with a discussion of the optimality (or lack thereof) of ALUM with respect to all the variables. Comprehensive experiments are presented to compare the performance of various algorithms on estimation ULD processes, showing that ALUM and its variants achieve better performance than previous algorithms.

summaries: In the context of sum-decomposable, smooth, and strongly convex potential functions the authors propose algorithms to improve the upper bound on the approximation error for estimating a Underdamped Langevin Diffusion (ULD) process and on the sampling error of the strongly-log-concave sampling problem. The continuous time ULD process has the strongly-log-concave distribution as its invariant distribution so approximating this process and sampling from this distribution are closely related. In the full gradient setting, the authors present a novel algorithm, ALUM which has optimal iteration complexity in $d$, $N$ and $\epsilon$ for the ULD approximation problem. Here, $N$ is the number of terms in the sum-decomposable function, $d$ is the input dimension and $\epsilon$ is the accuracy parameter. In comparison to the RMM algorithm (Lee et al), the dependence on $d$, $N$ and $\epsilon$ in the iteration complexity are matched, but the number of gradient evaluations are reduced by a factor of 2. Next the stochastic gradient oracle setting is considered: here the authors propose algorithms (SVRG-ALUM and SAGA-ALUM) and analyze the corresponding sampling/discretization error. Finally an improvement in the lower bound for the gradient complexity of the ULD approximation problem to $\Omega(N + d^{1/3}N^{2/3}/\epsilon^{2/3})$ is presented. This shows that the gradient complexity of the presented Variance reduction modifications of ALUM match the lower bound in $d$, $N$ and $n$ in the sum-decomposable setting with a stochastic gradient oracle.

summaries: The paper makes new contributions to the literature on sampling from a strongly log-concave smooth measure. This is done by approximating the Underdamped Langevin Diffusion (ULD), and so, a fortiori, the paper also contributes to this approximation problem.  The main result is a new algorithm to approximate the ULD, which is well suited to a stochastic setting, where the gradient of the potential is replaced by certain estimators.  In the full gradient (i.e. without using an estimator for the gradient) setting the new algorithm matches the asymptotic performance of previous works while cutting the number of Gradient queries by about half. In the stochastic setting, if the potential can be represented as a sum of $N$ smooth functions, then the authors also show how to modify their algorithm to accommodate estimators to the gradient and reduce the overall query and iteration complexity.  The results of the stochastic setting are also complemented by a lower bound which essentially show that the dependence on the parameters $N$, and $d$, the dimension are essentially optimal.

summaries: The paper studies the underdamped Langevin diffusion for strongly-convex potential consists of finite sum of smooth components. The authors propose ALUM, which achieves optimal asymptotic complexity in the full gradient setting, and only requires one gradient at each iteration compared to RMM. For the stochastic gradient setting, VR-ALUM methods are proposed which improve over the previous methods in gradient complexity. The authors also accompany with a lower bound showing that the proposed methods are optimal in related parameters for the approximation problem.

summaries: This paper proposes a randomized algorithm for approximating a trajectory of underdamped Langevin dynamics, similar to the RMM method cited in the paper. It further uses variance reduction method to improve the error over sum-decomposable problems. An information theoretic lower bound is established, which the proposed method matches.

abstractText: In the paper, we study the underdamped Langevin diffusion (ULD) with stronglyconvex potential consisting of finite summation of N smooth components, and propose an efficient discretization method, which requires O(N + d 1 3N 2 3 /ε 2 3 ) gradient evaluations to achieve ε-error (in √ E∥·∥2 distance) for approximating d-dimensional ULD. Moreover, we prove a lower bound of gradient complexity as Ω(N + d 1 3N 2 3 /ε 2 3 ), which indicates that our method is optimal in dependence of N , ε, and d. In particular, we apply our method to sample the strongly-log-concave distribution and obtain gradient complexity better than all existing gradient based sampling algorithms. Experimental results on both synthetic and real-world data show that our new method consistently outperforms the existing ULD approaches.

authors:
Andrew Gelman
John B Carlin
Hal S Stern
David B Dunson
Aki Vehtari
Donald B Rubin

title: Bayesian data analysis

venue: CRC press,

authors:
Deepak Ramachandran
Eyal Amir

title: Bayesian inverse reinforcement learning

venue: In IJCAI,

authors:
Peter E Kloeden
Eckhard Platen

title: Numerical solution of stochastic differential equations, volume 23

venue: Springer Science & Business Media,

authors:
Xiang Cheng
Niladri S. Chatterji
Peter L. Bartlett
Michael I. Jordan

title: Underdamped langevin mcmc: A non-asymptotic analysis

venue: In Proceedings of the 31st Conference On Learning Theory,

authors:
Ruoqi Shen
Yin Tat Lee

title: The randomized midpoint method for log-concave sampling

venue: In Advances in Neural Information Processing Systems,

authors:
Benedict Leimkuhler
Charles Matthews
Gabriel Stoltz

title: The computation of averages from equilibrium and nonequilibrium langevin molecular dynamics

venue: IMA Journal of Numerical Analysis,

authors:
GN Milstein
Michael V Tretyakov

title: Quasi-symplectic methods for langevin-type equations

venue: IMA journal of numerical analysis,

authors:
James Foster
Terry Lyons
Harald Oberhauser

title: The shifted ode method for underdamped langevin mcmc

venue: arXiv preprint arXiv:2101.03446,

authors:
Yu Cao
Jianfeng Lu
Lihan Wang

title: Complexity of randomized algorithms for underdamped langevin dynamics

venue: arXiv preprint arXiv:2003.09906,

authors:
Difan Zou
Pan Xu
Quanquan Gu

title: Stochastic variance-reduced hamilton monte carlo methods

venue: In Proceedings of the 35th International Conference on Machine Learning,

authors:
Zhize Li
Tianyi Zhang
Shuyu Cheng
Jun Zhu
Jian Li

title: Stochastic gradient hamiltonian monte carlo with variance reduction for bayesian inference

venue: Machine Learning,

authors:
Rie Johnson
Tong Zhang

title: Accelerating stochastic gradient descent using predictive variance reduction

venue: In Advances in neural information processing systems,

authors:
Aaron Defazio
Francis Bach
Simon Lacoste-Julien

title: Saga: A fast incremental gradient method with support for non-strongly convex composite objectives

authors:
Alain Durmus
Eric Moulines

title: High-dimensional bayesian inference via the unadjusted langevin

venue: algorithm. Bernoulli,

title: Further and stronger analogy between sampling and optimization: Langevin monte carlo and gradient descent

venue: In Conference on Learning Theory,

authors:
Arnak S Dalalyan
Lionel Riou-Durand

title: On sampling from a log-concave density using kinetic langevin diffusions

authors:
Niladri Chatterji
Nicolas Flammarion
Yian Ma
Peter Bartlett
Michael Jordan

title: On the theory of variance reduction for stochastic gradient monte carlo

authors:
Claude JP Bélisle
H Edwin Romeijn
Robert L Smith

title: Hit-and-run algorithms for generating multivariate distributions

venue: Mathematics of Operations Research,

authors:
László Lovász
Miklós Simonovits

title: Random walks in a convex body and an improved volume algorithm

venue: Random structures & algorithms,

authors:
Gareth O Roberts
Richard L Tweedie

title: Geometric convergence and central limit theorems for multidimensional hastings and metropolis algorithms

authors:
Kerrie L Mengersen
Richard L Tweedie

title: Rates of convergence of the hastings and metropolis algorithms

venue: Annals of Statistics,

authors:
László Lovász
Santosh Vempala

title: The geometry of logconcave functions and sampling algorithms

venue: Random Structures & Algorithms,

title: Exponential convergence of langevin distributions and their discrete approximations

authors:
Gareth O Roberts
Jeffrey S Rosenthal

title: Optimal scaling of discrete approximations to langevin diffusions

venue: Journal of the Royal Statistical Society: Series B (Statistical Methodology),

authors:
Raaz Dwivedi
Yuansi Chen
Martin J Wainwright
Bin Yu

title: Log-concave sampling: Metropolis-hastings algorithms are fast

title: Correlation functions and computer simulations

venue: Nuclear Physics B,

authors:
Wenlong Mou
Yi-An Ma
Martin J Wainwright
Peter L Bartlett
Michael I Jordan

title: High-order langevin diffusion yields an accelerated mcmc algorithm

venue: Journal of Machine Learning Research,

authors:
Yi-An Ma
Tianqi Chen
Emily B. Fox

title: A complete recipe for stochastic gradient mcmc

venue: In Proceedings of the 28th International Conference on Neural Information Processing Systems - Volume

authors:
Ilya Sutskever
James Martens
George Dahl
Geoffrey Hinton

title: On the importance of initialization and momentum in deep learning

venue: In International conference on machine learning,

authors:
Chih-Chung Chang
Chih-Jen Lin

title: Libsvm: a library for support vector machines

venue: ACM transactions on intelligent systems and technology (TIST),

authors:
Simon Duane
Anthony D Kennedy
Brian J Pendleton
Duncan Roweth

title: Hybrid monte carlo

venue: Physics letters B,

authors:
Yi-An Ma
Niladri Chatterji
Xiang Cheng
Nicolas Flammarion
Peter Bartlett
Michael I Jordan

title: Is there an analog of nesterov acceleration for mcmc

authors:
Arnak S Dalalyan
Avetik Karagulyan

title: User-friendly guarantees for the langevin monte carlo with inaccurate gradient

venue: Stochastic Processes and their Applications,

text:
1 3N 2 3 /ε 2 3 ) gradient evaluations to achieve ε-error (in √
E∥·∥22 distance) for approximating d-dimensional ULD. Moreover, we prove a lower bound of gradient complexity as Ω(N + d 1 3N 2 3 /ε 2 3 ), which indicates that our method is optimal in dependence of N , ε, and d. In particular, we apply our method to sample the strongly-log-concave distribution and obtain gradient complexity better than all existing gradient based sampling algorithms. Experimental results on both synthetic and real-world data show that our new method consistently outperforms the existing ULD approaches.

heading: 1 Introduction

text:
Sampling is an important research problem in statistics learning with many applications such as Bayesian inference [1], multi-arm bandit optimization [2], and reinforcement learning [3]. One of the fundamental problems in these applications is to sample from a high-dimensional stronglylog-concave distribution. Recently, several Markov chain Monte Carlo (MCMC) based methods were proposed to solve this problem based on underdamped Langevin diffusion (ULD). This continuous diffusion process converges to the target distribution exponentially fast. Thus, the methods approximating a ULD process could be used to sample from the target distribution within certain accuracy.
Multiple discretization methods have been proposed for approximating ULD. Among them, the Euler-Maruyama discretization [4] is the simplest one but generates the largest error. Recently the left point method (LPM) 1 [5] was introduced to fix the gradient term in ULD to be the gradient at k-th iteration, and then integrate the new linear stochastic differential equation (SDE) with a small time-interval. Subsequently, Shen and Lee [6] proposed randomized midpoint method (RMM) with smaller error. There are also discretization schemes based on splitting [7] or Runge-Kutta method [8, 9]. More recently, Cao et al. [10] derived an information-based complexity lower bound for simulating a d-dimensional ULD. Under the assumption that the full gradient oracle ∇f(x) is evaluated at most n times, they show a lower bound for worst-case error by perturbation analysis, which matches the discretization error upper bound of RMM in the dependence of d and n.
Although the ULD-MCMC methods with full gradient oracle are largely understood, many real-world applications involve summation form of potential function and large-scale data, which leads to the need of stochastic gradient methods. The vanilla stochastic gradient methods have been used to replace full gradient [5]. Albeit the computational cost for each iteration is reduced, the variance of
1Although this method is mostly just denoted as ULD-MCMC, we adopt the name LPM to distinguish it from other discretization methods. The name comes from the fact that gradient is evaluated at the left point of the time interval.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
stochastic gradient is much larger than the discretization error and therefore degenerates the overall performance. Previous works [11, 12] used stochastic variance reduced gradient (SVRG) [13] and SAGA [14] instead, but the gradient complexities of these methods are still worse than the full gradient RMM in terms of dependence on accuracy ε. Thus, there exists a natural question:
What is the optimal ULD-MCMC method with sum-decomposable potential?
In this paper, we focus on optimal dependence of dimension d, components number N and accuracy ε in gradient complexity for estimating a ULD process. We answer this question by two parts. We first provide a novel ULD-MCMC method and derive the corresponding complexity upper bound in Sections 4 and 5. After that, we analyze the worse case error and show that the lower bound matches the upper bound in Section 6. The major contributions of our paper can be summarized as follows.
1. We propose a new full gradient ULD-MCMC method, called as AcceLerated ULD method (ALUM), whose discretization error has the same order dependence on dimension d, step size h as RMM. Although RMM already has optimal asymptotic complexity in full gradient setting, ALUM is still of practical interest. Compared with RMM, which uses two gradient evaluations at each iteration, ALUM uses gradient less frequently and only requires one gradient at each iteration to achieve constant speedup.
2. We further propose VR-ALUM methods, including SVRG-ALUM and SAGA-ALUM, which utilize the unbiased variance reduction techniques in ALUM under sum-decomposable setting. We show that these methods achieve better gradient complexity than all existing gradient based MCMC approaches. These gradient complexities for sampling from a strongly-log-concave distribution are compared in Table 1.
3. We derive an information-based lower bound on worst-case error for estimating a ULD process with only gradient oracle and weighted Brownian motion oracle. We show that in order to achieve ε approximation accuracy, Ω(N + d 1 3N 2 3 ε− 2 3 ) single component gradient evaluations
are needed. This lower bound matches the upper bound for VR-ALUM in terms of dependence of dimension d, sample size N , and accuracy ε. Therefore, our VR-ALUM methods are indeed optimal for estimating a ULD process under sum-decomposable setting.
2In Zou et al. [11], the authors call their method SVR-HMC. However, their method is not based on Hamiltonian Monte Carlo (HMC), but based on ULD. Their method is just applying SVRG to replace full gradient in LPM.
3Further explanation and comparison is shown in Appendix A.2.

heading: 2 Related work

text:
Though we mainly use the gradient oracle in this paper, there is also a large body of sampling algorithms which leverage zeroth order potential value oracle. For example, Metropolis-Hastings accept-reject step could be used to ensure an MCMC converges to a stationary distribution equal to the target distribution [19, 20, 21, 22, 23] and the linear convergence can be obtained. For example, MALA [24, 21, 25], as Metropolis-adjusted ULA, can achieve ε error in total variation distance within O(d2 log(1/ε)) steps under certain initialization [26]. Note that purely gradient based sampling algorithms typically only converge sub-linearly because the stationary distribution is different from the target distribution.
Other diffusion processes could also be used for constructing MCMC. ULA, a discretization of Langevin diffusion (LD) [27, 24], is the first commonly used gradient-based MCMC. We will provide an overview for some LD-based and ULD-based methods in Appendix A.1. Later, higher order diffusion process is also discussed. Although the path is even smoother than ULD, currently it is not clear whether such higher order smoothness of path could be leveraged to accelerate the convergence without extra assumptions on potential function. Mou et al. [28] studied third-order diffusion, but their acceleration requires special structure or high order smoothness of potential. The more general diffusion process, which converges to the target distribution, was studied in Ma et al. [29].

heading: 3 Preliminary

text:
The problem of sampling from a strongly-log-concave distribution involves a probability density function p∗(x) defined on a real vector space Rd. A corresponding potential function f(x) = − log(p∗(x)) can be defined such that f(x) is strongly convex. One way to solve this sampling problem is constructing a Markov chain that converges to a stationary distribution that is the same as or similar to the target distribution. In Appendix A.1, we introduce several such Markov chains as discretization of certain continuous stochastic processes, and we roughly analyze discretization error of these methods with pointing out the bottleneck.
We use Õ(f) = O(f) logO(1)(f) to omit logarithm factor. ∥·∥2 means the Euclidean norm. We define a norm of random vector as ∥·∥L2 = √ E∥·∥22. Next, we list our assumptions on potential f(x).
Assumption 1 (Sum-decomposable). f(x) = ∑N
i=1 fi(x), where integer N is the sample size.
Assumption 2 (Smoothness). Each function fi is twice differentiable on Rd and there exists a constant L > 0, such that ∇2fi(x) ≼ LN I for any x ∈ R
d where I is the identity matrix. It can be easily verified that f(x) is L-smooth.
Assumption 3 (Strong Convexity). There exists a constant m > 0 such that: f(x) − f(y) ≥ ⟨∇f(y),x− y⟩+ m2 ∥x− y∥ 2 2. We define the condition number κ := L/m.
We finally define the 2-Wasserstein distance between distributions. For any pair of probability measures µ and ν on the same parameter space, a transference plan ζ between µ and ν is a joint distribution such that the marginal distributions on two sets of coordinates are µ and ν, respectively. We denote Γ(µ, ν) as the set of all transference plans, and define the 2-Wasserstein distance between µ and ν as follows:
W 22 (µ, ν) = inf ζ∈Γ(µ,ν)
∫ ∥x− y∥22dζ(x, y).

heading: 4 AcceLerated ULD-MCMC (ALUM) methods

text:
In this section, we propose a class of AcceLerated ULD-MCMC (ALUM) methods based on approximation of the continuous ULD process:
dXt = Vtdt, dVt = −∇f(Xt)dt− γVtdt+ √ 2γdBt. (1)
The gradient term in the above SDE could be highly non-linear such that closed form solution is not available. Thus, we propose an estimation of SDE solution at time point h which only requires
gradient evaluation at one single point:
X (o) h = X0 + ψ1(h)V0 − hψ1(h− ah)∇f(X (e) ah ) + ex,[0,h],
V (o) h = ψ0(h)V0 − hψ0(h− ah)∇f(X (e) ah ) + ev,[0,h],
X (e) ah = X0 + ψ1(ah)V0 + ex,[0,ah],
(2)
where the random variable a is uniform random variable in [0, 1], ψi(·) and noise terms ex,[0,h], ev,[0,h], ex,[0,ah] are defined as follows 4.
ψ0(h) = e −γh, ψ1(h) =
1 γ (1− e−γh),
ex,[0,h] = √ 2γ ∫ h 0 ψ1(h− s)dBs, ev,[0,h] = √ 2γ ∫ h 0 ψ0(h− s)dBs. (3)

heading: 4.1 Full gradient ALUM

text:
In this subsection, we propose a full gradient ALUM method, which is shown in Algorithm 1. Specifically, the update rule is obtained by setting x(o)k+1,v (o) k+1 to X (o) h ,V (o) h where the ULD starts from X0 = x (o) k ,V0 = v (o) k and h is a small step size.
Algorithm 1: Full gradient ALUM Method
Input: Initial point (x(o)0 ,v (o) 0 ), parameter γ, iteration number K, and step size h > 0. for k = 0 to K − 1 do
Randomly sample ak uniformly from [0, 1]; Generate ex,[0,h],k, ev,[0,h],k, ex,[0,akh],k according to Appendix A.5; x (e) k = x (o) k + ψ1(akh)v (o) k + ex,[0,akh],k; Calculate full gradient ∇f(x(e)k ); x (o) k+1 = x (o) k +ψ1(h)v (o) k −hψ1(h−akh)∇f(x (e) k )+ex,[0,h],k; v (o) k+1 = ψ0(h)v (o) k −hψ0(h−akh)∇f(x (e) k )+ev,[0,h],k;
end for Output: x(o)K . Next, we compare our method with both RMM and Nesterov’s accelerated gradient method (NAG) which has been shown optimal for gradient based optimization. Compared to RMM: By comparing the formula in (19) and (2), it is easy to see that ALUM and RMM differ in one single term −ψ2(ah)∇f(X0). Our motivation of dropping this term is to reduce the computations. By deleting this term from the algorithm, only one gradient evaluation at a randomized midpoint is needed at each iteration. This compares favorably to the RMM, which requires two gradient evaluations at each iteration.
However, could this method still estimate ULD accurately with only half gradient evaluations? We firmly answer this question with rigorous error analysis in Theorem 2. Roughly speaking, if we only consider the error’s dependence on step size h, the answer is yes. In this section, we only give an intuition of how that is possible.
In Appendix A.1, we show that RMM has very low bias and high variance. Therefore, the bottleneck is the variance, and increasing the bias slightly would not degenerate the overall performance much if the bias is not larger than the error introduced by variance. The dropped term −ψ2(ah)∇f(X0) has norm O(h2), and the coefficient of gradient in Vh is O(h). Thus, the bias introduced in single step is O(h3) when we use step size h. After accumulating for T/h = O(h−1) iterations, the bias is O(h2), which still has better dependence on h than the square root of variance O(h3/2).
With increased bias, the complexities of ALUM and some other variants5 still have the same dependence on ε as RMM, the highest order dependence on κ could deteriorate as discussed in Section 7.
Compared to Nesterov’s accelerated gradient (NAG) method: NAG is a gradient based optimization method [30]. Based on [31], NAG could be formulated as the following momentum method:
4ex,[0,ah] is obtained by simply substituting h with ah in the definition of ex,[0,h]. 5We can drop ex,[0,h],k and ex,[0,akh],k in (2) to derive other variants of ALUM. The bias of these variants increases to O(h3/2), which is still no larger than the square root of variance. The maximum order dependence on h is the same, thus the final iteration complexity has the same dependence on ε as RMM.
xk+1 = xk + ckxk + c ′ k∇f(x (e) k ), vk+1 = ckxk + c ′ k∇f(x (e) k ), x (e) k = xk + ckxk . (4)
We can see that the main differences between Algorithm 1 and (4) are coefficients and additional Gaussian noise terms. The coefficients in NAG are set in a deterministic way. However, the coefficients in ALUM could be random values from fixed distributions. The additional noise terms in ALUM come from the Brownian motion term in ULD and are necessary for sampling methods.
Despite the differences, both NAG and ALUM use one gradient at each iteration, and both of them take a big jump along current momentum direction to calculate the gradient instead of directly computing the gradient at current iterate.

heading: 4.2 Variance-reduced ALUM (VR-ALUM)

text:
In the subsection, we propose variance-reduced stochastic ALUM (VR-ALUM) methods based on two common unbiased 6 variance-reduced techniques: stochastic variance reduced gradient (SVRG) [13] and SAGA [14]. VR-ALUM comes from simply replacing all full gradient in Algorithm 1 with a gradient estimation ∇̃k. We show Algorithm 2 in the Appendix A.6 due to the limit of space. Next, we briefly introduce these two variance reduction techniques. SVRG utilizes the following gradient estimation (Bk is the batch of k-th iteration and b is the batch size): ∇̃SVRGk = N b ∑ i∈Bk ( ∇fi(x(e)∇̃k ) − ∇fi(x) ) + ∑N i=1 ∇fi(x). The full gradient ∑N i=1 ∇fi(x) and point x are updated after every τ evaluations of ∇̃SVRGk . We call the hyperparameter τ as epoch length.
SAGA estimates the gradient in the following way, where ϕik+1 is set as x (e)∇̃ k if and only if i ∈ Bk,
otherwise ϕik+1 = ϕ i k: ∇̃SAGAk = N b ∑ i∈Bk ( ∇fi(x(e)∇̃k ) − ∇fi(ϕik) ) + ∑N i=1 ∇fi(ϕik). SAGA does not re-compute but stores the latest gradient information ∇fi(ϕik) for each fi. Therefore, SAGA does not introduce extra gradient evaluation except for initialization, but has much larger storage requirements.
Both SVRG and SAGA are unbiased, which means EBk∇̃k = ∇f(x (e)∇̃ k ), where EBk means expectation over random batch at k-th iteration. Moreover, both SVRG and SAGA reduce the meansquared error of gradient estimation and satisfy bounded MSE property proposed in Appendix B.3.

heading: 5 Theoretical analysis

text: We provide non-asymptotic upper bounds on sampling error and discretization error for our methods, including full gradient ALUM and VR-ALUMs. The proof is shown in Appendix B. Throughout this section, we assume W2(p0, p∗) = O(1) √ d/m when deriving asymptotic results. This means the initialization is not too far away from the target distribution, and can be achieved under multiple setting as discussed in Appendix A.7.

heading: 5.1 Convergence analysis of full gradient ALUM

text:
Recall that ALUM can be used for solving two different but related problems: strongly-log-concave sampling and approximating the ULD. We show the upper bound for sampling error in Theorem 1 and the upper bound for approximation error in Theorem 2 separately. Theorem 1. Suppose Assumptions 1 to 3 hold. Given an initial distribution p0(x), we initialize ALUM with random x(o)0 based on probability p0 and random v (o) 0 from standard Gaussian distribution. Assume L = 1 and let γ = 2, pk be the distribution of x (o) k and p
∗ be the target distribution. Assume we use step size h ≤ m22 . After running the ALUM for k iterations, we have the following upper bound of sampling error in 2-Wasserstein distance:
W2(pk, p ∗) ≤ 2(1− mh
4 )kW2(p0, p ∗) + 12
√ h3d
m . (5)
6The reason we choose unbiased variance reduction instead of biased one is that bias accumulates quicker than variance, therefore generates higher overall error.
Theorem 2. With same assumptions and setup in Theorem 1, we consider Xkh which comes from a continuous ULD starting from the same position as x(o)0 with probability 1. We have the following upper bound of discretization error:
∥x(o)k −Xkh∥L2 ≤ 38 √ h3
m W2(p0, p
∗) + 31
√ h3d
m . (6)
Remark 1. The assumption L = 1 actually does not limit the availability of the algorithm. For any function f(x) with L ̸= 1, we can define another function f ′(x′) = f( 1√
L x′) which satisfies L = 1.
ALUM with step size h ≤ m22L and initial momentum from the standard Gaussian distribution can be used to sample x′ from the distribution p′(x′) ∝ exp(−f ′(x′)). The sample x from distribution p(x) ∝ exp(−f(x)) could be obtained by a transform x = 1√
L x′. This is essentially the same as
directly incorporating L into ULD process as in Appendix A.3.
Thus, we can derive the iteration complexity for both problems. For sampling problem, we define the iteration complexity as the number of iterations K needed to achieve W2(pK , p∗) ≤ ε √ d/m with certain step size h. For approximating the ULD at a given time point T , we define the iteration complexity as the number of iterations K needed to achieve ∥x(o)K − XT ∥L2 ≤ ε with step size h = T/K. The results are shown in Table 2. Detailed derivations can be found at Appendix A.8.

heading: 5.2 Convergence analysis of variance-reduced ALUM (VR-ALUM)

text:
We show the upper bound for sampling error in Theorem 3 and the upper bound for discretization error in Theorem 4 separately. Theorem 3. With same assumptions in Theorem 1, we use SVRG-ALUM with epoch length τ = ⌈N/b⌉ or SAGA-ALUM. We introduce an extra assumption h3 ≤ 12304cb
3mN−2. We have the following upper bound of sampling error in 2-Wasserstein distance.
W2(pk, p ∗) ≤ 2(1−mh
4 )kW2(p0, p ∗)+92
√ h3
m
√ c N
b 3 2
W2(p0, p ∗)+(12+57 √ c N
b 3 2
))
√ h3d
m . (7)
The constant c is defined as c = 1 for SVRG-ALUM and c = 2 for SAGA-ALUM. Theorem 4. With the same assumptions in Theorem 3, we consider Xkh which comes from a continuous ULD starting from the same position as x(o)∇̃0 with probability 1. We have the following upper bound of discretization error.
∥x(o)∇̃k −Xkh∥L2 ≤ (38 + 92 √ c N
b 3 2
)
√ h3
m W2(p0, p
∗) + (31 + 57 √ c N
b 3 2
)
√ h3d
m . (8)
Remark 2. Theorems 1 to 4 are specializations of more general results in Appendix B.1, where full gradient, SVRG and SAGA are unified under the framework of bounded MSE property that is defined in Appendix B.3. A unified approach not only simplifies the proof, but also indicates that our analysis could easily generalize to other gradient estimations that satisfy bounded MSE property.
Similar to the full gradient case, we derive the iteration complexity K in Appendix A.9. Moreover, we define gradient complexity as the number of single component gradient evaluation ∇fi(x) needed to achieve certain accuracy. We show the results in Table 3 and add the derivations in Appendix A.9. We finally simplify the result by only considering the dependence of d, N , b, and ε.
Corollary 1. When b ≤ O(N 23 ), the gradient complexity of SAGA-ALUM and SVRG-ALUM for sampling problem is Õ(N +N 2 3 ε− 2 3 ) and their gradient complexity for ULD approximation problem is O(N + d 1 3N 2 3 ε− 2 3 ).

heading: 6 Information-based complexity

text: We first declare our setup for the problem class and accessible information. Then we show informationbased lower bound for approximation error and oracle complexity. The proof is shown in Appendix C.

heading: 6.1 Setup

text:
Class of possible potential functions: We denote the class of possible potential functions as U = U(d,N,m,L) = {(f1, . . . , fN )|∇2fi(x) ≼ (L/N)I,mI ≼ ∇2 ∑N i=1 fi(x) and ∥Ep∗ [X]∥2 ≤√
d/m} where 0 < m < L and p∗(X) ∝ exp(− ∑N
i=1 fi(X)).
Besides Assumptions 1 to 3 used in the previous sections, we introduce a new assumption that the target distributions have mean in a ball of constant radius around the origin. This is necessary because the lower bound could be arbitrarily large if the mean is far away from the initialization point.
Solution mapping: We denote the probability space for Brownian motion in the ULD process as (M,Σ,P). Then, the true solution of ULD process starting from the origin 0 at time T can be denoted as the solution mapping XT : (ω,U) ∈ (M× U) 7→ XT (ω,U) ∈ Rd. Gradient oracle: For a given set of potential functions U = (f1, . . . , fN ), the single component gradient oracle is ΥU : (i,x) ∈ [N ]× Rd 7→ ∇fi(x) ∈ Rd. Brownian oracle: We assume the Brownian motion at a given time t > 0 could be evaluated with oracle Bt(ω) ∈ Rd for any event ω ∈ M. We further assume that the weighted Brownian motion is also admissible. B(θ)t (ω) = ∫ T 0 eθsdBs(ω).
Deterministic algorithm: A deterministic algorithm starts from empty information I0 = (). At i-th step, one oracle and corresponding parameters are picked by certain procedure. If the gradient oracle is picked, the algorithm will generate an index i and a point x. If the weighted gradient oracle is picked, the parameters are order θ for weighted Brownian motion and time t. The picked oracle and parameters are represented by ϕi(Ii). The evaluation result is represented by Υ(ϕ(Ii), ω, U) where ω ∈ M and U ∈ U . The picked oracle, parameters, and result will be stored as new information, hence Ii+1 = (Ii, ϕ(Ii),Υ(ϕ(Ii), ω, U)).
We consider deterministic algorithms that stop at n-th step. The final estimation is generated with a mapping Y (In) ∈ Rd. We expect Y (In) to be as close as possible to true solution XT (ω,U).
We denote a deterministic algorithm as a mapping A from M× U to Rd and A(ω,U) = Y (In) for some ϕ and Y . The family of all such algorithms is denoted by ADetn .
Randomized algorithm: We consider another probability space (M̃, Σ̃, P̃) as the source of randomness. A randomized algorithm A with n steps is a mapping from M× M̃× U to Rd such that A(·, ω̃, ·) ∈ ADetn for any ω̃ ∈ M̃. The family of all such randomized algorithms is denoted by An. Worst error of algorithms: We care about the following worst-case error for any possible algorithms:
e2A,U := inf A∈A sup U∈U Eω∈PEω̃∈P̃∥XT (ω,U)−A(ω, ω̃, U)∥ 2 2 (9)
We always assume T > 0 to avoid trivial cases.

heading: 6.2 Lower bounds

text:
We first provide lower bounds on worst-case estimation error.
Theorem 5. When n < N which means that gradient evaluation number is less than components number, we have e2An,U ≥ dC1, where C1 is positive and independent of d, N , and n.
The above theorem is based on the fact that no algorithm can accurately estimate the global minimum point of the sum of quadratic potentials ∑N i=1 fi(x) with information from only N − 1 components.
Theorem 6. When the gradient evaluation number n is a multiple of N , we have e2An,U ≥ dC2 N2
n3 , where C2 is positive and independent of d, N , and n.
The method we use to prove the above theorem is to show that an algorithm that uses a limited number of gradient evaluations cannot distinguish a class of perturbed quadratic functions.
Next, we use the above lower bounds on error to derive a lower bound on gradient complexity.
Corollary 2. For small enough target accuracy ε such that ε2 < dC1, in order to achieve eAn,U ≤ ε, the minimum number of single component gradient oracle evaluations is Ω(N + d 1 3N 2 3 ε− 2 3 ).
This lower bound matches the upper bound in Corollary 1 which indicates that gradient complexity of variance-reduced ALUM for estimating a ULD is optimal in the dependence of d, components number N , and approximation accuracy ε.

heading: 7 Optimality

text:
In this section, we discuss in what sense our ALUM is optimal (or not), and point out possible improvements for future work.
Optimal for approximating problem: ALUM can be used for two tasks: (1) estimating a ULD process, (2) sampling from a strongly-log-concave distribution. ALUM is only optimal for the first task, and is not necessarily optimal for the second one. For example, we see in Section 5 that there exists a logarithm factor in the gradient complexity for sampling problem. We believe some proper adaptive step size method could cancel that extra term.
Dependency on κ: ALUM is only optimal on the dependence of d, N , ε, but not optimal in κ dependence. Actually both our method and our analysis may not be optimized for the dependence of κ. Currently, there is no information-based lower bound with clear dependence on condition number κ, therefore no matter how good the dependence on κ is, it is not sufficient to say it is optimal. For sampling problem, it is believed that O(min(κ, d2)) is the “natural barrier” [6] for iterations needed for achieving W2 ≤ ε √ d/m with ε = 1/2. Currently, the best dependence on κ is achieved by RMM with O(κ7/6). For ALUM, the maximum dependence is O(κ2). This dependence is worse for two reasons. First, the analysis is not optimized for κ. Based on the analysis of bias and variance in Appendix A.10, we conjecture that the dependence could be improved to O(κ3/2). Second, the method is not optimized for κ. We save one extra gradient compared to RMM, at the price of a slightly increased bias. This degenerates the dependence of κ.
Despite the max dependence order being larger than 1, in the high-precision regime, when ε is small enough, the term κ/ε 2 3 is dominant, therefore, both full-gradient ALUM and full-gradient RMM achieve O(κ) dependence.
We note the complexity for VR-ALUM increases the dependence of κ to κ4/3 in high-precision regime. We conjecture that this is an artifact of error-based analysis compared to momentum-based analysis, and could be improved by some tighter analysis.
Assumptions: Finally, we point out that our method and analysis are based on Assumptions 1 to 3 and the gradient oracle. It is natural to obtain better algorithm by introducing new assumptions or new oracle.
Many higher order integrators [8, 9] leverage higher order smoothness assumption to reduce discretization error. Assumptions on structure of potential [28] have also been shown to make acceleration possible together with high order diffusion process. Apart from the gradient oracle, many Metropolisadjusted algorithms [24, 26] leverage zeroth-order oracle to achieve linear convergence. Second order oracle [17] has also been incorporated into estimating ULD.

heading: 8 Experiments

text:
In this section, we compare different algorithms for estimating a ULD process on both Gaussian model and Bayesian logistic regression model. For the Gaussian model, the potential is defined as:
fi(x) = 1
2N (di − x)⊤Σ−1(di − x) , (10)
where di and Σ is generated randomly to satisfy d = 5,N = 100,m = 1 and L = 10. For logistic regression model, the potential is:
fi(x) = m
2N ∥x∥22 + N∑ i=1 log(1 + exp(−yia⊤i x)) , (11)
where m is set such that κ = 104 and yi, ai are data points in australian dataset from LIBSVM [32]. We calculate the trajectory error which is defined as 1K ∑K i=1 √ ∥xk − x′k∥22 + ∥vk − v′k∥22, where xk and vk are generated by ALUM, VR-ALUMs or other algorithms. x′k and v ′ k come from a reference path that is very close to true solution. We specify how we generate this reference path in Appendix D.1. The reason for us to average the error along the path instead of just reporting error at the final iterate is that we need more data points to reduce the variance.
Figure 1 shows the error for ALUM, LPM, and RMM with full gradient, stochastic gradient, SVRG and SAGA.8 The detailed setup can be found in Appendix D.2. We summarize the messages in Figure 1 as follows:
• SAGA-ALUM achieves the best efficiency in the sense that with same number evaluations of single component gradient ∇fi(x), SAGA-ALUM has smaller discretization error than any other algorithms. Results on more dataset are shown in Appendix D.3.
• Variance reduced algorithms constantly outperform full gradient algorithm of the same type by a large margin.
• For full gradient method, the discretization error of ALUM and RMM has similar asymptotic dependence on gradient evaluation number, and they are better than LPM. This phenomenon will be shown clearer in Figure 2.
• ALUM achieves constant acceleration compared to RMM by saving one gradient evaluation per iteration.
7Most of the error bars are too small to be visible. The SG-ALUM and SG-LPM highly overlap. 8Currently there is no theoretical result for RMM with stochastic gradient, SVRG or SAGA. We also didn’t provide a theoretical result for ALUM with stochastic gradient. However, this doesn’t prevent us from evaluating them experimentally.
We also show the relationship between discretization error and step size in Figure 2. Only full gradient method are shown here due to the limit of space, and more results for VR-ALUMs are shown in Appendix D.4.
Our analysis in Theorem 2 gives upper bound of discretization error of ALUM as O(h3/2). Figure 2 shows that our analysis is tight. Moreover, the discretization error of ALUM is almost the same as RMM. Due to the fact that ALUM requires only half gradient evaluations per step than RMM, ALUM achieves better efficiency.
Next, we discuss the effect of batch size. According to our theory in Corollary 1, our method is not sensitive to this hyperparameter introduced by SAGA or SVRG when batch size is relatively small, as the gradient complexity remains the same for b = O(N2/3). We verify that in Figure 3, where the sampling efficiency is almost same for small batch sizes, and only deteriorate for very large batch size. We give further discussion on why our method is not sensitive to step size in Appendix D.5.
Finally, we apply ALUM and VR-ALUMs to sample from a target distribution in Appendix D.6.

heading: 9 Conclusion

text:
In this paper, we propose a class of MCMC methods for finite sum form of strongly-convex potential. Our methods are proven to be optimal in the sense that the discretization error has the best possible asymptotic dependence on dimension, number of potential summands, and number of gradient evaluations. Experiments on both synthetic and real data verify the superior performance of our algorithm. We also discuss possible improvements for future work.
Acknowledgments and Disclosure of Funding
This work was partially supported by NSF IIS 1845666, 1852606, 1838627, 1837956, 1956002, OIA 2040588.

title: Optimal Underdamped Langevin MCMC Method