Adapting to Continuous Covariate Shift via
Online Density Ratio Estimation
Yu-Jie Zhang1 , Zhen-Yu Zhang2 , Peng Zhao3 , Masashi Sugiyama2,1
1

The University of Tokyo, Chiba, Japan
2
RIKEN AIP, Tokyo, Japan
3
National Key Laboratory for Novel Software Technology, Nanjing University, Nanjing, China

Abstract
Dealing with distribution shifts is one of the central challenges for modern machinelearning. One fundamental situation is the covariate shift, where the input distributions of data change from the training to testing stages while the input-conditional
output distribution remains unchanged. In this paper, we initiate the study of amore challenging scenario — continuous covariate shift — in which the test data
appear sequentially, and their distributions can shift continuously. Our goal is toadaptively train the predictor such that its prediction risk accumulated over time canbe minimized. Starting with the importance-weighted learning, we theoreticallyshow the method works effectively if the time-varying density ratios of test andtrain inputs can be accurately estimated. However, existing density ratio estimation
methods would fail due to data scarcity at each time step. To this end, we propose
an online density ratio estimation method that can appropriately reuse historicalinformation. Our method is proven to perform well by enjoying a dynamic regretbound, which finally leads to an excess risk guarantee for the predictor. Empiricalresults also validate the effectiveness.

1

Introduction

How to deal with distribution shifts is one of the central challenges for modern machine learning [1, 2],which is also the key requirement of robust artificial intelligence in open and dynamic environments [3,4]. Covariate shift is a representative problem setup and has attracted much attention [5, 6], where theinput density D(x) varies from training data to test data but their input-conditional output D(y | x)
remains unchanged. On the one hand, covariate shift serves as a foundational setup for theoreticallyunderstanding general distribution shifts [7], and on the other hand, it can encompass many real-world
tasks such as brain-computer interface [8], speaker identification [9] and audio processing [10].
Existing works of coping with covariate shift mainly focused on the “one-step” adaptation, where
the learner aims to train a model well-performed on a fixed testing distribution. To mitigate the
distribution discrepancy, a common and classic solution is the importance-weighting framework [6],where one assigns an appropriate importance weight to each labeled training sample and then conducts
weighted empirical risk minimization. The importance weight, also known as the density ratio of test
and training inputs, usually needs to be estimated via a reasonable amount of unlabeled data sampledfrom the testing distribution [11, 12]. However, the one-step adaptation can be insufficient in many
real-world covariate shift scenarios, especially when data are accumulated in an online fashion suchthat testing environments continuously shift and only a few unlabeled samples are observed at eachtime. For instance, in the speaker identification task, the speech features vary over time due to session∗

Correspondence: Yu-Jie Zhang <yujie.zhang@ms.k.u-tokyo.ac.jp> and Peng Zhao <zhaop@lamda.nju.edu.cn>

37th Conference on Neural Information Processing Systems (NeurIPS 2023).

dependent variation, the recording environment change, and physical conditions and emotions of the
speakers [9]. Therefore, it is necessary to perform a prompt adaptation to the changes.
Motivated by such demands, we initiate the study of continuous covariate shift, where the learner’s
goal is to adaptively update the model to minimize the risk accumulated over time (see the formal
definition in Section 3). To achieve this goal, we propose the ACCOUS approach (Adapt to Continuous
Covariate shift with Unlabeled Stream) equipped with sound theoretical guarantees. Our approach isbased on the classic importance-weighting framework yet requires innovations to make it applicableto the continuous shift scenario. Indeed, we theoretically identify that the importance-weighting
works effectively if the cumulative estimation error of the time-varying density ratios of train andtest inputs can be reduced. However, at each time step, a direct application of one-step density ratio
estimation would lead to a high variance due to data scarcity; whereas reusing all previous data can behighly biased when testing distributions change dramatically. Thus, it is crucial to design an accurate
estimation of time-varying test-train density ratios by appropriately reusing historical information.
To combat the difficulty, we propose a generic reduction of the time-varying density ratio estimation
problem to the online convex optimization [13]: one can immediately obtain high-quality time-varying
density ratios by a suitable online process to optimize its dynamic regret over a certain sequence of loss
functions. Our reduction is based on the Bregman divergence density ratio matching framework [14]and applicable to various existing density ratio estimation models with specific configurations of thedivergence functions. To minimize the dynamic regret, our key algorithmic ingredient of this online
optimization is the online ensemble structure [15, 16], where a group of base-learners is maintained
to perform density ratio estimation with different lengths of historical data, and a meta-algorithm
is employed to combine the predictions. As such, we can properly reuse historical information
without knowing the cumulative intensity of the covariate shift. We further instantiate the reductionPT
framework with the logistic regression model [17]. Letting VT = t=2 ∥Dt (x) − Dt−1 (x)∥1 , wee 1/3 V 2/3 ) dynamic regret bound for the density ratio estimator, which finally leads toprove an O(TTe −1/3 V 1/3 ) averaged excess risk of the predictor trained by importance-weighted learning.
an O(TT
The rate can hardly be improved even if one receives labels of the testing stream after prediction (seemore elaborations below Theorem 3 and Appendix D.6). Finally, we conduct experiments to evaluate
our approach, and the empirical results validate the theoretical findings.
Technical Contributions. We note that online ensemble was also employed by [18] to handle
continuous label shift, another typical distribution shift assuming the change happens on the classprior Dt (y). However, their method crucially relies on the unbiasedness of weights in the risk
estimator, whereas density ratio estimators for covariate shift adaptation cannot satisfy the conditionin general. Instead of pursuing unbiasedness, we delicately design an online optimization process
to learn the density ratio estimator, which is versatile to be implemented with various models.
Our methodology is very general and might be of broader use, such as relaxing the unbiasedness
e 1/3 V 2/3 ) dynamic regret bound for therequirement in the continuous label shift. Besides, the O(TTlogistic regression density ratio estimator is obtained non-trivially. Our bound essentially holds
for online learning with exp-concave functions under noisy feedback. The only one achieving
this [19] holds in expectation and requires complicated analysis on the Karush-Kuhn-Tucker (KKT)
condition [20] of comparators. By contrast, our bound holds in high probability, and the analysis is
greatly simplified by virtue of exploiting the structure of comparators in our problem that they are
essentially the minimizers of expected functions (hence without analyzing the KKT condition).

2

Preliminaries

This section introduces the preliminaries and related work on the continuous covariate shift adaptation.
We discuss the related work on non-stationary online learning in Appendix B
One-Step Covariate Shift Adaptation. Let D0 (x, y) and D1 (x, y) be the training and test distributions. The “one-step” adaptation problem studies how to minimize the testing risk R1 (w) =0E(x,y)∼D1 (x,y) [ℓ(w⊤ x, y)] by the model trained with a labeled dataset S0 = {xi , yi }N
i=1 sampled
N1
from D0 (x, y) and unlabeled dataset S1 = {xi }i=1 sampled from D1 (x). We call this one-step
adaptation since the test distribution is fixed and a reasonable number of unlabeled data S1 is available.
Importance-Weighted ERM. A classic solution for the one-step covariate shift adaptation is theimportance-weighted empirical risk minimization (IWERM), which mitigates the distribution shift2e1 (w) = E(x,y)∼S [r∗ (x)ℓ(w⊤ x, y)], where r∗ (x) =by minimizing the weighted empirical risk R
110
e1 (w) is unbiased to R1 (w) and thus the learned
D1 (x)/D0 (x) is the importance weight. The risk R
model is consistent over the test distribution [5]. The importance weighted learning was studied fromthe lens of variance-bias trade off [5, 21], cross validation [22], model misspecification [23], and deepneural network implementation [24]. All those are conducted under the one-step adaptation scenario.
Density Ratio Estimation. The importance weighted estimation, also known as the density ratio
estimation (DRE) [12], aims to estimate the density ratio r1∗ (x) using datasets input of S0 and S1 .
Various methods were proposed with different statistical models [14, 17, 25, 26, 27, 28, 29]. All thosemethods focused on the one-step adaptation, and it is challenging to extend them to the continuous
shift due to the limited unlabeled data at each time step. The work [30] studied how to update thedensity ratio estimator with streaming data, but the ground-truth density ratio is assumed to be fixed.
Recently, the time-varying density ratio was investigated in [31]. Their problem setup fundamentallydiffers from ours as they assume sufficient online data at each iteration, while our challenge is dealing
with limited data per iteration. The work [32] proposed to learn the time-varying density ratios using
all past data at each time. Although learning theory insights are provided in the paper, it is still
unclear how the learned density ratios balance sample complexity with environmental shift intensity.
In contrast, our density ratio estimator updates in an online fashion with dynamic regret guarantees.
Continuous Distribution Shift. For broader continuous distribution shift problems, the study [33]focused on the label shift case and provided the first feasible solution. Then, the work [18] introduced
modern non-stationary online learning techniques to the problem, developing the first method withdynamic regret guarantees. As shown in Remark 2, our solution refines the previous methods [33, 18]
by decoupling density ratio estimation from predictor training, paving the way for a theoreticallygrounded method for continuous covariate shift. Similarly, the subsequent work [34] developed amethod for continuous label shift, which allows for training the predictor with various models by
separately estimating the label probability through an online regression oracle. Another research [35]studied the change detection for continuous covariate shift. However, they only focused on identifying
differences between the current and initial offline distributions, and it is still unclear how to update
the model adaptively with the online data. Moreover, the method [35] can only detect at a given time
granularity, while ours can perform the model update at each round.

3

Adapting to Continuous Covariate Shift

In this section, we formulate the problem setup of continuous covariate shift and then introduce ourapproach based on the IWERM framework and online density ratio estimation.
3.1

Problem Setup

There are two stages in continuous covariate shift. The first one is the offline initialization stage,0
where the learner can collect a reasonable number of label data S0 = {xn , yn }N
n=1 from the initial
distribution D0 . Then, we come to the online adaptation stage, where the unlabeled testing data
arrive sequentially, and the underlying distributions can continuously shift. Consider a T -round
online adaptation. At each round t ∈ [T ] ≜ {1, . . . , T }, the learner will receive an unlabeled datasett
St = {xn }N
n=1 sampled from the underlying distribution Dt . Without loss of generality, we consider
Nt = 1. We have the following continuous covariate shift condition.
Assumption 1 (Continuous Covariate Shift). For all x ∈ X in the feature space, y ∈ Y in the label
space, and any t ∈ [T ], we have
Dt (y | x) = D0 (y | x) and rt∗ (x) = Dt (x)/D0 (x) ≤ B < ∞.
We note that there are emerging discussions on the necessity of covariate shift adaptation [7, 23].
When there are infinite number of training samples, a well-specified large and over-parametrized
model can be trained to approximate D0 (y | x) and there is no need to perform covariate shift
adaptation. However, in the finite-sample cases, one would prefer to train a model with constraint
complexity to ensure its generalization ability, which leads to model misspecification and covariate
shift adaptation is indeed necessary. In this paper, we study the continuous covariate shift under amisspecified hypothesis space W, which generalizes the standard one-step covariate shift problem(see [6] and reference therein) to the continuous shift case. Our goal is to train a sequence of model3

{wt }Tt=1 that are comparable with the best model wt∗ ∈ arg minw∈W Rt (w) in the hypothesis spaceat each time. Consequently, we take the following average excess risk as the performance measure:
!
T
T
X
1 X
T
∗
b t }t=1 ) ≜b t) −
RT ({w
Rt (w
Rt (wt ) .
(1)
T t=1
t=1
We end this part by listing several common notations used throughout the paper. We denote by
R = maxx∈X ∥x∥2 the maximum norm of the input and by D = maxw1 ,w2 ∈W ∥w1 − w2 ∥2 thediameter of the hypothesis space W ⊆ Rd . The constant G = maxx∈X ,y∈Y,w∈W ∥∇ℓ(w⊤ x, y)∥2 is
the maximum gradient norm and L = maxx∈X ,y∈Y,w∈W |ℓ(w⊤ x, y)| is the upper bound of the lossevalues. We use the O(·)-notation
to hide dependence on logarithmic factors of T .
3.2

Importance Weighted ERM for Continuous Shift

Our algorithm is based on the importance weighted learning. Consider the scenario where the learnerhas a predefined rbt : X → [0, B] as an estimation for the true density ratio rt∗ (x) = Dt (x)/D0 (x)for each time t ∈ [T ]. Then, the predictor can be trained by the importance weighted ERM method:b t = arg minw∈W Ex∼S0 [bw
rt (x)ℓ(w⊤ x, y)].

(2)

The averaged excess risk of the predictor is closely related to the quality of the density ratio estimator.
Proposition 1. For any δ ∈ (0, 1], with probability at least 1 − δ, IWERM (2) with the estimator√PTb t }Tt=1 ) ≤ 2 t=1 Ex∼S0 |brbt (x) ensures RT ({wrt (x) − rt∗ (x)| /T + O(log(T /δ)/ N0 ) .
√In above, the O(log T / N0 ) term measures the generalization gap of the predictor since it istrained over the empirical data S0 instead of D0 . Such a rate is tight up to a logarithmic factor in T .
Indeed, considering a stationary environment (i.e., Dt (x, y) = D1 (x, y)) and an exact density ratioPTb t /T
estimation (i.e., rbt (x) = rt∗ (x)), Proposition 1 indicates that the averaged
model wT = t=1 w√
enjoys an excess risk bound of R1 (wT ) − R1 (w1∗ ) ≤ O(log T / N0 ), matching the lower boundfor importance weighted learning [36, Proposition 2] up to an O(log T ) factor.
PT
rt (x) − rt∗ (x)|]/T term, which is the averaged estimation
Our main focus lies on the t=1 Ex∼S0 [|berror of the density ratio estimator rbt to the time-varying ground truth rt∗ . To minimize the estimation
error, a fundamental challenge comes from the unknown non-stationarity exhibited in the online
environments — how to select a right amount of historical data to reuse for each iteration? Intuitively,
if the environments change slowly, it would be preferable to reuse all historical data to construct
the density ratio estimator rbt . However, when distribution shifts occur frequently, earlier data could
be useless even potentially harmful. In such scenarios, training the density ratio estimator with themost recent data would be a more rational strategy. Furthermore, even if we could capture a roughly
stationary period, it remains unclear how to update the density ratio estimator in an online mannerwith guarantees. To our knowledge, online density ratio estimation is underexplored in the literature,even for the stationary testing streaming, not to mention the more challenging non-stationary setup.
3.3

Online Density Ratio Estimation

Here, we present a generic reduction of the online density ratio estimation problem to a dynamic
regret minimization problem. This novel perspective paves our way for designing an algorithm intackling the non-stationarity of environments, as mentioned in Section 3.2.
Bregman Divergence Density Ratio Matching. Our reduction is based on the Bregman divergencedensity ratio matching [14], a general framework that unifies various existing DRE methods. Specifically, in the framework, the discrepancy between the ground-truth density ratio function rt∗ and any
density ratio function r is measured by the expected Bregman divergence over D0 (x) defined byEBψ (rt∗ ∥r) = Ex∼D0 (x) [Bψ (rt∗ (x)∥r(x))] ,(3)

where Bψ (a∥b) ≜ ψ(a) − ψ(b) − ∂ψ(b)(a − b) is the Bregman divergence and ψ : dom ψ → R isthe associated divergence function. The expected Bregman divergence can be equally rewritten asψ ∗ψ
EBψ (rt∗ ∥r) = Lψt (r) − Lt (rt ), where Lt is the loss for the density ratio function defined by




(4)
Lψ
t (r) = Ex∼D0 (x) ∂ψ(r(x))r(x) − ψ(r(x)) − Ex∼Dt (x) ∂ψ(r(x)) .
4

As a consequence, one can train the density ratio estimator by minimizing the loss Lψ
t (r). We havert∗ ∈ arg minr∈H Lψ
(r)
when
H
is
the
set
of
all
measurable
functions.
t
The Bregman divergence-based density ratio matching framework takes various DRE methods as
special cases. For instance, choosing the divergence function as ψLS (t) = (t − 1)2 /2 yields LSIF [28]and KMM [26]; choosing ψLR (t) = t log t − (t + 1) log(t + 1) leads to the logistic regression
method [17]; and one can also recover the UKL [37] and KLLEP [25] with ψKL (t) = t log t − t.
Online DRE via Dynamic Regret Minimization. For a single-round density ratio estimation, theBregman divergence-based framework suggests to train the density ratio estimator rbt by minimizing∗
the gap Lψrt ) − Lψ
to derive an estimator sequence that performs well over time,t (bt (rt ). Therefore,PTPT∗rt ) − t=1 Lψ
we utilize the cumulative loss gap t=1 Lψt (bt (rt ) as a performance measure. Indeed,
the cumulative loss gap serves as an upper bound of the estimation error of density ratio estimators.
Proposition 2. Let ψ be a µ-strongly convex function.
estimator sequenceq ForPany density ratio

PT
PT
T
ψ
1
∗
T
∗
rt (x) − rt (x)|] ≤ 2rt ) − t=1 Lψ{brt }t=1 , we have T t=1 Ex∼D0 [|bt (rt ) /(µT ).
t=1 Lt (bProposition 2 indicates that it suffices to optimize the cumulative loss gap to perform online density
ratio estimation. Such a measure quantifies the performance difference between the online algorithm(that yields the estimated ratio sequence {brt }Tt=1 ) and a sequence of time-varying comparators (the∗ Ttrue ratio sequence {rt }t=1 ). This is exactly the dynamic regret in online learning literature [38, 39].
As such, we have reduced the online density ratio estimation to a problem of dynamic regretT∗ T
minimization over loss functions {Lψt }t=1 against comparators {rt }t=1 .
Since the loss function involves the expectation over the underlying distribution D0 (x) (see thedefinition in (4)), we need a counterpart result with respect to empirical data S0 .
Theorem 1. Let ψ be a µ-strongly convex function satisfying t∂ 3 ψ(t) ≤ 0 and ∂ 3 ψ(t) ≤ 0 for allt ∈ dom ψ. Let Hθ = {x 7→ h(x, θ) | θ ∈ Θ} be a hypothesis space of density ratio functions
parameterized by a finite-dimensional bounded set Θ ≜ {θ ∈ Rd | ∥θ∥2 ≤ S} with a certain linkfunction h : X × Θ 7→ R. Denote by [z]+ ≜ max{z, 0}. Then, for any density ratio estimator
rbt ∈ Hθ , the empirical estimation error is bounded by
v
!
u
√

 T
T
T
X
 ∗
 u
4 X eψ
d
log(T
/δ)
1X
ψ
t
∗
e t (rt ) + O√
LEx∼S0 (x) |rt (x) − rbt (x)| ≤L (brt ) −,
T t=1µT t=1 tµ N0+
t=1provided that h(x, θ) is bounded for any x ∈ X and θ ∈ Θ and Lipschitz continuous. In the above,eψLt (r) = Ex∼S0 [∂ψ(r(x))r(x) − ψ(r(x))] − Ex∼Dt (x) [∂ψ(r(x))](5)

is the empirical approximation of the expected loss Lψt using the data S0 .
Remark 1 (assumptions on ψ). We imposed certain assumptions on the divergence function ψ in
Theorem 1. One can check these conditions hold for commonly used ψ in density ratio estimation,
including all the divergence functions mentioned earlier (ψLR and ψKL are strongly convex when the
inputs are upper bounded, which can be satisfied with suitable choices of Hθ ). Section 4 will present
an example with ψLR to show how the conditions are satisfied.
◁
Theorem 1 shows that we can immediately obtain a sequence of high-quality density ratio estimators
Teψ{rt∗ }Tt=1 by minimizing the dynamic regret with respect to {Lt }t=1 ,XT
XT∗eψ ∗ Teψ
eψ
Regd
Lrt ) −
L(6)t (bt (rt ).
T ({Lt , rt }t=1 ) =
t=1

t=1

eψ
One caveat is that the second term of L
t in the definition (5) requires the knowledge of underlyingbψ
distribution Dt , which is unavailable. Empirically, we can only observe Lt defined below, building
upon the empirical observations St ∼ Dt ,bψLt (r) = Ex∼S0 [∂ψ(r(x))r(x) − ψ(r(x))] − Ex∼St [∂ψ(r(x))] .

(7)

That said, we need to design an online optimization process to minimize the dynamic regret (6)
∗
eψ
bψ T
defined over Lt based on the observed loss {Lt }t=1 . Since the time-varying comparator rt in (6)5ψbψ
is not the minimizer of the observed loss Lt (but rather the minimizer of the expected loss Ltψb t will lead to a high estimation error. In thedefined in (4)), directly minimizing the empirical loss L
next section, we introduce how to optimize the dynamic regret with the online ensemble framework
developed in recent studies of non-stationary online convex optimization [15, 16, 19].
Remark 2 (comparison with previous work). For the continuous label shift [18], the online ensemble
framework was employed to train the predictor wt . However, the previous method crucially relies on
the construction of an unbiased importance ratio estimator satisfying EDt [brt (x)] = rt∗ (x). Such afavorable property is hard to be satisfied in the covariate shift case. For example, in the Bregmanbψ
divergence density ratio matching framework, one can only observe an unbiased loss Lt (r), whose∗bψ
minimizer rbt = arg minr L
(r)
is
not
unbiased
to
the
true
value
r
in
general.
To
thisend, wettdisentangle the model training and importance weight estimation process in this paper. Our reduction
holds for the general Bergman divergence matching framework and thus can be initiated with different
models (not necessarily unbiased). It is possible to extend our framework to continuous label shiftproblem with other importance weight estimators besides the unbiased one [40] used in [18].
◁4Instantiation: Logistic Regression Modelbψ
When the loss function L
t is non-convex, it is generally intractable to conduct the online optimization,
regardless of minimizing the standard regret or the strengthened dynamic regret. Fortunately, theattained loss functions are convex or enjoy even stronger curvature with the properly chosen hypothesis
space and divergence function. In this section, we instantiate our framework with the logistic
regression model. A corresponding online DRE method is presented with dynamic regret guarantees.
Example 1 (Logistic regression model). Consider the functionψ = ψLR ≜ t log t − (t + 1) log(t + 1)
and the hypothesis space HθLR = {x 7→ exp −θ⊤ ϕ(x) | ∥θ∥2 ≤ S}. Here, ϕ : X 7→ Rd represents
a specific basis function with bounded norm ∥ϕ(x)∥2 ≤ R, which could be, for instance, the featurebψ
representation extractor from a deep neural network. Then, the loss function Lt (θ) as per (7) becomes


⊤
⊤
1
bψ
L
ES0 [log(1 + e−ϕ(x) θ )] + ESt [log(1 + eϕ(x) θ )] .
t (θ) =
2Let β = exp(SR). Then, the output of any density ratio function r ∈ Hθ is bounded by [1/β, β]. It
can be validated that ψLR is 1/(β + β 2 )-strongly convex and satisfies the condition ∂ 3 ψLR (z) ≤ 0
required by Theorem 1. Besides, the logistic loss is a 1/(1(1 + β))-expconcave function and R2 /2smooth function, presenting a favorable function properties for online convex optimization [13].
We note that the other two divergence functions ψLS and ψKL also exhibit desirable function properties
as ψLR . When choosing the hypothesis space Hθ = {x 7→ θ⊤ ϕ(x) | θ ∈ Θ}, the methods recover
the the uLISF [28] and UKL [37] density ratio estimators equipped with the generalized linear model.
Our analysis is also applicable in the two cases. More details are provided in Appendix D.7.
4.1

Density Ratio Estimation via Online Ensemble

In this part, we introduce our online ensemble method for online DRE with the logistic regression
model. Since the logistic regression loss is exp-concave, [19] shows that one can employ thefollowing-the-leading-history (FLH) method [41] to minimize the dynamic regret. However, a caveatb t is only an empirical estimation of Le t established on few online data St .
is that the observed loss L
The result of [19] only implies an expected bound for our problem.
We twisted the FLH algorithm to achieve a high probability bound. As shown in Figure 1, ouralgorithm maintains multiple base-learners, each learning over different intervals of the time-horizon,
and then employs a meta-learner to aggregate their predictions. We modify the meta-learner in FLH
from Hedge [42] to Adapt-ML-Prod [43], which ensures that the meta-learner can track each baselearner on the associated time interval with high probability. This strategy allows us to selectively
reuse the historical information to handle the non-stationary environments. We introduce ingredientsof our algorithms as follows. A more detailed algorithm descriptions and comparison with the
dynamic regret minimization literature in OCO can be found in Appendices D.1 and B respectively.
Base-learner. Our algorithm runs multiple base-learners that are active on different intervals of
the time horizon. Let C = {Ii = [si , ei ] ⊆ [T ]} be a interval set. Each base-learner Ei will only
6

t =1

t =2

t =3

t =4

t =5

t =6

t =7

t =8

t =9
update θbt,i by (8)

...
...
...
...

active base-learner
update pt,i by (9)meta-learner: θbt =

Pi∈At

pt,i θbt,iFigure 1: An illustration of our online ensemble method, where we employ a meta-learner to aggregate thepredictions from base-learners running over different intervals of the time horizon.

update and submit the model θbt,i to the meta-learner if t ∈ Ii . We use the online Newton step
(ONS) for the model update. For an exp-concave loss function, the method ensures the base-learner’s
model is comparable with the best fixed model over the interval up to a logarithmic factor in T [44].
Specifically, for each base-learner Ei running on the interval [si , ei ], ONS updates the model by

A 
b b
θbt+1,i = ΠΘt,i θbt,i − γA−1
(8)t,i ∇Lt (θt,i ) , ∀t ∈ [si , ei ] ,Ptb τ (θbτ,i )∇Lb s (θbτ,i )⊤ . In the above, the projection function
where the matrix At,i = λI + τ =si ∇LA

is defined as ΠΘt,i [θ1 ] = arg minθ∈Θ ∥θ − θ1 ∥At,i and Θ = {θ ∈ Rd | ∥θ∥2 ≤ S} is the parameter
space. The constant λ > 0 and η > 0 are the regularizer parameter and step size to be specified latter.
Meta-learner. As shown by Figure 1, the meta-learner’s role is to aggregate the models producedby the base-learners using a weighting scheme. At each iteration t, the meta-learner will maintain aweight pt,i for each “active” base-learner Ei , defined as the one whose associated interval contains t.
We update the weights for active base-learners based on the Adapt-ML-Prod method [43]. Specifically,for every active base-learner Ei , our algorithm maintains a “potential” vt,i ∈ R+ at iteration t, whichreflects the historical performance of the base-leaner until time t. Then, denoting by At the index setof the active base-learners at time t, their weights and the output model are obtained byX
pt,i ∝ εt−1,i vt−1,i for all i ∈ At and θbt+1 =
pt+1,i θbt+1,i .
(9)i∈At+1In the above εt,i > 0 is a step size that can be automatically tuned along the learning process. A
noteworthy ingredient of our algorithm is the construction of the potential vt,i , using the linearlizedb t (θbt ), θbt − θbt,i ⟩/(SR). Then, the generalization gap between Lb t and L
e t can be controlled
loss ⟨∇Lwith the negative term introduced by the exp-concave loss function. As a result, we can establish
a high probability bound to ensure that θbt is competitive with any base-learner’s model θbt,i on thecorresponding interval Ii . The detailed configurations of vt,i and εt,i are deferred to Appendix D.1.
Schedule of Intervals. As shown in Figure 1, we specify the intervals
with the geometS
ric covering scheme [41], where the interval set is defined as C = k∈N∪{0} Ck with Ck =
[i · 2k , (i + 1) · 2k − 1] | i ∈ N and i · 2k ≤ T . One can check that |C| is at most T , and thenumber of active base-learner is bounded as |At | ≤ ⌈log t⌉. Thus, we only need to maintain at mostO(log t) base-learners at time t. The intervals specified by the geometric covering are informative tocapture the non-stationarity of the environments, leading to the following dynamic regret guarantees.
4.2

Theoretical Guarantees

This part presents the theoretical guarantees of our method. The estimator rbt (x) = exp(−ϕ(x)⊤ θbt )established on θbt returned by the online ensemble method (9) achieves the dynamic regret guarantee.
Theorem 2. Assume the true density ratio rt∗ (x) = Dt (x)/D0 (x) is contained in the hypothesis
space as rt∗ ∈ HθLR ≜ {x 7→ exp(−ϕ(x)⊤ θ) | θ ∈ Θ} for any t ∈ [T ]. Then, with probability atleast 1 − δ, the dynamic regret of the density ratio estimator rbt (x) = exp(−ϕ(x)⊤ θbt ) is bounded by


2
1
∗ T
e
e
3 V 3 , 1} + T /N
Regd
({
L
,
r
}
)
≤
O
max{T
,
t
0
T
t t=1
T
when the parameters are set as γ = 3(1 + β) and λ = 1. In the above, VT =Dt−1 (x)∥1 measures the variation of input densities.
7PT

t=2 ∥Dt (x) −Theorem 2 imposes a realizable assumption for the true density ratio function such that rt∗ ∈ HθLR .
Such an assumption is required since we train the density ratio estimator on a given hypothesis space
while the true density ratio function could be arbitrary. We note that the realizability assumption isfrequently used in the analysis for the density ratio estimation problem [29, 45] and other related
topics, including active learning [46] and contextual bandit [47], where the density ratio (or density)
estimation is required. A possible direction to relax the assumption is to consider a richer function
class, e.g., a neural network or a non-parametric model. We leave this as a future work.
12

e 3 V 3 ) rate in our bound exhibits the same rate as the minimax optimal dynamic regret
The O(TTbound for the squared loss function with the noisy feedback [48]. Since the squared loss function
enjoys even stronger curvature than the logistic loss, our result is hard to be improved.
To achievethis fast-rate result, the key was to show a “squared” formulation O max{1, |I|VI2 } of the dynamic
regret for the base-leaner on each interval I (Lemma 3 in Appendix D.2). [19] achieved this by acomplicated analysis with the KKT condition to capture the structure of the comparators, while wegreatly simplified the analysis by exploiting the structure that the comparator θt∗ is essentially theminimizers of expected functions Lt (hence avoiding analyzing the KKT condition). Our analysis isapplicable to the case where the minimiers lie in the interior of the decision set and only requires thesmoothness of the loss function, which can be of independent interest in online convex optimization.
Averaged Excess Risk Bound for Continuous Covariate Shift. After obtaining the density ratio
estimator, we can train the predictive model by IWERM (2), leading to the following guarantee.
PT
Theorem 3. Under the same condition as Theorem 2 and letting VT = t=2 ∥Dt (x) − Dt−1 (x)∥1 ,running IWERM (2) with the estimated density ratio function rbt (x) = exp(−ϕ(x)⊤ θbt ) yields


1
1

e N − 2 + max T − 13 V 3 , T − 12 .
b t }Tt=1 ) ≤ O
RT ({w
0T
b t adapts to the environment with a converged average
Theorem 3 shows that our learned predictor w
excess risk when compared with the per-round best predictor wt∗ . By the discussion below Proposie −1/2 ) generalization error over the offline data S0 is hard to improve. We focus ontion 1, the O(N0the error in the online learning part. When the environment is nearly-stationary, i.e., VT ≤ O(T −1/2 ),e −1/2 +T −1/2 ) average excess risk, matching the same rate for the one-stepour bound implies an O(N0adaptation [49], even if the unlabeled data appear sequentially and the comparator wt∗ could change
1/3
over time. When the environments shift quickly, the O(T −1/3 VT ) rate still exhibits a diminishingerror for covariate shift adaptation once VT = o(T ). Indeed, our result has the same rate as that forcontinuous label shift [18] (with a slightly different definition of VT though). In Appendix D.6, we
1/3
further provide evidence to show that the O(max{T −1/3 VT , T −1/2 ) rate can hardly be improved,
even if one can receive labels of the testing stream after prediction.
Discussion on Assumptions. We end this section by a discussion on possible future directions torelax the linear model assumption for DRE and covariate shift assumption. Since our analysis is basedon the online convex optimization framework, we employed the (generalized) linear model for DRE
to ensure the convexity. To go beyond the linear model while still having theoretical guarantees, one
might extend the online ensemble framework to learn within the Reproducing Kernel Hilbert Space,leveraging advances in online kernel learning [50]. As for the step towards handling the general
distribution shift, it is possible to consider the sparse joint shift model [51] where both covariate andlabel distribution could shift. Besides, studying how to handle the joint distribution shift with few
online labeled data also presents an interesting future direction. Our research on the covariate shift
might serve as a basic step towards addressing more complex real-world distribution shifts.

5

Experiments

Setups. We generate continuous covariate shift by combining two fixed distributions with a timevarying mixture proportion αt ∈ [0, 1]. Specifically, given two fixed distributions D′ (x) and D′′ (x),we generate samples from Dt (x) = (1 − αt )D′ (x) + αt D′′ (x) at each round. The mixture proportion
αt shifts in four patterns: in Sin Shift and Squ Shift, αt changes periodically, following sine andsquare waves; in Lin Shift, the environment changes slowly from α1 = 1 to αT = 0 linearly over T
rounds, while in Ber Shift, the proportion αt flips quickly between 0 and 1 with a certain probability.
For parameterization in the Accous implementation, we set R by directly calculating the data norm
8

6.3 6.4 6.5 6.7 7.1 8.2 10.918.8 4.2 4.4 4.4 4.3 4.1 4.2 3.4
5.3 5.4 4.5 4.7 5.9 6.5 8.5 13.820.9 4.4 4.3 4.4 3.3 4.9 3.3
6.2 5.3 5.6 5.7 5.9 6.5 7.6 8.5 11.521.8 5.5 2.7 3.4 3.1 2.5

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14

25
Lin

Squ

Sin

Ber

Interval length (|Ii| = 2i)

0.6

1

0.55
0.8
0.5
0.45

0.6

0.4
0.4

0.35
0.3

0.2

Estimator Error

6.6 6.8 7.1 7.7 8.3 13.521.7 1.8 0.6 0.1 5.9 7.6 2.4 1.8 8.1

20.0
17.5
15.0
12.5
10.0
7.5
5.0
2.5

Average Error (%)

30

7.1 7.4 8.0 9.1 11.622.010.3 6.9 3.2 4.0 3.0 0.9 0.8 3.8 1.9

64

32

16

Distribution Change Interval (M)35

8.4 8.9 9.7 10.914.4 8.4 5.8 1.4 7.1 4.8 5.1 4.7 3.5 2.8 4.1

512 256 128

Average Error (%)

40

0.25
0.2
100

200

300

400

500

600

700

800

0
900

Time

Figure 2: Performance compari- Figure 3: Weight(%) heatmap of Figure 4: Average error and estison on the four kinds of shifts.
base-learners in squared shift.
mator loss in squared shift.

and set S = d/2 for all experiments. We repeat all experiments five times and evaluate the algorithms
by the average classification error over the unlabeled test stream for 10,000 rounds.
Contenders. We compare our method with six algorithms, which can be divided into three groups.
The first is a baseline approach that predicts directly using the model trained on initial offline data
(FIX). The second group consists of the one-step covariate shift adaptation methods that do notreuse historical information: DANN [52] handles the shifts by learning an invariant feature space, IWKMM [26], IW-KEIEP [25] and IW-uLSIF [28] equip the one-step IWERM method (2) with different
density ratio estimators. The third is OLRE, which serves as an ablation study for our Accous
algorithm. The OLRE algorithm estimates the density ratio by running an ONS (8) with all historical
data and then performs the IWERM at each round. All algorithm parameters are set by their default.
In the following, we focus on the results of empirical studies. Detailed configurations for the datasets,
simulated shifts, and setups for the Accous algorithm and contenders can be found in Appendix A.
5.1

Illustrations on Synthetic Data

Average Classification Error Comparison. We summarize the comparison results on synthetic
data with different types of covariate shifts with all contenders in Figure 2. The Accous algorithm
outperforms almost all other methods in the four shift patterns. We observe that the offline model FIX
performs poorly compared to the online methods. The best result of DANN, IW-KMM, IW-KLIEPand IW-uLSIF is called One-Step for comparison. By reusing historical data, the Accous algorithm
achieves a lower average classification error compared to these One-Step competitors that use only
one round of online data. The Accous also outperforms OLRE, especially when the environmentschange relatively quickly (Squ, Sin, and Ber). These results justify that performing density ratio
estimation with either one round of online data or all historical data is inappropriate for non-stationary
environments. Moreover, when the number of unlabeled data per round changes from Nt = 5
to Nt = 1, the one-step methods suffer severe performance degradation. In contrast, our Accous
algorithm still performs relatively well, highlighting the need for selective reuse of historical data.
Effectiveness of Online Density Ratio Estimation. We then take a closer look at the core componentof our algorithm, the online density ratio estimator, in the Squ shift where the distribution of theonline data shifts every M rounds. Figure 3 shows the weight assignment for base-learners withdifferent interval lengths, averaged over their active period. The result shows that our meta-learner
successfully assigns the largest weight to the base-learner whose interval length matches the switching
period M , which ensures that the right amount of historical data is reused. Figure 4 also shows theaverage error of Accous and OLRE over 1, 000 iterations. The results show that Accous (red line)
can quickly adapt to new distributions after the covariate shift occurs. In contrast, OLRE (blue line)
struggles because it uses all the historical data generated from different distributions. Furthermore,b t (θbt ) by the green line, which is quickly minimized
we denote the loss of the density ratio estimator L
after the covariate shift occurs. The similar tendency between the average error (red line) and theloss of the estimator (green line) confirms our theoretical finding in Theorem 1: we can minimize theexcess risk of the predictor by minimizing the dynamic regret of the density ratio estimator.
5.2

Comparison on Real-world Data

Performance Comparison on Benchmarks. In Table 1 and Table 2, we present the results of
the average classification error comparison with the contenders on four benchmark datasets. Our
9Table 1: Average classification error (%) of different algorithms on various real-world datasets with Lin
and Ber shifts. We report the mean accuracy and standard deviation over five runs. The best algorithms areemphasized in bold. The number of data in each round is set as Nt = 5.
Lin

Diabetes
Breast
MNIST-SVHN
CIFAR10-CINIC10

Ber

FIX

DANN

KMM

KLIEP

uLSIF

LR

OLRE

Accous

FIX

DANN

KMM

KLIEP

uLSIF

LR

OLRE

Accous

39.55
±6.58
10.25
± 3.00
8.19
±1.39
29.24
±0.51

46.88
±7.09
16.03
± 4.79
13.37
±0.97
31.59
±2.18

43.67
±6.23
8.26
± 0.31
6.24
±1.52
29.60
±1.43

37.85
±6.15
6.90
± 0.77
7.52
±0.33
29.33
±0.97

38.32
±6.00
11.07
± 2.34
6.38
±0.76
29.18
±1.50

39.17
±7.10
7.33
± 1.59
6.81
±0.84
29.75
±0.48

37.40
±6.56
3.59
± 0.25
6.65
±0.93
28.17
±1.24

35.92
±6.44
4.60
± 0.16
5.93
±0.72
27.80
±0.83

40.90
±6.90
11.68
± 4.16
11.33
±0.83
33.16
±0.94

49.03
±6.54
15.34
± 5.23
16.54
±1.10
42.50
±1.29

46.80
±6.67
10.12
± 0.60
10.81
±0.63
33.49
±0.83

42.98
±6.52
7.13
± 1.30
11.50
±0.97
31.44
±1.57

43.22
±6.78
9.01
± 1.91
12.19
±1.20
32.98
±0.67

40.28
±6.25
8.14
± 2.02
11.18
±0.50
33.03
±1.26

42.05
±6.22
5.27
± 0.44
13.42
±0.49
31.22
±1.01

39.14
±4.51
3.89
± 0.37
10.31
±0.51
30.69
±0.62

Table 2: Average classification error (%) of different algorithms on various real-world datasets with Squ
and Sin shifts. We report the mean accuracy and standard deviation over five runs. The best algorithms areemphasized in bold. The number of data in each round is set as Nt = 5.
Squ

Diabetes
Breast
MNIST-SVHN
CIFAR10-CINIC10

Sin

FIX

DANN

KMM

KLIEP

uLSIF

LR

OLRE

Accous

FIX

DANN

KMM

KLIEP

uLSIF

LR

OLRE

Accous

37.85
±6.67
9.92
±0.23
11.28
±1.22
34.75
±0.29

45.69
±6.51
12.54
±1.55
13.98
±0.97
38.13
±1.58

43.57
±6.97
8.53
±0.30
10.21
±0.53
34.08
±0.85

37.00
±6.80
6.52
±0.91
11.01
±1.41
34.27
±1.24

38.03
±6.45
7.92
±1.39
10.49
±0.67
33.80
±1.73

37.83
±5.87
7.39
±0.83
11.72
±1.89
33.56
±1.59

39.17
±6.67
4.67
±0.23
13.68
±0.59
32.61
±0.69

36.45
±5.36
3.68
±0.58
9.54
±0.40
31.08
±0.60

37.30
±7.02
8.68
±0.56
9.98
±1.29
36.44
±0.47

44.96
±6.79
14.34
±1.78
16.73
±1.03
42.53
±1.09

41.38
±6.46
7.18
±0.39
8.91
±0.68
37.00
±0.67

36.40
±6.93
6.11
±0.92
9.41
±1.01
36.58
±1.52

37.23
±6.88
9.59
±1.02
9.30
±0.53
36.53
±1.02

36.99
±6.21
7.50
±0.51
9.00
±1.77
36.21
±1.31

38.58
±7.45
4.16
±0.42
15.38
±0.88
38.24
±1.86

35.17
±6.44
3.29
±0.43
7.74
±0.74
35.42
±1.22

algorithm Accous outperforms other contenders in the four types of shifts, whether in a slowly
evolving environment (Lin) or in relatively non-stationary environments (Squ,Sin,Ber). The
OLRE algorithm uses all historical data, but does not necessarily outperform contenders that use only
current data, suggesting that we should selectively reuse the historical data.
Average Error (%)

Case Study on a Real-life Application. We conducted empirical
0.45
studies on a real-world application using the yearbook dataset [53],0.4which contains high school front-facing yearbook photos from 1930
to 2013. The photo distributions shift along years due to the changing
0.35
social norms, fashion styles, and population demographics. Consequently, it is essential to adapt to these changes. We used photos
0.3
before 1945 as the offline labeled dataset and generated the online
unlabeled data stream using photos after 1945, arranged chrono0.25
100 200 300 400 500 600 700 800 900
logically with 10 photos per round. Figure 5 shows the averageTimeerror curves. The OLER algorithm performs well when the distriFigure
5:
Averageerror on Yearbutions do not change too much, but the method starts to suffer frombook
datasetwithreal-life
covarilarge prediction errors after 400 iterations. In contrast, our method
ate shift.
outperforms all competitors by selectively reusing historical data.

6

Conclusion

This paper initiated the study of the continuous covariate shift with the goal of minimizing the
cumulative excess risk of the predictors. We showed that the importance-weighted ERM method
works effectively given high-quality density ratio estimators (Proposition 1), whose estimation canbe cast as a dynamic regret minimization problem (Theorem 1). Instantiating with the logistic
regression density ratio estimation model, we proposed a novel online that can adaptively reuse
the historical data and enjoy a tight dynamic regret bound (Theorem 2). The regret bound finallye −1/3 V 1/3 ) averaged excess risk guarantee for the predictor. Experiments validate theimplies an O(TT
effectiveness of the proposed method and demonstrate the need for selective reuse of historical data.

Acknowledgments
YJZ and MS were supported by the Institute for AI and Beyond, UTokyo. YJZ was also supported byTodai Fellowship. Peng Zhao was supported by NSFC (62206125) and JiangsuSF (BK20220776).
The authors thank Yong Bai and Yu-Yang Qian for their assistance in the experiments.
10References[1] Joaquin Quinonero-Candela, Masashi Sugiyama, Anton Schwaighofer, and Neil D Lawrence.
Dataset Shift in Machine Learning. MIT Press, 2008.
[2] Yoshua Bengio, Yann Lecun, and Geoffrey Hinton. Deep learning for AI. Communications of
the ACM, 64(7):58–65, 2021.
[3] Thomas G. Dietterich. Steps toward robust artificial intelligence. AI Magazine, 38(3):3–24,
2017.
[4] Zhi-Hua Zhou. Open-environment machine learning. National Science Review, 9(8), 07 2022.
nwac123.
[5] Hidetoshi Shimodaira. Improving predictive inference under covariate shift by weighting the
log-likelihood function. Journal of Statistical Planning and Inference, 90(2):227–244, 2000.
[6] Masashi Sugiyama and Motoaki Kawanabe. Machine Learning in Non-Stationary Environments- Introduction to Covariate Shift Adaptation. MIT Press, 2012.
[7] Samory Kpotufe and Guillaume Martinet. Marginal singularity and the benefits of labels incovariate-shift. The Annals of Statistics, 49(6):3299–3323, 2021.
[8] Yan Li, Hiroyuki Kambara, Yasuharu Koike, and Masashi Sugiyama. Application of covariate
shift adaptation techniques in brain-computer interfaces. IEEE Transactions on Biomedical
Engineering, 57(6):1318–1324, 2010.
[9] Makoto Yamada, Masashi Sugiyama, and Tomoko Matsui. Semi-supervised speaker identification under covariate shift. Signal Processing, 90(8):2353–2361, 2010.
[10] Ali Hassan, Robert I. Damper, and Mahesan Niranjan. On acoustic emotion recognition: Compensating for covariate shift. IEEE Transactions on Audio, Speech and Language Processing,21(7):1458–1468, 2013.
[11] Jiayuan Huang, Alexander J. Smola, Arthur Gretton, Karsten M. Borgwardt, and Bernhard
Schölkopf. Correcting sample selection bias by unlabeled data. In Advances in Neural Information Processing Systems 19 (NIPS), pages 601–608, 2006.
[12] Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori. Density Ratio Estimation in Machine
Learning. Cambridge University Press, 2012.
[13] Elad Hazan. Introduction to Online Convex Optimization. Foundations and Trends in Optimization, 2(3-4):157–325, 2016.
[14] Masashi Sugiyama, Taiji Suzuki, and Takafumi Kanamori. Density-ratio matching under the
Bregman divergence: a unified framework of density-ratio estimation. Annals of the Institute of
Statistical Mathematics, 64(5):1009–1044, 2012.
[15] Lijun Zhang, Shiyin Lu, and Zhi-Hua Zhou. Adaptive online learning in dynamic environments.
In Advances in Neural Information Processing Systems 31 (NeurIPS), pages 1330–1340, 2018.
[16] Peng Zhao, Yu-Jie Zhang, Lijun Zhang, and Zhi-Hua Zhou. Dynamic regret of convex andsmooth functions. In Advances in Neural Information Processing Systems 33 (NeurIPS), pages
12510–12520, 2020.
[17] Steffen Bickel, Michael Brückner, and Tobias Scheffer. Discriminative learning under covariate
shift. Journal of Machine Learning Research, 10:2137–2155, 2009.
[18] Yong Bai, Yu-Jie Zhang, Peng Zhao, Masashi Sugiyama, and Zhi-Hua Zhou. Adapting to online
label shift with provable guarantees. In Advances in Neural Information Processing Systems 35(NeurIPS), pages 29960–29974, 2022.
[19] Dheeraj Baby and Yu-Xiang Wang. Optimal dynamic regret in exp-concave online learning. In
Proceedings of the 34th Conference on Learning Theory (COLT), pages 359–409, 2021.
[20] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex Optimization. Cambridge
University Press, 2004.
[21] Makoto Yamada, Taiji Suzuki, Takafumi Kanamori, Hirotaka Hachiya, and Masashi Sugiyama.
Relative density-ratio estimation for robust distribution comparison. Neural Computation,
25(5):1324–1370, 2013.
11

[22] Masashi Sugiyama, Matthias Krauledat, and Klaus-Robert Müller. Covariate shift adaptation byimportance weighted cross validation. Journal of Machine Learning Research, 8:985–1005,
2007.
[23] Junfeng Wen, Chun-Nam Yu, and Russell Greiner. Robust learning under uncertain test
distributions: Relating covariate shift to model misspecification. In Proceedings of the 31th
International Conference on Machine Learning (ICML), pages 631–639, 2014.
[24] Tongtong Fang, Nan Lu, Gang Niu, and Masashi Sugiyama. Rethinking importance weighting
for deep learning under distribution shift. In Advances in Neural Information Processing
Systems 33 (NeurIPS), page 11996–12007, 2020.
[25] Masashi Sugiyama, Taiji Suzuki, Shinichi Nakajima, Hisashi Kashima, Paul von Bünau, andMotoaki Kawanabe. Direct importance estimation for covariate shift adaptation. Annals of theInstitute of Statistical Mathematics, 60(4):699–746, 2008.
[26] Arthur Gretton, Alex Smola, Jiayuan Huang, Marcel Schmittfull, Karsten Borgwardt, andBernhard Schölkopf. Covariate shift by kernel mean matching. Dataset Shift in Machine
Learning, 3(4):5, 2009.
[27] Masashi Sugiyama. Superfast-trainable multi-class probabilistic classifier by least-squares
posterior fitting. IEICE Transactions on Information and Systems, 93-D(10):2690–2701, 2010.
[28] Takafumi Kanamori, Shohei Hido, and Masashi Sugiyama. A least-squares approach to direct
importance estimation. Journal of Machine Learning Research, 10:1391–1445, 2009.
[29] Masahiro Kato and Takeshi Teshima. Non-negative Bregman divergence minimization for deep
direct density ratio estimation. In Marina Meila and Tong Zhang, editors, Proceedings of the
38th International Conference on Machine Learning (ICML), pages 5320–5333, 2021.
[30] M. C. du Plessis, H. Shiino, and M. Sugiyama. Online direct density-ratio estimation applied toinlier-based outlier detection. Neural Computation, 27(9):1899–1914, 2015.
[31] Yu Chen, Song Liu, Tom Diethe, and Peter A. Flach. Continual density ratio estimation in an
online setting. ArXiv preprint, arXiv:2103.05276, 2021.
[32] Rasool Fakoor, Jonas Mueller, Zachary C. Lipton, Pratik Chaudhari, and Alexander J. Smola.
Time-varying propensity score to bridge the gap between the past and present. ArXiv preprint,
arXiv: 2210.01422, 2022.
[33] Ruihan Wu, Chuan Guo, Yi Su, and Kilian Q. Weinberger. Online adaptation to label distribution
shift. In Advances in Neural Information Processing Systems 34 (NeurIPS), pages 11340–11351,
2021.
[34] Dheeraj Baby, Saurabh Garg, Tzu-Ching Yen, Sivaraman Balakrishnan, Zachary Chase Lipton,
and Yu-Xiang Wang. Online label shift: Optimal dynamic regret meets practical algorithms.
ArXiv preprint, arXiv:2103.05276, 2023.
[35] Sooyong Jang, Sangdon Park, Insup Lee, and Osbert Bastani. Sequential covariate shiftdetection using classifier two-sample tests. In Proceedings of the 39th International Conference
on Machine Learning (ICML), pages 9845–9880, 2022.
[36] Corinna Cortes, Yishay Mansour, and Mehryar Mohri. Learning bounds for importanceweighting. In Advances in Neural Information Processing Systems 23 (NIPS), pages 442–450,
2010.
[37] XuanLong Nguyen, Martin J. Wainwright, and Michael I. Jordan. Estimating divergence
functionals and the likelihood ratio by penalized convex risk minimization. In Advances in
Neural Information Processing Systems 20 (NeurIPS), pages 1089–1096, 2007.
[38] Martin Zinkevich. Online convex programming and generalized infinitesimal gradient ascent. In
Proceedings of the 20th International Conference on Machine Learning (ICML), pages 928–936,
2003.
[39] Omar Besbes, Yonatan Gur, and Assaf J. Zeevi. Non-stationary stochastic optimization. Operations Research, 63(5):1227–1244, 2015.
[40] Zachary C. Lipton, Yu-Xiang Wang, and Alexander J. Smola. Detecting and correcting forlabel shift with black box predictors. In Proceedings of the 35th International Conference on
Machine Learning (ICML), pages 3128–3136, 2018.
12

[41] Elad Hazan and C. Seshadhri. Efficient learning algorithms for changing environments. In
Proceedings of the 26th International Conference on Machine Learning (ICML), pages 393–400,
2009.
[42] Yoav Freund and Robert E. Schapire. A decision-theoretic generalization of on-line learningand an application to boosting. Journal of Computer and System Sciences, 55(1):119–139,
1997.
[43] Pierre Gaillard, Gilles Stoltz, and Tim van Erven. A second-order bound with excess losses. In
Proceedings of The 27th Conference on Learning Theory (COLT), pages 176–196, 2014.
[44] Elad Hazan, Amit Agarwal, and Satyen Kale. Logarithmic regret algorithms for online convex
optimization. Machine Learning, 69(2-3):169–192, 2007.
[45] Takafumi Kanamori, Taiji Suzuki, and Masashi Sugiyama. Statistical analysis of kernel-based
least-squares density-ratio estimation. Machine Learning, 86(3):335–367, 2012.
[46] Akshay Krishnamurthy, Alekh Agarwal, Tzu-Kuo Huang, Hal Daumé III, and John Langford.
Active learning for cost-sensitive classification. In Proceedings of the 34th International
Conference on Machine Learning (ICML), pages 1915–1924, 2017.
[47] Dylan J. Foster and Alexander Rakhlin. Beyond UCB: optimal and efficient contextual bandits
with regression oracles. In Proceedings of the 37th International Conference on Machine
Learning (ICML), pages 3199–3210, 2020.
[48] Dheeraj Baby and Yu-Xiang Wang. Online forecasting of total-variation-bounded sequences. In
Advances in Neural Information Processing Systems 32 (NeurIPS), pages 11071–11081, 2019.
[49] Yaoliang Yu and Csaba Szepesvári. Analysis of kernel mean matching under covariate shift.
In Proceedings of the 29th International Conference on Machine Learning (ICML), pages
1147–1154, 2012.
[50] Daniele Calandriello, Alessandro Lazaric, and Michal Valko. Efficient second-order online
kernel learning with adaptive embedding. In Advances in Neural Information Processing
Systems 30 (NeurIPS), pages 6140–6150, 2017.
[51] Lingjiao Chen, Matei Zaharia, and James Y. Zou. Estimating and explaining model performance
when both covariates and labels shift. In Advances in Neural Information Processing Systems
35 (NeurIPS), pages 11467–11479, 2022.
[52] Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François
Laviolette, Mario Marchand, and Victor S. Lempitsky. Domain-adversarial training of neuralnetworks. Journal of Machine Learning Research, 17:59:1–59:35, 2016.
[53] Huaxiu Yao, Caroline Choi, Bochuan Cao, Yoonho Lee, Pang Wei Koh, and Chelsea Finn.
Wild-time: A benchmark of in-the-wild distribution shift over time. In Advances in Neural
Information Processing Systems 35 (NeurIPS) Datasets and Benchmarks Track, pages 10309–
10324, 2022.
[54] Shiry Ginosar, Kate Rakelly, Sarah M Sachs, Brian Yin, Crystal Lee, Philipp Krähenbühl, andAlexei A Efros. A century of portraits: A visual historical record of american high school
yearbooks. IEEE Transactions on Computational Imaging, 3(3):421–431, 2017.
[55] András György and Csaba Szepesvári. Shifting regret, mirror descent, and matrices. In
Proceedings of the 33nd International Conference on Machine Learning (ICML), pages 2943–
2951, 2016.
[56] Eric C. Hall and Rebecca Willett. Dynamical models and tracking regret in online convex
programming. In Proceedings of the 30th International Conference on Machine Learning(ICML), pages 579–587, 2013.
[57] Ali Jadbabaie, Alexander Rakhlin, Shahin Shahrampour, and Karthik Sridharan. Online
optimization : Competing with dynamic comparators. In Proceedings of the 18th International
Conference on Artificial Intelligence and Statistics (AISTATS), pages 398–406, 2015.
[58] Aryan Mokhtari, Shahin Shahrampour, Ali Jadbabaie, and Alejandro Ribeiro. Online optimization in dynamic environments: Improved regret rates for strongly convex problems. In
Proceedings of the 55th IEEE Conference on Decision and Control (CDC), pages 7195–7201,
2016.
13

[59] Tianbao Yang, Lijun Zhang, Rong Jin, and Jinfeng Yi. Tracking slowly moving clairvoyant:
Optimal dynamic regret of online learning with true and noisy gradient. In Proceedings of the
33rd International Conference on Machine Learning (ICML), pages 449–457, 2016.
[60] Lijun Zhang, Tianbao Yang, Rong Jin, and Zhi-Hua Zhou. Dynamic regret of strongly adaptive
methods. In Proceedings of the 35th International Conference on Machine Learning (ICML),
pages 5877–5886, 2018.
[61] Yu-Jie Zhang, Peng Zhao, and Zhi-Hua Zhou. A simple online algorithm for competing withdynamic comparators. In Proceedings of the 36th Conference on Uncertainty in Artificial
Intelligence (UAI), pages 390–399, 2020.
[62] Dheeraj Baby and Yu-Xiang Wang. Adaptive online estimation of piecewise polynomial trends.
In Advances in Neural Information Processing Systems 33 (NeurIPS), pages 20462–20472,
2020.
[63] Peng Zhao and Lijun Zhang. Improved analysis for dynamic regret of strongly convex andsmooth functions. In Proceedings of the 3rd Conference on Learning for Dynamics and Control(L4DC), pages 48–59, 2021.
[64] Zhi-Hua Zhou. Ensemble Methods: Foundations and Algorithms. Chapman & Hall/CRC Press,
2012.
[65] Shai Shalev-Shwartz and Shai Ben-David. Understanding Machine Learning: From Theory to
Algorithms. Cambridge University Press, 2014.
[66] Gilles Pisier. The Volume of Convex Bodies and Banach Space Geometry, volume 94. Cambridge
University Press, 1999.
[67] Peng Zhao, Yu-Jie Zhang, Lijun Zhang, and Zhi-Hua Zhou. Adaptivity and non-stationarity:
Problem-dependent dynamic regret for online convex optimization. ArXiv preprint, arXiv:
2002.07729, 2021.
[68] Mehrdad Mahdavi, Lijun Zhang, and Rong Jin. Lower and upper bounds on the generalizationof stochastic exponentially concave optimization. In Proceedings of The 28th Conference on
Learning Theory (COLT), pages 1305–1320, 2015.
[69] Joel A. Tropp. An introduction to matrix concentration inequalities. Foundations and Trends in
Machine Learning, 8(1-2):1–230, 2015.
[70] Kwang-Sung Jun, Francesco Orabona, Stephen Wright, and Rebecca Willett. Improved strongly
adaptive online learning using coin betting. In Proceedings of the 20th International Conferenceon Artificial Intelligence and Statistics (AISTATS), pages 943–951, 2017.
[71] Yoav Freund, Robert E. Schapire, Yoram Singer, and Manfred K. Warmuth. Using and combining predictors that specialize. In Proceedings of the 29th Annual ACM Symposium on the
Theory of Computing (STOC), pages 334–343, 1997.
[72] Nicolò Cesa-Bianchi and Gábor Lugosi. Prediction, Learning, and Games. Cambridge
University Press, 2006.
[73] Haipeng Luo and Robert E. Schapire. Achieving all with no parameters: AdaNormalHedge. In
Proceedings of the 28th Annual Conference Computational Learning Theory (COLT), pages
1286–1304, 2015.
[74] Dmitry Adamskiy, Wouter M. Koolen, Alexey V. Chernov, and Vladimir Vovk. A closer look at
adaptive regret. Journal of Machine Learning Research, 17:23:1–23:21, 2016.
[75] Lijun Zhang, Shiyin Lu, and Tianbao Yang. Minimizing dynamic regret and adaptive regret
simultaneously. In Proceedings of the 23rd International Conference on Artificial Intelligenceand Statistics (AISTATS), pages 309–319, 2020.

14

A

Omitted Details for Experiments

In this section, we first provide additional details that were omitted in section 5. We then provide thedescriptions of our experimental setups, which include details of the datasets used, simulations ofcontinuous covariate shifts, descriptions of the competing algorithms, and the parameter configurationof our proposed Accous algorithm. After that, we supplement the detailed numerical results obtained
from the synthetic experiments.
A.1

Experiment Setup

Continuous Covariate Shift Simulation. To simulate the continuous covariate shift, we generate
offline and online samples from different marginal distributions and annotate them with the samelabeling function. Given two different marginal distributions D′ (x) and D′′ (x), the offline dataset
S0 is sampled from the distribution D0 (x) = (1 − α0 )D′ (x) + α0 D′′ (x), where α0 ∈ [0, 1] is a
mixture proportion coefficient. The online dataset St received at each iteration t is sampled from the
distribution Dt (x) = (1 − αt )D′ (x) + αt D′′ (x), where αt ∈ [0, 1] is the time-varying coefficient.
By changing the coefficient αt for t ≥ 1, we can simulate different types of continuous covariate
shifts with different shift patterns.
• Linear Shift (Lin): We set αt = t/T to simulate a gradually changing environment.
• Square Shift (Squ): In this case, the parameter αt switches√between 1 and 0 every
M rounds, where 2M is the periodic length. We set M = Θ( T ) to simulate a fastchanging environment with periodic patterns.
iπ• Sine Shift (Sin): We set αt = sin M, where i = t mod M and M is the periodic√
length. We set M = Θ( T ) to simulate a fast-changing environment with periodic patterns.

• Bernoulli Shift (Ber): In this case, we keep the αt =√αt−1 with probabilityp ∈ [0, 1] and otherwise set αt = 1 − αt−1 . We set p = Θ(1/ T ) to simulate a fastchanging environment but without periodic patterns.
Datasets. For the synthetic dataset, we generate D′ (x) and D′′ (x) using the multinomial Gaussian
distribution, wherein each class is generated from a Gaussian distribution. For D′ (x), the positiveclass data is generated from N (x; µp1 , Σ), and the negative class data is generated from N (x; µn1 , Σ).
Here, µp1 is set to −1.2 in R12 , and µn1 is set to −0.8 in R12 . As for D′′ (x), its positive class data is
generated from N (x; µp2 , Σ), and the negative class data is generated from N (x; µn2 , Σ). Here, µp2 is
specified as 1.2 in R12 , and µn2 as 0.8 in R12 . The covariance matrix is set as Σ = Id in R12×12 ,and the prior probability for the positive class data is fixed at 0.5. Besides, we generated offline datawith α0 = 0.9 and introduced a time-varying coefficient, αt , to simulate continuous covariate shifts.
The benchmark datasets are generated from real-world datasets.
• Diabetes (Dia) is a tabular UCI dataset. Each sample is an 8-dimensional vector with abinary label. We split the dataset into two subsets and choose the split such that the classifier
trained on one dataset generalizes poorly to the other. We choose one of them as the offline
dataset and combine them with a time-varying coefficient to generate the online data.
• Breast (Bre) is a tabular UCI dataset. Each sample is a 9-dimensional vector with a binary
label. This dataset is generated similarly to the Diabetes dataset.
• MNIST-SVHN (M-S) contains the digital images collected from handwriting or natural scene
images. A feature extractor pre-trained on this dataset produces a 512-dimensional vector
with a 10-class label for each sample. We choose 10% MNIST data and 90% SVHN data as
the offline dataset and generate the online data by mixing the MNIST and SVHN datasetswith a time-varying proportion coefficient αt .
• CIFAR10-CINIC10 (C-C) contains real-life subjective images collected from different
sources. This dataset is generated similarly to the MNIST-SVHN dataset.
We also conduct experiments on a real-world dataset where the underlying data distributions exhibit
a continuous shift.
15• Yearbook [53]: This dataset contains the 37,921 frontal-facing American high school
yearbook photos from 1930 to 2013 [54]. Each photo is a 32 × 32 × 1 grayscale image
associated with a binary label representing the student’s gender. This dataset captures thereal-life changing of social norms, fashion styles, and population demographics. The offline
data is generated with photos before 1945, the online data is generated from 1946, and thedata size is 10 per round.
Contenders. We compare our proposed approach with six contenders, including:• FIX is the baseline method that predicts with the classifier that are pre-trained on offline data
without any updates in the online adaptation stage.
• is a general domain adaptation algorithm that mitigates the discrepancy between the training
and test distributions by learning an invariant feature representation. We adapt this method
to the online setting by running DANN at each round t ∈ [T ] based on the offline data S0
and the online data St available at that round.
• IW-KMM is an importance weighting learning method that first estimates the density ratiosby the Kernel Mean Matching (KMM) method [26] and then performs the IWERM to train
the model. At each iteration t ∈ [T ], the density ratios are estimated with the offline dataset
S0 and the online dataset St by the KMM. Then the model is trained by running the IWERM
method on the offline data S0 .
• IW-KLIEP is also an importance weighting learning method, where the density ratios are
estimated by the KLIEP method [25]. At each iteration t ∈ [T ], after the density ratios havebeen estimated by KLIEP, the learner performs the IWERM over the offline data S0 .
• IW-uLSIF is an importance weighting learning method where the density ratio estimatedby the uLSIF method is used [28]. At each iteration t ∈ [T ], the learner first estimates thedensity ratios by the uLSIF and then performs the IWERM on the offline data S0 .
• LR serves as an ablation study for the proposed Accous algorithm. At each iteration t ∈ [T ],the learner first estimates the density ratios by an ONS with a logistic regression model andthen performs the IWERM on the offline data S0 .
• OLRE also serves as an ablation study for the proposed Accous algorithm. It estimates thedensity ratios in an on-line manner by running an ONS with a logistic regression model withall the historical data and then performing the IWERM each round.
Accous Setup. We implement the Accous algorithm using a logistic regression model. Before
running the Accous algorithm, two parameters must be defined. They are the maximum norm of
the feature R and the maximum norm of the density functions S. We set R by directly calculating
the data norm. The covariate shift level determines the parameter S, which is difficult to predict inadvance. We found that setting S = d/2 is a good choice in the experiments. For stable performance,
we set a threshold of 100 and truncate the large weights after the density ratio estimation in eachround. We also ignore the base models with interval lengths 1, 2 in the density ratio estimation.
For the synthetic, diabetes, and breast tabular datasets, the linear model is used. For the MNISTSVHN, CIFAR10-CINIC10, and Yearbook datasets, the deep model is used for density ratio estimation and predictive model training. The deep models consist of a backbone and a linear layer. In all
experiments, only the linear layer of the neural network is updated, while the backbone parameters
trained with offline data remain fixed. Backbone settings for different datasets are presented below.
• MNIST-SVHN and CIFAR10-CINIC10: The backbone is a pre-trained ResNet34 backbone
from torchvision with its weights initialized by training on ImageNet. Both the weight
estimator and the classifier share the same backbone and its corresponding parameters.
• Yearbook: Our network adopts a similar structure and updating process as in the MNISTSVHN and CIFAR10-CINIC10 experiments, following the configuration described in [53].
The backbone of our network consists of a CNN model with the default settings.
Computational Resources. We run experiments with two Xeon Gold 6248R processors (24 cores,
3.0GHz base, 4.0GHz boost), eight Tesla V100S GPUs (32GB video memory each), 768GB RAM,
all managed by the Ubuntu 20.04 operating system.
16

Table 3: Average error (%) over the entire time horizons for different algorithms under various simulated shifts.
The best algorithms are emphasized in bold.
Nt = 1

Lin
Squ
Sin
Ber

A.2

Nt = 5

FIX

DANN

KMM

KLIEP

uLSIF

LR

OLRE

Accous

FIX

DANN

KMM

KLIEP

uLSIF

LR

OLRE

Accous

32.85
±0.43
33.19
±0.56
32.80
±0.75
34.11
±0.45

46.93
±2.40
40.59
±3.28
42.65
±4.49
46.20
±4.12

37.97
±0.56
35.50
±0.33
36.25
±0.37
36.49
±0.42

36.26
±0.82
35.13
±0.61
35.83
±0.65
37.07
±0.67

38.34
±0.40
36.44
±0.79
37.01
±0.38
36.28
±0.36

36.51
±1.12
36.21
±0.46
37.42
±0.73
37.42
±0.85

27.13
±0.74
34.82
±0.30
33.15
±0.50
34.91
±0.47

27.69
±0.49
29.30
±0.82
31.08
±0.58
29.55
±0.46

32.59
±0.53
33.37
±0.67
32.66
±0.80
33.89
±0.67

40.33
±2.99
41.24
±4.01
39.94
±2.94
45.33
±2.88

32.68
±0.86
31.52
±0.46
34.15
±0.45
32.89
±0.36

34.13
±0.30
30.94
±0.53
34.47
±0.72
33.52
±0.78

31.70
±0.46
32.07
±0.36
33.93
±0.37
32.01
±0.59

32.21
±0.72
31.63
±1.08
34.62
±0.56
34.13
±0.72

27.34
±0.81
34.46
±0.86
31.85
±0.28
32.57
±0.40

27.11
±0.69
28.46
±0.54
30.57
±0.59
28.60
±0.36

Supplementary Numerical Results for Section 5.1

In Section 5.1 we summarize the results of the synthetic experiments, and here in Table 3, we
supplement their detailed numerical results. Overall, Accous outperforms almost all other methods
in the four shift patterns. Compared to the competitors that use only one round of online data, ourproposed method achieves a lower average error. The general domain adaptation method DANN
is also inferior in the online setting due to the limited amount of data per round, which makes it
difficult to learn a good representation. The proposed Accous algorithm also outperforms the OLRE
algorithm, which reuses all historical data. The above evidence demonstrates the need for selective
reuse of historical data and validates the effectiveness of our algorithm.

B

More Related Work

B.1

Related Work on Non-stationary Online Learning

Online convex optimization (OCO) [13] is a powerful paradigm to handle sequential prediction
problems. To adapt to the continuously shifting environments, we exploit the online ensemble
technique [15, 16] developed in non-stationary online learning literature to hedge the uncertainty. In
this part, we will introduce the related works and discuss the difference between our algorithm andthe previous work.
Specifically, OCO paradigm considers a T -round game between a learner and the environment. At
every round t, the learner makes a prediction θt ∈ Θ, and, at the same time, the environments reveal
the convex loss function Lt : Θ → R. Then, the learner suffers from the loss Lt (θt ) and observes
certain information about Lt to update the model for the next iteration. A classical performancemeasure for OCO framework in non-stationary environments is the dynamic regret,
e ∗ T
Regd
T ({Lt , θt }t=1 ) =

T
X
t=1

Lt (θt ) −TX

Lt (θt∗ ),

(10)

t=1

which compares the learner’s prediction with the function minimizer θt∗ at every iteration. Over the
decades, a variety of online learning algorithms have been proposed to optimize this dynamic regret
measure [38, 39, 55, 56, 57, 58, 59, 60, 61, 62, 63] in the online convex optimization problem.
Worst-case Dynamic Regret. Many works focus on the full information setting [57, 58, 59, 60, 63],
where the learner can observe the loss functions entirely. We call the dynamic regret in such a caseas “worse-case” dynamic regret in the sense that the comparator θt∗ is exactly the minimizer of theobserved function Lt . According to different complexity measures for quantifying the environmental
non-stationarity, there are two lines of results. The first is the path length of comparators VTθ =PT∗∗
t=2 ∥θt − θt−1 ∥2 measuring the variation of the function minimizers; and the second one is thePTvariation of the function values VTL = t=2 maxθ∈Θ |Lt (θ) − Lt−1 (θ)|. For the full information
setting, the prior art [63] showed that a simple greedy strategy can already be already well-behaved andachieves an O(min{VTθ , VTL }) dynamic regret bound that minimizes the two kinds of non-stationarity
measure simultaneously.
Universal Dynamic Regret. Competing to the sequences of observed functions’ minimizers can besometimes too pessimistic for the online optimizing, which may lead to a severe overfitting. A more
appropriate performance measure is the universal dynamic regret [38], which has drawn more andmore attention in recent years. Universal dynamic regret compares the learner’s prediction with an17

arbitrary comparator sequence {νt }Tt=1 ,TeRegd
T ({Lt , νt }t=1 ) =

T
X
t=1

Lt (θt ) −TX

Lt (νt ).

(11)

t=1

PT
Let VTν ={νt }Tt=1 . [38]t=1 ∥νt − νt−1 ∥2 be the path length of the comparator sequence√νshowed that the online gradient descent algorithm achieves an O((1 + VpT ) T ) universal dynamic
regret bound for convex function, whereas there is still a gap to the Ω( T (1 + VTν )) lower boundpestablished by [15], who further proposed an algorithm attaining an O( T (1 + VTν )) dynamic regretbound and thus closed the gap. Their key algorithmic ingredient is the online ensemble structure [64],which hedges the non-stationarity of the environments by employing a meta-learner to ensemble
the predictions
from a group of base-learners.
pp When, the loss function is convex and smooth,
the O( T (1 + VTν )) can be improved to O( (1 + min{FT , GT } + VTν )(1 + VTν )) [16], where
PT
PT
FT = t=1 Lt (νt ) and GT = t=2 maxθ∈Θ ∥∇Lt (θ) − ∇Lt−1 (θ)∥22 are two problem-dependent
quantities capturing different kinds of easiness of the environments. For the exp-concave and smoothfunctions, [19] showed an O(T 1/3 (VTν )2/3 ) universal dynamic regret by an improper algorithm (i.e.,the output decisions can be slightly out of the feasible domain, while comparators have to locate in thedomain), which is proven to be minimax optimal for exp-concave functions. Note that the worst-case
dynamic regret (10) is clearly a special case of the universal dynamic regret defined in (11), so auniversal dynamic regret upper bound directly implies an upper bound for the worst-case dynamic
regret, by substituting comparators {νt }Tt=1 as functions’ minimizer {θt∗ }Tt=1 . This method yields apath-length type bound. Furthermore, we mention that it is also possible to obtain a function-variation
type worst-case dynamic regret bound from the universal dynamic regret guarantee, which is revealedin [61, Appendix A.2].
An Intermediate Dynamic Regret. There is also an intermediate case between the worst-case anduniversal dynamic regret. We refer to it the “intermediate” dynamic regret in this paper. In thisintermediate case, the online learner still aims to optimize the worst-case dynamic regret of theb t (θ) at each iteration. Since theloss function Lt (θ), but she can only observe a noisy feedback L∗
b t (θ), such a kind ofcomparator θt = arg minθ∈Θ Lt (θ) is not the minimizer of the observed loss L
problem is more challenging than the worst-case dynamic regret minimization problem given theuncertainty of the comparator sequence. [39] showed a restarted online gradient descentq algorithm can
1/3
L 2/3
achieve an O(T (VT ) ) dynamic regret bound for convex function and an O( T VTL ) bound forstrongly convex function under the noisy feedback. However, their algorithm requires the knowledgePTof the non-stationarity measure VTL = t=2 maxθ∈Θ |Lt (θ) − Lt−1 (θ)| defined over Lt , whichis generally unknown. [48] showed that an O(T 1/3 (VTθ )2/3 ) bound is achievable with the noisy
feedback when Lt is a one-dimensional squared loss. The proposed method is based on the de-noisingtechnique and are free from the knowledge of VTθ . Evidently, this intermediate dynamic regret is stilla feasible realization of the universal dynamic regret, so one can achieve the intermediate dynamic
regret bound by the universal dynamic regret. However, this reduction will only provide an expectedguarantee. More detailed discussions are provided below.
1/3
e ∗ T
Comparison of Our Results and Earlier Works. Our Regd
(VTθ )2/3 )T ({Lt , θt }t=1 ) ≤ O(Thigh-probability dynamic regret bound for logistic regression (Theorem 2) falls into the intermediate
case, where we aim to minimize the worst-case dynamic regret bound in terms of Lt with noisy
feedback. The closest related work is by [19], who also exploited an online ensemble structure andT
1/3eachieved a universal dynamic regret in the form of Regd
(VTν )2/3 ) universal
T ({Lt , νt }t=1 ) ≤ O(T∗bdynamic regret by taking Lt as the input. Further choosing νt = θt and taking the expectation over
1/3
St for each round, their result implies an expectation bound E[Regd(VTθ )2/3 ). One mainT ] ≤ O(Tadvantage of our result is that the bound holds with high probability. To achieve the high-probability
bound, we twist the meta-learner from Hedge [42] to Adapt-ML-Prod [43], which helps controlb t and L
e t with the negative term introduced by the exp-concavity
the generalization error between Lof the loss function. Besides, similar to [19], our analysis crucially relies on providing a squared
formulation of the dynamic regret O(max{1, |I|(VTθ )2 }) on each interval I. However, to obtain sucha result, [19] crucially rely on a complicated analysis on the KKT condition to depict the structure ofcomparators. We greatly simplify the analysis by exploiting the fact that in our case the comparator
θt∗ is actually the minimizer of expected functions Lt , hence avoiding analyzing KKT condition. Our

18

analysis is applicable to the case when the minimizers lie in the interior of the feasible domain andrequire the smoothness of the online functions only (which does not rely on the particular propertiesof logistic loss). This result might be of independent interest in other online learning studies.

C

Omitted Details for Section 3

This section provides the proofs for Section 3 in the main paper.
C.1

Proof of Proposition 1b t at each time t and
Proof of Proposition 1. We first provide an excess risk bound for the model wethen combine them for T rounds. For each time t ∈ [T ], let Rt (w) = E(x,y)∼S0 rt∗ (x)ℓ(w⊤ x, y)bt (w) = E(x,y)∼S rbt (x)ℓ(w⊤ x, y)be the risk built upon the true density ratio function rt∗ and R0
be that built upon the empirical density ratio function rbt . We can decompose the excess risk of theb t astrained model w
et (w
et (w∗ ) + R
et (wt∗ ) − Rt (w∗ ) .
et (wb t) − Rb t ) − Rt (wt∗ ) = Rt (wb t) − Rb t) + R
Rt (w
t
t
{z
} |
{z
} |
{z
}
|
term (a)term (b)term (c)For term (a), by a standard generalization error analysis as shown in Lemma 2, with probability atet (wb t ) and Rb t ) is bounded byleast 1 − δ/(2T ), the gap between Rt (ws
2BGDR
2 ln((8T )/δ)et (wb t) − Rb t) ≤ √term (a) = Rt (w
.
+ 4BL
N0
N0
Then, we proceed to bound term (b).
et (wet (w∗ )b t) − Rterm (b) = R
t
et (w
bt (w
b t (w
bt (wt∗ ) + R
bt (wt∗ ) − Ret (wt∗ )b t) − Rb t) + Rb t) − R
=R
et (w
bt (w
bt (wt∗ ) − Ret (wt∗ )b t) − Rb t) + R
≤R
 ∗



b t⊤ x, y) + E(x,y)∼S0 (b= E(x,y)∼S0 (rt (x) − rbt (x))ℓ(wrt (x) − rt∗ (x))ℓ((wt∗ )⊤ x, y)≤ 2LEx∼S0 [|rt∗ (x) − rbt (x)|],b t and the second inequality is by the condition
where the first inequality is due to the optimality of w
maxx∈W,y∈Y,w∈W |ℓ(w⊤ x, y)| ≤ L.
Finally, as for term (c), since wt∗ is independent of S0 , a direct application of the Hoeffding’s
inequality [65, Lemma B.6] implies,set (w∗ ) − Rt (w∗ ) ≤ BL 2 ln((8T )/δ)
term (c) = Rtt
N0with probability at least 1 − δ/(2T ). Combining the upper bounds for term (a), term (b) and term (c)and taking the summation over T rounds, we gets
TX
2BGDRT
2 ln((8T )/δ)T√b t }t=1 ) ≤+ 5BLT+ 2L
Ex∼S0 [|rt∗ (x) − rbt (x)|]RT ({wN0
N0
t=1with probability at least 1 − δ, which completes the proof.
C.2

Proof of Proposition 2

Proof of Proposition 2. We begin with converting the estimation error to the squared formulation,
T
X
t=1

Ex∼D0 (x) [|rt∗ (x) − rbt (x)|] ≤

T q
X
t=1

19

Ex∼D0 (x) [|rt∗ (x) − rbt (x)|2 ]

v
u
u
≤ tT

T
X

!
Ex∼D0 (x) [|rt∗ (x) − rbt (x)|2 ]

,

(12)

t=1

where the first inequality is by Jensen’s inequality and the last inequality is due to the Cauchy-Schwarz
inequality. Then, we can bound the squared density ratio estimation error asµEx∼D0 (x) |b
rt (x) − rt∗ (x)|2 ≤ Ex∼D0 (x) [ψ (rt∗ (x)) − ψ (brt (x)) − ∂ψ (brt (x)) (rt∗ (x) − rbt (x))]
2∗= EBψ (rt∗ ∥brt ) = Lψrt ) − Lψ(13)t (bt (rt ),
where the first inequality is due to the strong convexity of the function ψ. The last equality is by thedefinition of the loss function Lψ
t . Combining (12) and (13), we have
v
!
u
T
T
T
u 2T X
X
X
ψ
ψ ∗
∗
t
Ex∼D0 (x) [|brt (x) − rt (x)|] ≤L (brt ) −Lt (rt ) .
µ t=1 t
t=1
t=1

C.3

Proof of Theorem 1

Proof of Theorem 1. Following the same arguments in obtaining (12), the empirical cumulative
estimation error of the density ratio estimators can be bounded as
v
!
u
T
T
h
i
u
X
X
(14)
E
|r∗ (x) − rb (x)|2 .
E[|r∗ (x) − rb (x)|] ≤ tT
x∼S0 (x)ttx∼S0 (x)

t

t

t=1

t=1

Then, we proceed to bound the squared cumulative estimation error by the regret of the loss functioneψLt (r) defined in (5). For notation simplicity, we define
f1 (z) = ∂ψ(z)z − ψ(z)andf2 (z) = ∂ψ(z).

eψIn such a case, the loss function Lt (r) = Ex∼S0 [f1 (r(x))] − Ex∼Dt (x) [f2 (r(x))] and the regretψ
e t can be written as
over L∗eψeψLrt ) − Lrt (x)) − f1 (rt∗ (x))] − Ex∼Dt [f2 (brt (x)) − f2 (rt∗ (x))].
t (bt (rt ) = Ex∼S0 [f1 (b(15)

For the first term of R.H.S. of (15), by the Taylor’s Theorem, for each x ∈ X , we have
Ex∼S0 [f1 (brt (x)) − f1 (rt∗ (x))]
∂ 2 f1 (ξx )∗∗∗2= Ex∼S0 ∂f1 (rt (x))(brt (x) − rt (x)) +(b
rt (x) − rt (x))
2


 µ


≥ Ex∼S0 ∂ 2 ψ rt∗ (x) rt∗ (x) rbt (x) − rt∗ (x) + Ex∼S0 (brt (x) − rt∗ (x))2 .
2(16)

In above ξx ∈ dom ψ is a certain point on the line collecting rbt (x) and rt∗ (x). The second inequality
is due to the definition of f1 such that,∂ 2 f1 (ξx ) = ∂ 2 ψ(ξx ) − ∂ 3 ψ(ξx )ξx ≥ µ,where the inequality holds due to the µ-strong convexity of ψ(z) and the condition ξx ∂ 3 ψ(ξx ) ≤ 0
for all ξx ∈ domψ.
Then, we proceed to bound the second term of R.H.S. of (15),Ex∼Dt [f2 (brt (x)) − f2 (rt∗ (x))]
 2




∂ f2 (ξx )2= Ex∼Dt ∂f2 rt∗ (x) (brt (x) − rt∗ (x)) + Ex∼Dt(brt (x) − rt∗ (x))
2
 2


 2


∗
∗
rt (x) − rt (x)) = Ex∼D0 ∂ ψ rt∗ (x) rt∗ (x) (brt (x) − rt∗ (x)) ,≤ Ex∼Dt ∂ ψ rt (x) (b
20

(17)

In above ξx ∈ dom ψ is a certain point on the line connecting rbt (x) and rt∗ (x). The first inequality is
due to the definition of f2 and the condition on ψ, such that ∂ 2 f2 (ξx ) = ∂ 3 ψ(ξx ) ≤ 0. Plugging (16)and (17) into (15), and rearranging the terms, we arrive


µ
∗
eψ
eψ
Ex∼S0 (brt (x) − rt∗ (x))2 ≤ Lrt ) − L(18)t (bt (rt ) + Ut ,2where the term Ut is defined as






Ut = Ex∼D0 ∂ 2 ψ rt∗ (x) rt∗ (x) (brt (x) − rt∗ (x)) − Ex∼S0 ∂ 2 ψ rt∗ (x) rt∗ (x) rbt (x) − rt∗ (x) .
It remains to analyze the generalization gap Ut for every iteration. The main challenge is that thelearned model rbt depends on the initial dataset S0 and the generalization error is related to thecomplexity of the hypothesis of the density ratio estimator. In the following lemma, we use thecovering number [65, Chapter 27] to measure the complexity of the hypothesis space. It is possible
to provide an upper bound of Ut with other advanced complexity measure.
Lemma 1. Let Hr ≜ {r : X → R | ∥r∥∞ ≤ Br′ } be a hypothesis space of the density ratio function
whose covering number is N (Hr , ε, ∥·∥∞ ) . Then, for all model r ∈ Hr and δ ∈ (0, 1), we have thefollowing upper bound






Ut = Ex∼D0 ∂ 2 ψ rt∗ (x) rt∗ (x) (r(x) − rt∗ (x)) − Ex∼S0 ∂ 2 ψ rt∗ (x) rt∗ (x) r(x) − rt∗ (x)Cr (µ) log (2N (Hr , ε, ∥·∥∞ )/δ) µµε2≤+ ES0 (r(x) − rt∗ (x))2 + 2Ar ε +,N044with probability at least 1 − δ for any ε > 0, where Cr (µ) = 4Ar Br /3 + 16A2r /µ + 3µBr2 =O(max{µ, 1/µ}) with constants Ar = |maxx∈X ∂ 2 ψ(rt∗ (x))rt∗ (x)|, and Br = max{B, Br′ }.
Combining (18) and Lemma 1 with the hypothesis space Hr = Hθ ≜ {x 7→ h(x, θ) | θ ∈ Θ} andrearranging the terms yields,µ
Cr (µ) log (2N (Hθ , ε, ∥·∥∞ )/δ)
µε2
∗
eψ
eψ
Ex∼S0 (brt (x) − rt∗ (x))2 ≤ Lrt ) − L.
+ 2Ar ε +t (bt (rt ) +4
N04
Further choosing ε = 1/T and taking the summation of the T iterations, we have
T



µX
Ex∼S0 (brt (x) − rt∗ (x))2
4 t=1
≤

T
X
t=1

eψ
L
rt ) −t (bTX

Cr (µ)T log (2T N (Hθ , /T , ∥·∥∞ )/δ)
∗
eψ
L
+ Dr
t (rt ) +
N0
t=1
1

(19)

where Dr = 2Ar + µ/(4T ) = O(1) is a constant. Then, we proceed to bound the covering numberof the hypothesis space Hθ . Let θ, θ′ ∈ Θ be the corresponding parameters for the two density ratio
function r, r′ ∈ Hθ . Then, we can show that for any ∥θ − θ′ ∥2 ≤ ε, we have∥r − r′ ∥∞ = max |h(x, θ) − h(x, θ′ )| ≤ Gh ∥θ − θ′ ∥2 ,x∈X

where Gh = maxx∈X ,θ∈Θ ∥∇h(x, θ)∥2 is the Lipschitz continuity constant. Then, we can check thatthe covering number of Hθ in terms of ∥·∥∞ can be bounded by that of Θ in terms of ∥·∥2 as
N (Hθ , 1/T , ∥·∥∞ ) ≤ N (Θ, 1/(Gh T ), ∥·∥2 ) ≤ (3SGh T )d .

(20)

In the above, the last inequality holds because the parameter space Θ is essentially a L2 -ball withradius S, whose covering number is bounded by (3S/ε)d [66]. Then, plugging (20) into (19), we
arriveTT
 4dT C (µ) log (6SG T /δ) 4DX
4 X  eψrhr∗eψEx∼S0 (brt (x) − rt∗ (x))2 ≤Lt (brt ) − L
(r
)
+
+
.
t
t
µ
µN
µ
0
t=1
t=1
Then, we complete the proof by plugging this inequality into (14),T
X

Ex∼S0 (x) [|brt (x) − rt∗ (x)|]

t=1

21

v
!
u
T
T
u 4T XX
4dT 2 Cr (µ) log (6SGh T /δ) 4T Dr
ψ
ψ ∗
t
Lt (b
rt ) −+≤Lt (rt ) +
µ t=1
µN0
µ
t=1
v
!
( T
)
u
√


T
u 4T
X ψ
X
T d log(T /δ)
1
ψ ∗
t
√
=
.
max
·
Lt (b
rt ) −Lt (rt ), 0 + O max 1,
µ
µ
N0
t=1
t=1

C.4

Useful Lemmaset (w) = E(x,y)∼S rt∗ (x)ℓ(w⊤ x, y) . With probability at least 1 − δ over theLemma 2. Let R
0drawn of S0 , all predictions w ∈ W satisfys
2BGDR
2 ln(4/δ)et (w)| ≤ √+ 4BL|Rt (w) − R.
N0
N0et (w) is bounded by BL for all w ∈ W, the standard
Proof of Lemma 2. Since the risk function R
analysis on generalization bound with Rademacher complexity (e.g. [65, Theorem 26.5]) shows thatsb S (Lt ) + 4BL 2 ln(4/δ)
et (w)| ≤ 2R
|Rt (w) − R0N0
for any w ∈ W. In above, the function spaceas Lt = {(x, y) 7→ rit∗ (x)ℓ(w⊤ x, y) |h Lt is definedPN
0b S (Lt ) = 1 Eσ∼{±1}N0 sup
w ∈ W} and R
σi rt (xn )ℓ(w⊤ xn , yn ) is the empirical
0

N0

w∈W

n=1

Rademacher complexity of Lt .
Note that the function ϕn : R → R defined as ϕn (z) = rt∗ (xn )ℓ(z, yn ) is BG-Lipschitz, such that|ϕn (w1⊤ xn ) − ϕn (w2⊤ xn )| ≤ BG|(w1 − w2 )⊤ xn | for any w1 , w2 ∈ W. Then, according to the
Talagrand’s lemma [65, Lemma 26.9], we haveb S (Lt ) ≤ BGR
b S (W),fR00
f = {x 7→ w⊤ x | w ∈ W} is a hypothesis space for the linear functions. By [65, Lemma
where W
26.10], the Rademacher complexity of the linear function class is bounded byDRb S (W)f ≤ D maxx√n ∈S0 ∥xn ∥2 ≤ √
R
.
0
N0
N0
We complete the proof by combining the above displayed inequalities.
Proof of Lemma 1. We begin with the analysis for the generalization gap for a fixed model r ∈ Hr .
For notation simplicity, we denote byfr (x) = ∂ 2 ψ rt∗ (x) rt∗ (x) (r(x) − rt∗ (x)) ,and let the random variable Zi = Ex∼D0 [fr (x)] − fr (xi ) for any xi ∈ S0 . Since xi is i.i.d. sampledfrom D0 and fr is independent of S0 , we have E[Zi ] = 0 and |Zi | ≤ 2Ar Br with the constants
Ar = maxx∈X |∂ 2 ψ(rt∗ (x))rt∗ (x)| and Br = max{B, maxx∈X ,r∈Hr |r(x)|}. Besides, the varianceof Zi is bounded byE[Zi2 ] = Ex∼D0 [(fr (x))2 ] − (Ex∼D0 [fr (x)])2 ≤ Ex∼D0 (fr (x))2 ,Then, by the Bernstein’s inequality [65, Lemma B.9], with probability at least 1 − δ, we haves
4Ar Br log(1/δ)
2 log(1/δ)ED0 [(fr (x))2 ])Ex∼D0 [fr (x)] − Ex∼S0 [fr (x)] ≤
+
3N0
N0


2
4Ar Br
16Ar log(1/δ)
µ
2
≤
+
+
ED [(fr (x)) ]3µN0
8A2r 022

≤

4Ar Br
16A2r
+
3
µ





log(1/δ) µ+ ED0 (r(x) − rt∗ (x))2 ,
N0
8
(21)

where the first inequality is due to the AM-GM inequality. The second inequality is due to the factthat |fr (x)| ≤ |∂ 2 ψ(rt∗ (x))rt∗ (x)| · |r(x) − rt∗ (x)| ≤ Ar |r(x) − rt∗ (x)| for any x ∈ X .


∗
2
Then,
 since r∗ is a2 fixed model in Hr , we can further bound ED0 (r(x) − rt (x)) by
ES0 (r(x) − rt (x)) by applying the Bernstein’s inequality to the random variable
Yi =

∗
2
2
∗
2
2
(r(xi ) − rt (xi )) /(4Br ). Specifically, after checking E[Yi ] = ED0 (r(x) − rt (x)) /(4Br ) andYi ∈ [0, 1], a direct application of [65, Lemma B.10] implies
ED0 (r(x) − rt∗ (x))2
s
2


16B
log(1/δ)8Br2 ES0 [(r(x) − rt∗ (x))2 ] log(1/δ)
r≤ ES0 (r(x) − rt∗ (x))2 +
+
N0
N0
≤



24Br2 log(1/δ)+ ES0 (r(x) − rt∗ (x))2 ,
N0

(22)

where the last inequality is due to the AM-GM inequality. Plugging (22) into (21), for a fixed r ∈ Hr ,we haveEx∼D0 [fr (x)] − Ex∼S0 [fr (x)]




16A2r
log(2/δ) µ
4Ar Br2≤++ 3µBr+ ES0 (r(x) − rt∗ (x))2 ,3µN08with probability at least 1 − δ.
Then, based on the notion of covering number, we show a union bound of Ut over all r ∈ Hr . Let
N (Hr , ε, ∥·∥∞ ) be the ε-net of Hr such that for any r ∈ Hr one can find a r′ ∈ N (Hr , ε, ∥·∥∞ )satisfying ∥r − r′ ∥∞ ≤ ε. The covering number N (Hr , ε, ∥·∥∞ ) is defined as the minimal cardinalityof an ε-net of Hr . By taking the union bound over all r′ ∈ N (Hr , ε, ∥·∥∞ ), with probability at least
1 − δ, the following holds for any r′ ∈ N (Hr , ε, ∥·∥∞ ),Ex∼D0 [fr′ (x)] − Ex∼S0 [fr′ (x)]∞)Cr (µ) log 2N (Hr ,ε,∥·∥


δ
µ
≤
+ ES0 (r′ (x) − rt∗ (x))2 ,
N0
8

(23)

where Cr (µ) = 4Ar Br /3 + 16A2r /µ + 3µBr2 = O(max{µ, 1/µ}). Then, for any model r ∈ Hr ,we can show
Ex∼D0 [fr (x)] − Ex∼S0 [fr (x)]= Ex∼D0 [fr′ (x)] − Ex∼S0 [fr′ (x)] + Ex∼D0 [fr (x) − fr′ (x)] − Ex∼S0 [fr (x) − fr′ (x)]≤ Ex∼D0 [fr′ (x)] − Ex∼S0 [fr′ (x)] + 2εAr(23) Cr (µ) log (2N (Hr , ε, ∥·∥∞ )/δ)µ≤
+ ES0 (r′ (x) − rt∗ (x))2 + 2Ar εN08Cr (µ) log (2N (Hr , ε, ∥·∥∞ )/δ) µµε2≤+ ES0 (r(x) − rt∗ (x))2 + 2Ar ε +,N044where the first inequality is due to the property of the ε-net such thatEx∼D0 [fr (x) − fr′ (x)] ≤ Ex∼D0 [|∂ 2 ψ(rt∗ (x))r′ (x)| · |r(x) − rt∗ (x)|] ≤ Ar ∥r − r′ ∥∞ ≤ Ar ε.
The last inequality is by

 µ

 µ


µ
ES0 (r′ (x) − rt∗ (x))2 ≤ ES0 (r(x) − rt∗ (x))2 + ES0 (r(x) − r′ (x))2
8
4
4

 µ

 µε2
µ
µ
∗
2
′ 2
≤ ES0 (r(x) − rt (x)) + ∥r − r ∥∞ ≤ ES0 (r(x) − rt∗ (x))2 +
.
4
4
4
4
Thus, we complete the proof.
23

Algorithm 1 Base-learner Ei
Input: Active interval Ii = [si , ei ], regularization parameter λ and step size η, feasible domain Θ.
1: Initialize the model θbsi ,i as any point in the domain Θ and the matrix Asi −1 = λId .
2: for t ∈ Ii = [si , ei ] do
3:
Submit the model θbt,i to the meta-learner.
b t (θbt,i )∇Lb t (θbt,i )⊤ .
4:
Update the matrix by At,i = At−1,i + ∇L
5:Update the model for next iteration by

A 
b b
θbt+1,i = ΠΘt,i θbt,i − γA−1t,i ∇Lt (θt,i ) ,Awhere ΠΘt,i [ν] = arg minθ∈Θ ∥θ − ν∥At,i is the projection operator.
6: end forD

Omitted Details for Section 4This section presents the omitted details for Section 4. We will first introduce the algorithm details
for our online ensemble algorithm, characterized by its meta-base structure, in Appendix D.1. Then,
we present the theoretical guarantees alongside their proofs for the base-learner (Lemma 3) and metalearner (Lemma 10) in Appendix D.2 and Appendix D.3 respectively. These elements collectively
contribute to the dynamic regret guarantee of the overall algorithm. The dynamic regret of the online
ensemble algorithm (Theorem 2) is presented in Appendix D.4. We provide the proof for the average
excess risk (Theorem 3) in Appendix D.5. Further discussions regarding the tightness of the boundand alternative choices for the Bregman divergence function are located in Appendices D.6 and D.7.
D.1

Omitted Algorithm Details

Here, we introduce the algorithmic details for the online ensemble method omitted in Section 4.1.
Interval Schedules. We run multiple based learners on the geometric covering:[C=
Ck and Ck = [i · 2k , (i + 1) · 2k − 1] | i ∈ N and i · 2k ≤ T
k∈N∪{0}(24)

and then employ a meta-learner to combine them by weights. For any interval Ii = [si , ei ] ∈ C,
which starts at time t = si and end at t = ei , we denote by Ei the base-learner running over Ii .
Base-learner. The base-learner Ei updates his/her model with the ONS algorithm [44] in theassociated active period Ii . The algorithmic details of the base-learner are summarized in Algorithm 1.
At every iteration t ∈ Ii , the base-learner Ei will first submit his/her model θbt,i to the meta-learnerand then update the model for the next iteration by Line 4 and Line 5 as shown in Algorithm 1.
Meta-learner. The meta-learner aggregates the predictions of each active base-learner by a weighted
combination scheme. To obtain the weight pt,i for each active base-learner Ei , we maintain threeterms mt,i , vt,i and εt,i for each base-learner Ei when he/she is active, i.e., t ∈ Ii . Let gt (θ) =b t (θbt ), θ − θbt ⟩ be a linearized loss for the original function L
b t (θ). The first term is defined by⟨∇L
mt,i ≜b t (θbt ), θbt − θbt,i ⟩⟨∇L
gt (θbt ) − gt (θbt,i )
=
,
SR
SR

which measures the performance gap between the meta-learner’s prediction θbt and the i-th baselearner’s prediction θbt,i over the linearized loss gt (θ). The second term vt,i is the “potential” of thebase-learner Ei , which measures the historical performance of Ei over the interval. When invoking
the base-learner Ei at t = si , we initialize vsi −1,i = 1/K and the term is updated byεt,j εt−1,jvt,i = vt−1,i · 1 + εt−1,i mt,i, ∀t ∈ Ii .
(25)
The last term εt,i is the learning rate, which is set as
( s
)
1
ln K
,
.
εt,i = minPt2
1 + τ =si m2τ,i
24

(26)

Algorithm 2 Meta-learner
Input: Interval set for base-learners C.
1: Iitialize the set for active base-learner A0 = ∅ and let K = |C|.
2: for t = 1, . . . , T do
3:
for Ii = [si , ei ] ∈ C do
4:
if si is equal to t then5:
Create a base-leaner Ei by Algorithm
1 and initialize the potential vsi −1,i = 1/K and√step size εsi −1,i = min{1/2, ln K}.
6:
Add the index to the active set At = At−1 ∪ {i}.
7:else if ei + 1 is equal to t then8:
Remove the index from the active set At = At−1 /{i}.
9:end if
10:
end for
11:
Receive the θbi,t from active base-learners Ei ∈ At and update their weights pt,i by (27).
P
12:Submit the model by θbt = i∈At pt,i θbt,i .
13:
for i ∈ At do
14:
Update the step size εt+1,i for all current active base-learner by (26).
15:
Update the potential vt+1,i for all current active base-learner by (25).
16:
end for17: end for

After obtaining the vt,i and εt,i for each active base-learner, the weight for the next iteration t + 1 isupdated by
εt,i vt,i, for all i ∈ At+1 .
i∈At+1 εt,i vt,i

pt+1,i = P(27)

The algorithm is summarized in Algorithm 2. At the beginning of each iteration t, the algorithm will
update the index set for the active base-learner (Line 3-Line 10). Then, the meta-learner’s prediction
is generated by the weighted combination rule (Line 11). After that, the meta-learner will update thepotential vt+1,i and the step εt+1,i for the next iteration.
D.2

Regret Guarantee for Base-Algorithm (Lemma 3)

By suitable configurations of the parameters γ, λ > 0, we show that the ONS algorithm enjoys thefollowing dynamic regret guarantees,
Lemma 3. Let θt∗ = arg minθ∈Rd Lt (θ) and all the minimizer θt∗ ’s lie in the interior of Θ. Set
γ = 6(1 + β) and λ = 1. Then, with probability at least 1 − δ, the ONS running on the interval
Ii = [si , ei ] ⊆ [T ] ensures


X
X


e t (θbt,i ) −e t (θ∗ ) ≤ O dβ log |Ii |/δ + |Ii | V θ 2 + |Ii | log(dT /δ)
L
L
t
Ii
N0
t∈Ii

where VIθi =
D.2.1

t∈Ii

Pei

∗
∗
t=si +1 ∥θt − θt−1 ∥2 is the path length measuring the fluctuation of the minimizers.

Main Proof

This section provides the proof of the ONS algorithm. Since the analysis of the single base-learner is
independent from the online ensemble structure, for notation simplicity, we omit the subscribe of
the base-learner in the analysis. Specifically, we abbreviate the prediction θbt,i , matrix At,i and theinterval Ii as θbt , matrix At and the interval I respectively.
Proof of Lemma 3. We begin the proof of the dynamic regret with a lemma on the static regret of
the ONS algorithm, whose proof is proved in Appendix D.2.2. For notation simplicity, we omit thesubscribes of the base-learner, where the prediction θt,i , matrix At,i and interval Ii is abbreviated as
θt , At and I, respectively.
25

Lemma 4. Set γ = 6(1 + β) and λ = 1. With probability at least 1 − 2δ, for any δ ∈ (0, 1], theONS algorithm (8) running on the interval I = [s, e] ⊆ [T ] ensuresXX
e t (θbt ) −e t (θ) ≤ C (1) (δ) = O (dβ log(|I|/δ))
L
L
I
t∈I

t∈I

√
(1)
for any θ ∈ Θ, where CI (δ) = 12dβ ′ log

1 + |I|R2

S2
6S 2 R2′+ 3β′ + (48β +β ′ ) log

2|I|+1+δq
2 2and β ′ = β + 1.
(4SR + S βR′ ) log 2|I|+1
δ2
Then, we can decompose the dynamic regret as follows,
X
X
X
X
X
X
e t (θbt ) −e t (θ∗ ) =e t (θbt ) −e t (θ̄∗ ) +e t (θ̄∗ ) −e t (θ∗ )
L
L
L
L
L
L
t

t∈I

I

t∈I

t∈I

I

t∈I

(1)

≤ CI (δ) +Xt∈Ie t (θ̄I∗ ) −
L

t∈I

X

t

t∈I

e t (θt∗ ),
L

(28)

t∈I

∗
t∈I θt is an averaged prediction of the function minimizers over T rounds. Under
the condition all function minimizers lie in the interior of Θ, we have θ̄I∗ ∈ Θ. As a consequence, the1where θ̄I∗ = |I|P

last inequality holds due to Lemma 4.
Then, we proceed to bound the last two terms of the R.H.S. of the above inequality. Since thefirst-order derivative of the sigmoid function σ ′ (z) = exp(−z)/(1 + exp(−z))2 ≤ 1/4 for any
e t is a R2 /4-smooth function such thatz ∈ R. Let Id ∈ Rd×d be the identity matrix. we can check L
e t (θ) = 1 Ex∼S [σ ′ (x⊤ θ)ϕ(x)ϕ(x)⊤ ] + 1 Ex∼D [σ ′ (x⊤ θ)ϕ(x)ϕ(x)⊤ ] ≼ R2 Id for any θ ∈ Θ.
∇2 L
0
t
2
2
4
Then, with probability at least 1 − 2δ, we have2e t (θ̄I∗ ) − L
e t (θt∗ ) ≤ ⟨∇Le t (θt∗ ), θ̄I∗ − θt∗ ⟩ + R ∥θ̄I∗ − θt∗ ∥22
L
8
2
R
8
∗
2
e t (θt )∥2 +≤ 2 ∥∇L∥θ̄∗ − θt∗ ∥22
R
4 I


16 ln((d + 1)/δ) R2 ∗log(d/δ)∗ 2∥θ̄ − θt ∥2 + o
≤
+
N0
4 I
N0


2
16 ln((d + 1)/δ) R
log(d/δ)θ 2(V ) + o
≤
+
(29)
N0
4 I
N0
Pe
∗
where VIθ = t=s+1 ∥θt−1− θt∗ ∥2 with s being the start time of the interval I and e being the ending
time. In above, the first inequality is due to the smoothness of the loss function and the second
inequality comes from the AM-GM inequality. The third inequality is due to Lemma 5. The lastinequality can be obtained by the definition of θ̄I∗ such that∥θ̄I∗ − θt∗ ∥22 =

!2

1 X ∗
θτ − θt∗
|I|
τ ∈I

≤

1 X ∗
∥θτ − θt∗ ∥2|I|!2≤ (VIθ )2 .

τ ∈I

2

Plugging (29) into (28) and taking the summation over T rounds, we obtainXX
e t (θbt ) −e t (θ)
L
L
t∈I

t∈I

16|I| ln((d + 1)|I|/δ) R2(1)≤ CI (δ) ++|I|(VIθ )2 + oN04|I| ln(d|I|/δ)θ 2= O dβ log(|I|/δ) ++ |I|(VI )N0with probability at least 1 − 3δ. We have completed the proof.
26

|I| log(d|I|/δ)N0

D.2.2

Useful Lemmas

Lemma 4. Set γ = 6(1 + β) and λ = 1. With probability at least 1 − 2δ, for any δ ∈ (0, 1], theONS algorithm (8) running on the interval I ⊆ [T ] ensures,XX
e t (θbt ) −e t (θ) ≤ C (1) (δ) = O (dβ log(|I|/δ))
L
L
I
t∈I

t∈I

√

2|I|+1(1)S2
6S 2 R2′
for any θ ∈ Θ, where CI (δ) = 12dβ ′ log 1 + |I|R2 + 3β)log+
′ + (48β +
′
β
δ
q
2 2
2|I|+1(4SR + S βR′ ) log δ2 and β ′ = β + 1.
Proof of Lemma 4. The main difference between the proof of Lemma 4 and the standard analysis for
the ONS algorithm (e.g., proof of [13, Theorem 4.3]) is that the later one is for the full informatione t . However, in
online learning setting, in which the learner can exactly observe the loss function Lbthe continuous covariate shift problem, we can only observe the empirical loss Lt , which poses thee t and the empirical observation
challenge to bound a generalization gap between the expected loss L
b t (e.g. term (b) in (35)). We manage to show such a gap can be bounded by O(log T ) with theLself-normalized concentration inequality for the martingale different sequence.
We start with the analysis of the instantaneous regret. By the Taylor’s theorem, we haveb t (ξt )(θbt − θ),b t (θbt ) − L
b t (θ) = ⟨∇Lb t (θbt ), θbt − θ⟩ − 1 (θbt − θ)⊤ ∇2 L
(30)
L
2where ξt ∈ Rd is a certain point on the line connecting θ and θbt . In above, the second order derivative
e t can be further lower bounded by∇2 L
1
1Ex∼S0 [σ ′ (x⊤ ξt )ϕ(x)ϕ(x)⊤ ] + Ex∼St [σ ′ (x⊤ ξt )ϕ(x)ϕ(x)⊤ ]22

1
Ex∼S0 [ϕ(x)ϕ(x)⊤ ] + Ex∼St [ϕ(x)ϕ(x)⊤ ]≽SR
4(1 + e )1b t (θbt )∇Lb t (θbt )⊤ ,≽∇L
2(1 + β)b t (ξt ) =
∇2 L

(31)

where σ ′ (z) = exp(−z)/(1 + exp(−z))2 is the first-order derivative of the sigmoid function
σ(z) = 1/(1 + exp(−z)) and A ≽ B indicates A − B is a positive semi-definite matrix for any
e t and the second inequality is due
matrix A, B ∈ Rd×d . The first equality is due to the definition of L
to Lemma 6.
Plugging (31) into (30), we haveb t (θbt ) − L
b t (θ) ≤ ⟨∇Lb t (θbt ), θbt − θ⟩ −
L


2
1
b t (θbt )⊤ (θbt − θ) .
∇L
4(1 + β)(32)

Let Et [·] = ESt ∼Dt [· | S0 , S1 , . . . , St−1 ] be the expectation taken over the draw of St conditioned onb t (θbt )] = Le t (θbt ) andthe randomness until round t − 1. Since θbt is independent of St , we have Et [L
b t (θbt )] = ∇Le t (θbt ). Taking the expectation over both sides of (32) yields
Et [∇Le t (θbt ) − L
e t (θ)Le t (θbt ), θbt − θ⟩ −
≤ ⟨∇L

h
 i
1
b t (θbt )⊤ (θbt − θ) 2
Et ∇L
4(1 + β)b t (θbt ), θbt − θ⟩ + ⟨∇L
e t (θbt ) − ∇Lb t (θbt ), θbt − θ⟩ −
= ⟨∇L

h
 i
1
b t (θbt )⊤ (θbt − θ) 2
Et ∇L
4(1 + β)(33)

where the first term of the R.H.S. can be further bounded byb t (θbt ), θbt − θ⟩⟨∇Lb t (θbt ), θbt+1 − θ⟩ + ⟨∇L
b t (θbt ), θbt − θbt+1 ⟩
= ⟨∇L
27


1  b
b t (θbt ), θbt − θbt+1 ⟩∥θt − θ∥2At − ∥θbt+1 − θ∥2At − ∥θbt − θbt+1 ∥2At + ⟨∇L
2γ

1  b
b t (θbt )∥2 −1 + 1 ∥θbt − θbt+1 ∥2A≤∥θt − θ∥2At − ∥θbt+1 − θ∥2At − ∥θbt − θbt+1 ∥2At + 2γ∥∇L
t
At
2γ
2γ


b t (θbt )∥2 −1 + 1 ∥θbt − θ∥2 − ∥θbt+1 − θ∥2= 2γ∥∇L
At
At
At
2γ


b t (θbt )∥2 −1 + 1 ∥θbt − θ∥2A − ∥θbt+1 − θ∥2A + 1 (∇Lb t (θbt )⊤ (θbt − θ))2 .
= 2γ∥∇L
(34)
t−1
t
At
2γ
2γ
In above, the first inequality is due to the update rule of the ONS algorithm and the Bregman proximal
inequality [67, Lemma 5]. The second inequality is due to the AM-GM inequality. The last equalityb t (θbt )⊤ (θbt − θ))2 .
is due the definition of At such that ∥θbt − θ∥2At − ∥θbt − θ∥2At−1 = (∇L≤

Plugging (34) into (33) and taking the summation over the interval I = [s, e], we haveXX
e t (θbt ) −e t (θ)
L
L
t∈I

t∈I

X
b t (θbt )∥2 −1 + 1 ∥θbs − θ∥2 + 1
b t (θbt )⊤ (θbt − θ))2
∥∇L
(∇L
As
At
2γ
2γ
t∈I
t∈I
|
{z
}

≤ 2γ

X

term (a)+Xe t (θbt ) − ∇Lb t (θbt ), θbt − θ⟩ −
⟨∇L

(35)

X h
 i
1
b t (θbt )⊤ (θbt − θ) 2 .
Et ∇L
4(1 + β)
t∈I

t∈I

{z

|

}

term (b)Now, we proceed to bound term (a) and term (b), respectively. Firstly, according to Lemma 7, we canshow term (a) is bounded by


|I|R2
term (a) ≤ 2dγ log 1 +.
λThen, we handle term (b) with the Bernstein-type self-normalized concentration inequality [68,Theorem 4] . Specifically, let Et [·] = ESt ∼Dt [· | S0 , S1 , . . . , St−1 ] be the expectation taken overe t (θbt ) −the draw of St conditioned on the randomness until round t − 1. Denoting by Zt = ⟨∇Lb t (θbt ), θbt −θ⟩, it is easy to check {Zt }et=s is a martingale difference sequence such that Et [Zt ] = 0∇Land |Zt | ≤ 4SR. Lemma 8 with the choice ν = 1/(48(1 + β)) indicates,

2 
X
2
1
⊤ b
⊤ b
b
b
b
b
∇Lt (θt ) (θt − θ) + Et ∇Lt (θt ) (θt − θ)term (b) ≤24(1 + β)t∈Irp
2|I| + 1
2|I| + 1+ 48(1 + β) log+ 4SR log
.
(36)
δ
δ2
Combining the upper bounds for term (a) and term (b) with (35), we haveXX
e t (θbt ) −e t (θ)
L
L
t∈I

t∈I

r
p


2|I| + 1
|I|R2
2λS 2
2|I| + 1≤ 2dγ log 1 +++ 48(1 + β) log+ 4SR log
λ
γ
δ
δ2

X
h
X
 i
1
1
5
b t (θbt )⊤ (θbt − θ))2 −b t (θbt )⊤ (θbt − θ) 2
+
+
(∇L
Et ∇L
2γ
24(1 + β)24(1 + β)
t∈I
t∈I
|
{z
}
term (c)(37)
Then, setting γ = 6(1 + β), with probability at least 1 − δ, we can bound term (c) by!
2 i
1
3X b b ⊤ b
5X h b b ⊤ b
2term (c) ≤∇Lt (θt ) (θt − θ)) −Et ∇Lt (θt ) (θt − θ)
6(1 + β) 4
4
t∈I

t∈I

28

6S 2 R2
≤
log
1+β

p

2|I| + 1 S 2 R2
+
δ
1+β

r
log

2|I| + 1,
δ2

where the last inequality is due to Lemma 9. Combining the upper bound for term (c) and (37) andsetting λ = 1, we haveXX
e t (θbt ) −e t (θ)
L
L
t∈I

t∈I

S2
3β + 3
p



r
2 2
2|I| + 1
6S R
2|I| + 1
S 2 R2+ 48(β + 1) +log
log
+ 4SR +
1+β
δ
β+1
δ2
=O(dβ log(|I|/δ)),≤ 12d(1 + β) log 1 + |I|R2 +with probability at least 1 − 2δ, which completes the proof.
Lemma 5. When all the minimizers θt∗ = arg minθ∈Rd Lt (θ) lie in the interior of Θ, then we havelog(d/δ)2R2 ln((d + 1)/δ)
∗ 2
e
+o
.
∥∇Lt (θt )∥2 =N0
N0
Proof of Lemma 5. Since all the minimizers θt∗ = arg minθ∈Rd Lt (θ) lie in the interior of Θ, thene t (θt∗ )∥2 can be rewritten as
we have ∇Lt (θt∗ ) = 0 and the term ∥∇L2e t (θt∗ )∥22 = ∥∇Le t (θt∗ ) − ∇Lt (θt∗ )∥22∥∇L
1
= ∥Ex∼S0 [(σ(x⊤ θt∗ ) − 1)ϕ(x)] − Ex∼D0 [(σ(x⊤ θt∗ ) − 1)ϕ(x)]∥22
4
2
N0
1 1 X
=
Zn ,
4 N0 n=1

(38)

2

⊤ ∗
⊤ ∗
0
where {Zn }N
n=1 with Zn = (σ(ϕ(xn ) θt ) − 1)ϕ(xn ) − Ex∼D0 [(σ(ϕ(x) θt ) − 1)ϕ(x)] are N0
independent random vectors. We can check that Exn ∼D0 [Zn ] = 0 and ∥Zn ∥2 ≤ 2R for any xn ∈ X .
Then by the Hoeffding’s inequality for the random vector [69, Theorem 6.1.1], with probability atleast 1 − δ, we haves
N0
2R ln((d + 1)δ)1 X
8R2 ln((d + 1)/δ)
Zn ≤
+
.
(39)
N0 n=1
3N0
N0
2


2
e t (θt∗ )∥2 = 2R ln((d+1)/δ) + o log(d/δ) .
Plugging (39) into (38) obtains ∥∇L2
N0
N0

b t (θbt ) = 1 Ex∼S [(σ(ϕ(x)⊤ θbt ) − 1)ϕ(x)] + 1 Ex∼S [σ(ϕ(x)⊤ θbt )ϕ(x)]. We have
Lemma 6. Let ∇L0t22b t (θbt )∇Lb t (θbt )⊤ ≼ Ex∼S [xx⊤ ] + Ex∼S [ϕ(x)(ϕ(x)⊤ ].
2∇L
0t
Proof of Lemma 6. For notation simplicity, we denote by at = Ex∼S0 [(σ(x⊤ θbt ) − 1)x] and bt =b t (θbt ), we have
Ex∼St [σ(x⊤ θbt )x]. Then, according to the definition of ∇L11(at + bt )(at + bt )⊤ ≼ (at at⊤ + bt b⊤t ),42
where the last inequality is due to Lemma 16. We can further bound the R.H.S. of the above inequality
by


⊤
at at⊤ = Ex∼S0 [(σ(ϕ(x)⊤ θbt ) − 1)ϕ(x)] Ex∼S0 [(σ(ϕ(x)⊤ θbt ) − 1)ϕ(x)]b t (θbt )∇Lb t (θbt )⊤ =∇L

≼ Ex∼S0 [(σ(x⊤ θbt ) − 1)2 ϕ(x)ϕ(x)⊤ ] ≼ Ex∼S0 [ϕ(x)ϕ(x)⊤ ].
29

In above the first inequality is due to the fact that E[a]E[a]⊤ ≼ E[aa⊤ ] for any random vectora ∈ Rd . The second inequality comes from (σ(x⊤ θbt ) − 1)2 ≤ 1 since the output value of the⊤sigmoid function σ is bounded in (0, 1). A similar arguments shows bt b⊤t ≤ Ex∼St [ϕ(x)ϕ(x) ].
Combining it with above displayed inequalities completes the proof.
Lemma 7. The predictions returned by the ONS algorithm satisfies,


X
|I|R2
2
b
b
∥∇Lt (θt )∥A−1 ≤ d log 1 +
.
t
λ
t∈I

Proof of Lemma 7. The proof follows the standard arguments in the analysis of the online Newton
step algorithm [13, Theorem 4.3]. We present the proof here for self-containedness.
XX
b t (θbt )∥2 −1 =b t (θbt )⊤ A−1
b b
∥∇L
∇L
t ∇Lt (θt )
A
t

t∈I

t∈I

=

X

=

X

≤

X



⊤
b
b
b
b
trace A−1
∇
L
(
θ
)∇
L
(
θ
)
t t
t t
t

t∈I

trace A−1t (At − At−1 )



t∈I

t∈I

log

|At |
|Ae |
= log
,
|At−1 |
|As |where the last inequality is due to the fact that trace(A−1 (A − B)) ≤ log(|A|/|B|) for any positive
definite matrix satisfying A ≽ B ≻ 0 [13, Lemma 4.5].
P
b t (θbt )∇Lb t (θbt )⊤ .
By definition, we have As = λId ∈ Rd×d and we have Ae = λId + t∈I ∇L
2dbb
Furthermore, since ∥∇Lt (θt )∥2 ≤ R, we have |Ae | ≤ (|I|R + λ) . In such a case, we can boundthe above inequality byX(|I|R2 + λ)d
|I|R2
|Ae |
2
b
b
∥∇Lt (θt )∥A−1 ≤ log≤ log≤ d log 1 +
t
|As |
λd
λ
t∈I

and finish the proof.
e t (θbt ) − ∇Lb t (θbt ), θbt − θ⟩. Then, with probability at least 1 − δ, we have
Lemma 8. Let Zt = ⟨∇L
2
2 
X
X
X 
⊤ b
⊤ b
b
b
b
b
Zt ≤ 2ν
Et ∇Lt (θt ) (θt − θ)
∇Lt (θt ) (θt − θ) + 2ν
t∈I

t∈I

t∈I

r
p
2|I| + 1
1
2|I| + 1+ log+ 4SR logνδ
δ2
for any δ ∈ (0, 1) and ν > 0.
Proof of Lemma 8. We can check that {Zt }et=s is a martingale difference sequence such that Et [Zt ] =e t (θbt ) − ∇Lb t (θbt )∥2 ∥θbt − θ∥2 ≤ 4SR,
0 since θbt is independent of St . Further noting that |Zt | ≤ ∥∇La direct application of Lemma 17 shows thatrp
XX
X   1
2|I| + 1
2|I| + 1
2
2
Zt ≤ ν
Zt + ν
Et Zt + log+ 4SR log
,
(40)
ν
δ
δ2
t∈I

t∈I

t∈I

with probability at least 1 − δ for any ν > 0. Then, for any t ∈ [T ], we can further bound the term
Zt2 and Et [Zt2 ] bye t (θbt ) − ∇Lb t (θbt ), θbt − θ⟩2 ≤ 2⟨∇Le t (θbt ), θbt − θ⟩2 + 2⟨∇Lb t (θbt ), θbt − θ⟩2
Zt2 = ⟨∇Landb t (θbt ), θbt − θ⟩2 ] − ⟨∇Le t (θbt ), θbt − θ⟩2 ,Et [Zt2 ] = Et [⟨∇L
30

(41)

where the first inequality is due to the fact (a − b)2 ≤ 2a2 + 2b2 for any a, b ∈ R. Then, combining
the above two displayed inequalities and taking the summation over T iterations, we haveX
X
X
X
b t (θbt ), θbt − θ⟩2 + 2Et [⟨∇Lb t (θbt ), θbt − θ⟩2 ].
Zt2 +Et [Zt2 ] ≤Zt2 + 2Et [Zt2 ] ≤ 2⟨∇L
t∈I

t∈I

t∈I

t∈I

(42)
We complete the proof by plugging (42) into (40).
Lemma 9. With probability at least 1 − δ, we havei
3X b b ⊤ b
5X h b b ⊤ b(∇Lt (θt ) (θt − θ))2 −Et (∇Lt (θt ) (θt − θ))2 ≤ ΛT (δ)
4
4
t∈I

t∈I

√
for any δ ∈ (0, 1), where ΛI (δ) = 32S 2 R2 log

2|I|+1+ 4S 2 R2δqlog 2|I|+1
δ2

Proof of Lemma 9. We can check {Yt }et=s is a martingale difference sequence such that |Yt | ≤4S 2 R2 . As a consequence, a direct application of Lemma 17 implies,
i
X
X
X h
b t (θbt )⊤ (θbt − θ))2b t (θbt )⊤ (θbt − θ))2 −Yt =(∇L
Et (∇L
t∈I

≤ 2ν

t∈I

X

t∈I


b t (θbt )⊤ (θbt − θ))4 + Et [(∇Lb t (θbt )⊤ (θbt − θ))4 ]
(∇L

t∈I

r
p
2|I| + 1
1
2|I| + 1
2 2+ log+ 4S R log
ν
δ
δ2


X
1
b t (θbt )⊤ (θbt − θ))2 + Et [(∇Lb t (θbt )⊤ (θbt − θ))2 ]
≤
(∇L
4
t∈I
r
p
2|I| + 1
2|I| + 1
2 2
2 2+ 4S R log,+ 32S R logδδ2
where the first inequality follows the proof of Lemma 8 and the second inequaliuty is due tob t (θbt )⊤ (θbt − θ))2 ≤ 4S 2 R2 .
the parameter setting ν = 1/(32S 2 R2 ) and the condition that (∇L
Rearranging the above inequlaity, we havei
3X b b ⊤ b
5X h b b ⊤ b(∇Lt (θt ) (θt − θ))2 −Et (∇Lt (θt ) (θt − θ))2
4
4
t∈I
t∈I
r
p
2|I| + 1
2|I| + 1
2 2
2 2≤ 32S R log+ 4S R log,δδ2which completes the proof.
D.3

Regret Guarantee for Meta-Algorithm (Lemma 10)

We have the following guarantee on the meta-algorithm, which ensures that for any interval Ii ∈ C,
the final prediction is comparable with the prediction of the base-learner Ei with an O(log T ) cost.
Lemma 10. For any interval Ii ∈ C and δ ∈ (0, 1), with probability at least 1 − δ, the final prediction
returned by our meta-learner satisfiesXXe t (θbt ) −e t (θbt,i ) ≤ O β log(T /δ) + β ln K ,
L
L
t∈Ii

t∈Ii

where K = |C| is the number of the intervals contained in the set C.
31

D.3.1Main Proof

This part presents the proof of Lemma 10. Our analysis for the meta-learner is based on that of thestrongly adaptive algorithm [70], where the problem of tracking base-learners on the corresponding
intervals is converted to a sleeping expert problem [71]. The main challenge of our problem is thatthe feedback function is noisy. Thus, we require to provide tight analysis on the generalization gap
b t and L
e t without ruining the O(log T ) regret guarantee. By twisting the meta-algorithm
between Lfrom Hedge [42] used in the previous work to Adapt-ML-Prod [43], we show that the generalization
gap can be cancelled by the negative term introduced by the exp-concavity of the loss functions (seehow to bound term (b) in (48)).
Proof of Lemma 10. To prove Lemma 10, we show that the final prediction θbt generated by ruleP
θbt = Ei ∈At pt,i θbt,i with the weight update procedure (25), (26) and (27) is identical to an algorithm
learning in the standard prediction with expert advice (PEA) setting [72], where all base-learners
submit their predictions at each time t ∈ [T ].
Lemma 11. For any base-learner Ei , let Ei′ be a surrogate of Ei , whose prediction θ̄t,i = θbt,i is thesame as that of Ei for any t ∈ Ii and θ̄t,i = θbt for other iterations. Then, for any time t ∈ [T ], thesurrogate meta-learner predicts asXε̄t−1,i v̄t−1,iθ̄t =p̄t,i θ̄t,i and p̄t,i = P,(43)
i∈[K] ε̄t−1,i v̄t−1,ii∈[K]where the potential v̄t,i and the learning rate ε̄t,i is defined asε̄t,i ε̄t−1,ifor all t ∈ [T ],v̄t,i = v̄t−1,i · 1 + ε̄t−1,i m̄t,i

(44)

and
(
ε̄t,i = min

1
,
2

s
1+

)

ln K
Pt

2
τ =1 m̄τ,jfor all t ∈ [T ](45)

b t (θ̄t ), θ̄t − θ̄t,i ⟩, v̄0,i = 1/K and ε̄0,i = min{1/2, ln K} for all i ∈ [K]. Then, wewith m̄t,i = ⟨Lhave θ̄t = θbt for any t ∈ [T ].
We have the following guarantee for the final prediction generated by the surrogate algorithm definedby (43).
Lemma 12. With probability at least 1 − 2δ, for any δ ∈ (0, 1), the prediction returned by thesurrogate meta-algorithm defined as (43) ensures,
T
T


X
X
e t (θ̄t ) −e t (θ̄t,j ) ≤ C (2) (δ) = O β log(T /δ) + β ln K
L
L
T
t=1

t=1
(2)

Ei′ ,

16(1+β)S 2 R2 (C ′ )2

1
T
for any base-algorithmi ∈ [K], where CT (δ) =+ 6(1+β)
+
ln K


q


√
2 2
R
S 2 R22CT′ + 24(1 + β) + 11S
log 2Tδ +1 + 4SR + 2(1+β)
log 2Tδ+1
and CT′ = 3 ln K +
2
1+β

K
ln 1 + 2e(1 + ln(T + 1)) = O(ln K + ln ln T ).

Then, according to the relationship between Ei and Ei′ , we can bound the regret of the final prediction
θbt with respect to that of any base-learner Ei over the interval Ii byXt∈Ii

e t (θbt ) −
L

X
t∈Ii

e t (θbt,j ) =
L

T
X
t=1

e t (θ̄t ) −LTX

e t (θ̄t,j ) ≤ C (2) (δ) = O(β log(T /δ) + β ln K).
L
T

t=1

In above, the first equality is due to Lemma 11 and the definition of the surrogate base-learner Ei′
such that θ̄t,i = θbt,i for any t ∈ Ii and θ̄t,i = θbt for any t ∈ [T ]/Ii . The first inequality is due toLemma 12. In the last equality, we treat double logarithmic factors in T as a constant, following
previous studies [73, 74]. We have completed the proof.
32

D.3.2

Useful Lemma

Proof of Lemma 11. We show that the predictions θbt and θ̄t by the two algorithm are exactly the samefor any t ∈ [T ] by induction. At iteration t ∈ [T ], we have the following two induction hypotheses:• IH1: the final predictions satisfy θbτ = θ̄τ for all τ ≤ t.
• IH2: any base-algorithms Ei ∈ Aτ satisfy pτ,i = p̄τ,i /(P

i∈Aτ p̄τ,i ) for all τ ≤ t.

Base case: For the base case of t = 1, there is only one active base-algorithm. Denote by i1 theindex of the active base-algorithm, we have θb1 = θb1,i1 , p1,i1 = 1 and A1 = {i1 }. Moreover, for thePsurrogate algorithm, we have p̄1,i = 1/K for any i ∈ [K] and θ̄1 = ( i∈[K]/{i1 } θb1 + θb1,i1 )/|K| =P
θb1 . We can check that θb1 = θ̄1 and p1,i = p̄1,i /
p̄1,i for any i ∈ A1 .
i∈A1Induction step: For the induction step, we first showp̄t+1,ipt+1,i = P
i∈At+1 p̄t+1,i(46)

for all i ∈ At+1 . By the definition of p̄t+1,i , it is sufficient to prove
ε̄t,i v̄t,ipt+1,i = P
for any i ∈ At+1 .
i∈At+1 ε̄t,i v̄t,i(47)

For the base-algorithm Ei active at iteration t + 1 and his/her corresponding surrogate base-algorithm
Ei′ , we can decompose the time horizon until t + 1 into two intervals: asleep interval Iislp = [1, si − 1]and active part Iiact = [si , t + 1]. For any time stamp belonging to the asleep interval τ ∈ Iislp , wecan check that the loss for the base-algorithm Ei′ satisfiesb τ (θ̄τ ), θ̄τ − θ̄τ,i ⟩ = ⟨Lb τ (θbτ ), θbτ − θbτ ⟩ = 0,m̄τ,i = ⟨L
where the first equality is by the definition of m̄τ,i . The second equality is due to IH1 such thatθ̄τ = θbτ for any τ ≤ t and the definition of the surrogate base-learner Ei′ such that the predictionθ̄τ,i = θbτ on the asleep interval. For any time stamp belonging to the active interval τ ∈ Iiact , similarargument shows thatb τ (θ̄τ ), θ̄τ − θ̄τ,i ⟩ = ⟨Lb τ (θbτ ), θbτ − θbτ,i ⟩ = mτ,i ,m̄τ,i = ⟨L
Thus, by definition, we can draw the conclusion that v̄t,i = vt,i and ε̄t,i = εt,i for any active
base-algorithm with index i ∈ At+1 , which would finally lead to (47).
Then, we show θbt+1 = θ̄t+1 based on IH1 and IH2. We haveXXθ̄t+1 =
p̄t+1,i θbt+1,i +p̄t+1,i θbt+1
i∈At+1i∈[K]/At+1




X

= 

p̄t+1,i  

i∈At+1X

pt+1,i θbt+1,i  + 1 −i∈At+1




= 


X


X

p̄t+1,i  θbt+1

i∈At+1Xp̄t+1,i  θbt + 1 −i∈At+1p̄t+1,i  θbt+1 = θbt ,i∈At+1

where the first equality is due to (46) and the second inequality is due to the definition of θbt+1 . We
have finished the induction steps and completed the proof.
Proof of Lemma 12. Since θbt and θbt,i are both independent of St , using the same arguments for
obtaining (33), we can obtain thatTX
t=1

e t (θbt ) −
L

T
X
t=1

e t (θbt,i ) ≤
L

T
T
X
X
b
b
b
e t (θbt ) − ∇Lb t (θbt ), θbt − θbt,i ⟩b⟨∇Lt (θt ), θt − θt,i ⟩ +
⟨∇L
t=1

|

t=1

{z

term (a)

33

}

|

{z

term (b)

}

−

T
X

h
 i
1
b t (θbt )⊤ (θbt − θbt,i ) 2
Et ∇L2(β + 1)t=1

(48)

For term (a), we have
v
u
T 
2
X
2SRCT′ ut1 +b t (θbt )⊤ (θbt − θbt,i ) + 2CT′∇L
term (a) ≤ √
ln K
t=1
T 
2
X
1
1
24(β + 1)S 2 R2 (CT′ )2b t (θbt )⊤ (θbt − θbt,i ) + 2CT′ ,
+
+
∇L
ln K
6(β + 1) 6(β + 1) t=1(49)

K
for any i ∈ [K], where CT′ = 3 ln K + ln 1 + 2e
(1 + ln(T + 1)) = O(ln K + ln ln T ). In above,
the first inequality is a direct application of [43, Corollary 4] with the linearized loss Mt (θ) =b t (θbt ), θ⟩/(2SR) + 1/2 and the second inequality is due to the AM-GM inequality.
⟨L≤

For term (b), by Lemma 8 with ν = 6(1 + β) and a similar argument in (36), we have

T 
2 
X

1
b t (θbt )⊤ (θbt − θbt,i ) 2 + Et ∇L
b t (θbt )⊤ (θbt − θbt,i )
term (b) ≤∇L
12(1 + β) t=1
r√
2T + 1
2T + 1+ 4SR log(50)
+ 24(1 + β) log
δ
δ2
Combining the upper bound for term (a) and term (b) with (48) yields,TX

e t (θbt ) −
L

t=1

T
X

e t (θbt,i )
L

t=1
T 
X


2
2 
1
5
⊤ b
⊤ b
b
b
b
b
b
b
≤
∇Lt (θt ) (θt − θt,i ) −Et ∇Lt (θt ) (θt − θt,i )4(1 + β) t=1
12(1 + β)|{z}term (c)r√
2T + 124(1 + β)S1
2T + 1′+++ 2CT + 24(1 + β) log+ 4SR log.
ln K
6(1 + β)δδ2Then, according to Lemma 9, we can further bound term (c) by

T 
2
2 
X
1
5
b t (θbt )⊤ (θbt − θbt,i ) −b t (θbt )⊤ (θbt − θbt,i )
term (c) =∇L
Et ∇L
4(1 + β) t=1
12(1 + β)rp
2|I| + 1
11S 2 R2
4S 2 R2
2T + 1
≤
log
+
log
,
1+β
δ
3(1 + β)δ2which implies
2

T
X
t=1

R (CT′ )22

e t (θbt ) −LTX

e t (θbt,i )L

t=124(1 + β)S 2 R2 (CT′ )2
1≤++ 2CT′ + C1′ log
ln K
6(1 + β)= O(β log(T /δ) + β ln K).

√

2T + 1
+ C2′
δ

r
log

2T + 1
δ2

In above, C1′ = 24(1 + β) + 11S 2 R2 /(1 + β) and C2′ = 4SR + 4S 2 R2 /(3(1 + β)). In the lastequality, we treat double logarithmic factors in T as a constant, following previous studies [73, 74].
We have completed the proof.
D.4

Proof of Theorem 2This part presents the proof for Theorem 2, which is based on the dynamic regret for the base-learner(Lemma 3) and that for the meta-learner (Lemma 10) presented in Appendix D.1.
34

D.4.1Main Proof

Proof of Theorem 2. When choosing the Bregman divergence function as ψ(t) = ψLR (t) ≜ t log t −(t + 1) log(t + 1) and the hypothesis space HθLR = {x 7→ exp(−θ⊤ ϕ(x)) | ∥θ∥2 ≤ S}. The expected
loss function Lt (θ) as defined by (4) is equal to⊤⊤1
Lt (θ) =ED0 [log(1 + e−ϕ(x) θ )] + EDt [log(1 + eϕ(x) θ )] .
2
Let θt∗ ≜ arg minθ∈Rd Lt (θ). Since the ground-truth density ratio rt∗ (x) = Dt (x)/D0 (x) function
is contained in the hypothesis space as Dt (x)/D0 (x) ∈ HθLR . One can show that Dt (x)/D0 (x) =exp(−(θt∗ )⊤ ϕ(x)) and θt∗ ∈ Θ = {θ ∈ Rd | ∥θ∥2 ≤ S}.
Then, for the base-learner Ei running on any interval Ii , Lemma 3 in Appendix D.2 shows that its
prediction θbt,i is comparable with the best model at every iteration θt∗ with an O(|Ii |(VIθi )2 ) cost:


X
X


e t (θbt,i ) −e t (θt∗ ) ≤ O dβ log |Ii |/δ + |Ii | VIθ 2 + |Ii | log(dT /δ) ,
L
L
(51)
i
N0
t∈Ii

t∈Ii

where VIθi =

Pei

∗
∗
t=si +1 ∥θt − θt−1 ∥2 is the path length measuring the fluctuation of the minimizers.

Besides, Lemma 10 in Appendix D.3 ensures that the prediction θbt of the meta-learner is comparablewith the prediction θbt,i of each base-learner Ei on the associated interval Ii with an O(log T ).
XX
e t (θbt ) −e t (θbt,i ) ≤ O β log(T /δ) + β ln K .
L
L
(52)
t∈Ii

t∈Ii

A direct combination of (51) and (52) shows that the meta-learner is able to track the best prediction
θt∗ on any interval Ii ∈ C. Actually, by exploiting the structure of the geometric covering (24), it canbe shown that such a property can be extended to arbitrary interval I ∈ [T ] as the following lemma,
whose proof is deferred to Appendix D.4.2.
Lemma 13. With probability at least 1 − 5δ for any δ ∈ (0, 1), for any interval [s, e] ∈ [T ], themeta-learner’s prediction θbt returned by (9) ensures,XX
e t (θbt ) −e t (θ∗ )
L
L
t
t∈I
2

t∈I




 16|I| ln((d+1) log|I|T/δ)R|I| log(dT /δ)(1)|I|(VIθ )2 + 2 log|I| CT (δ ′ ) + CT2 (δ ′ ) +
+o
4
N0
N0
P
(1)
e
∗
where VIθ = t=s+1 ∥θt−1− θt∗ ∥2 and δ ′ = δ/2 log|I|. The coefficient CT (δ) = O (dβ log(T /δ))

(2)and CT (δ) = O β log(T /δ) is defined in Lemma 4 and Lemma 12, respectively.
≤Then, we can use the bin partition [19, Lemma 30] to convert the dynamic regret on each interval to a
dynamic regret over T iterations. For self-containedness, we restate the lemma here.
Lemma 14 (Lemma 30 of [19]). There exists a partition P of the time horizon [T ] into M =O(max{dT 1/3 UTθ )2/3 , 1} intervals ,i.e., {Ii = [si , ei ]}M
i=1 such that for any interval Ii ∈ P,
p
Pei
∗
∗
∗
UIθi ≤ Dmax / |Ii |, where UIθi =∥θ−θt−1 ∥1 and Dmax = maxt∈[T ] ∥θt ∥∞ ≤t=si +1 t
maxt∈[T ] ∥θt∗ ∥2 ≤ S.
Then, we can decompose the overall dynamic regret of our method into those over the key partition
P. Specifically, with probability at least 1 − 5δ, we haveTX

e t (θbt ) −
L

t=1

T
X

e t (θt∗ )
L

t=1

!
=

X

X

Ii ∈P

t∈Ii

e t (θbt ) −LXe t (θt∗ )
L

Ii

35

≤

X 16|Ii | ln

X R2



2(d+1)T log|Ii |M/δ

|Ii |(VIθi )2 +
4
N0
Ii ∈P
Ii ∈P





X
(1)
+2
log|Ii | CT δ/2M log|Ii | + CT2 δ/2M log|Ii |
Ii ∈P
2

≤

2

16T ln

MS R
(1)
(2)
+ 2M log T CT δ/2M log T + CT δ/2M log T +4{z} || {z } |










2(d+1)T log T M/δ


,
}

N0
{z

term (b)term (a)term (c)where the second last inequality is due to Lemma 13 and a union bound taken over all intervals√
in the key partition. The last inequality is due to Lemma 14 such that VIθi ≤ UIθi ≤ S/ Ii
for any interval Ii ∈ P. Lemma 14 shows that we can bound the partition number by M =O(max{dT 1/3 (UTθ )2/3 ), 1}) = O(max{d4/3 T 1/3 (VTθ )2/3 , 1}). Then, we can further show that
term (a) ≤ O d4/3 max{T 1/3 (VTθ )2/3 , 1} ,term (b) ≤ O d4/3 β log T log(T /δ) max{T 1/3 (VTθ )2/3 , 1} ,T ln(dT /δ)term (c) ≤ O,
N0which implies that the overall dynamic regret is bounded byTX

e t (θbt ) −
L

t=1

T
X




1
4
T ln(dT /δ)
θ 23
∗
e
3
3
Lt (θt ) ≤ O d β log T log T /δ · max{T (VT ) , 1} +
.
N0
t=1

e t (θt∗ ) = L
e t (rt∗ ) under the realizable assumption andHence, we complete the proof by showing L
PT∗
converting the variation of the model parameter VTθ = t=2 ∥θt∗ − θt−1
∥2 to the variation of thePTfeature distribution VT = t=2 ∥Dt (x) − Dt−1 (x)∥1 by Lemma 15.
D.4.2

Useful Lemmas

Lemma 13. With probability at least 1 − 5δ for any δ ∈ (0, 1), for any interval [s, e] ∈ [T ], themeta-learner’s prediction θbt returned by (9) ensures,XX
e t (θbt ) −e t (θ∗ )
L
L
t
t∈I

t∈I



 16|I| ln((d+1) log|I|T/δ)
R2|I| log(dT /δ)
(1) ′
θ 2
2 ′
≤
|I|(VI ) + 2 log|I| CT (δ ) + CT (δ ) +
+o
,
4
N0
N0
Pe
(1)
∗
where VIθ = t=s+1 ∥θt−1− θt∗ ∥2 and δ ′ = δ/2 log|I|. The coefficient CT (δ) = O (dβ log(T /δ))

(2)and CT (δ) = O β log(T /δ) is defined in Lemma 4 and Lemma 12, respectively.
Proof of Lemma 13. A direct combination of Lemma 3 and 10 shows that, with probability at least
1 − 5δ, we have the following dynamic regret guarantee on any interval Ii = [si , ei ] ∈ C.
X
X
e t (θbt ) −e t (θt∗ )
L
L
t∈Ii

t∈Ii
2

16|Ii | ln((d + 1)|Ii |/δ)|Ii |(VIθi )2 +
+o
4
N0

R
(1)
(2)
≤ CIi (δ) + CT (δ) +|Ii | log(d|Ii |/δ)N0

, (53)(1)∗− θt∗ ∥2 . The coefficient CIi = O (dβ log(|Ii |/δ)) is defined in
t=si +1 ∥θt−1

(2)
Lemma 4 and CT (δ) = O β log(T /δ) is defined in Lemma 12.

where VIθi =

Pei

36

For any interval I = [s, e] ⊆ [T ], [75, Lemma 3] showed that it can be partitioned into two sequencesof disjoint and consecutive intervals, denoted by I−p , . . . , I0 ∈ C and I1 , . . . , Iq ∈ C, such that
|I−i |/|I−i+1 | ≤ 1/2, ∀i ≥ 1

and|Ii |/|Ii−1 | ≤ 1/2, ∀i ≥ 2.

(54)

Then, with probability at least 1 − 5δ, we can decompose the regret over I = [s, e] ⊆ [T ] as
!
q
X
X
X
X
X
∗
∗
b
b
e
e
e
e
Lt (θt ) −Lt (θt ) =Lt (θt ) −Lt (θt )
t∈I

i=−p

t∈I

t∈Ii

t∈Ii
q

R2 X
|Ii |(VIθi )2
4 i=−pqqX
X
16|Ii | ln((d+1)(p+q)T/δ)|Ii | log(dT (p + q)/δ)
+
+
o
.
N0
N0
i=−p
i=−p
(1)

(2)

≤ (p + q)(CT (δ/p+q) + CT (δ/p+q)) +where the last inequality is due to (53) and a union bound taken over all intervals. Then, due to thepartition of the intervals (54), we can check that p + q ≤ 2 log|I|. Then, we can further bound theabove displayed inequality by


2
X
X
e t (θbt ) −e t (θt∗ ) ≤ 2 log|I| C (1) (δ/(2 log |I|)) + C (2) (δ/(2 log |I|)) + R |I|(VIθ )2
L
L
T
T
4
t∈I
t∈I


|I| log(dT log|I|/δ)
16|I| ln(2(d+1) log|I|T/δ)
+
+o
N0
N0


2
R
|I| log(dT log |I|/δ)θ 2
T
log|I|=|I|(VI ) + O dβ log|I| log(
/δ) +
,
4
N0
Pq
where the first inequality is due to the fact that i=−p VI2i ≤ VI2 . This completes the proof.
Lemma 15. Under the condition that Dt (x)/D0 (x) ∈ HθLR ≜ {exp(−θ⊤ ϕ(x)) | ∥θ∥2 ≤ S} forany t ∈ [T ] and Ex∼D0 [|ϕ(x)⊤ a|] ≥ α∥a∥2 for any a ∈ Rd with a certain α > 0. Then, we have∗∥θt−1
− θt∗ ∥2 ≤ β∥Dt−1 (x) − Dt (x)∥1 /α.
Proof of Lemma 15. Under the realizable assumption such that Dt (x)/D0 (x) ∈ HθLR . One canverify that exp(−ϕ(x)⊤ θt∗ ) = Dt (x)/D0 (x). Then, we can convert the variation of the model∗∥θt−1− θt∗ ∥2 to the variation of the distribution ∥Dt−1 (x) − Dt (x)∥1 .


1
∗
∗
∥θt−1
− θt∗ ∥2 ≤ Ex∼D0 |ϕ(x)⊤ (θt−1− θt∗ )|
α



⊤ ∗
⊤ ∗ 
Dt−1 (x)βββDt (x)−= ∥Dt (x) − Dt−1 (x)∥1 ,≤ Ex∼D0 |eϕ(x) θt−1 − eϕ(x) θt | ≤ Ex∼D0αα
D0 (x)
D0 (x)
αwhere the first inequality is by the condition Ex∼D0 [|ϕ(x)⊤ a|] ≥ α∥a∥2 and the second is by the∗mean value theorem and the boundedness of |ϕ(x)⊤ (θt−1− θt∗ )| ≤ 2SR for any x ∈ X . The lastinequality is due to the realizable assumption.

D.5

Proof of Theorem 3

Proof of Theorem 3. Then, we can bound the overall excess risk over the model trained by IWERM (1) by a combination of Proposition 1 and Theorem 1 as
!
T
T
T
X


1 X
B · CT (δ) 2L X∗√
b t) −Rt (w
Rt (wt ) ≤+Ex∼S0 |brt (x) − rt∗ (x)|
T t=1
T
N
0
t=1
t=1
v
u  T
√ !

T
X
X
d
B · CT (δ) 4Lβ u
log(T
/δ)
t
∗
e t (θbt ) −e t (θt ) + O
√
≤ √
+
T
L
L
,
(55)
T
N0
N
0
+
t=1
t=1
37

p
where the first inequality is due to Proposition 1 and CT (δ) = 2GDR + 5L 2 log(8T /δ) =O(log(T /δ)). The second inequality is due to Theorem 1 with the choice of the logistic regression
model ψ(t) = t log t−(1+t) log(1+t) and rbt (x) = exp(−ϕ(x)⊤ θbt ) and the realizability assumptionrt∗ (x) = exp(−ϕ(x)⊤ θt∗ ). Then, by Theorem 2, we have


T
T
X
X
2

1
4
T ln(dT /δ)
∗
3
e
e
3
3
Lt (b
rt ) −Lt (rt ) ≤ O d log T log T /δ · max{T VT , 1} +
,
N0
t=1
t=1
which implies
v "
s
#
u


T
T
X
u X
1 √

2
2
ln(dT /δ)
∗
3
tT
b
e
e
3
3
Lt (θt ) −≤ O d log T /δ · max{T VT , T } + T.
Lt (θt )
N0
t=1
t=1
+

(56)
We complete the proof by substituting (56) it into (55).
D.6

Discussion on Tightness of the Bound

This part illustrates the tightness of our bound with the case where the labels of testing samples areavailable after the prediction.
Continuous Shift with Labeled Feedback. We consider a T -round online learning process. At
b t ∈ W. At the same time, the environments
iteration t, the learner submits her prediction w
pick the data pair (xt , yt ) ∼ Dt . Then, the learner updates her model with the loss function
bt (w) = ℓ(w⊤ xt , yt ) and obtains the prediction w
b t+1 ∈ W for the next iteration. The goalRof the learner is to minimize the cumulative excess risk against the sequence of optimal models
wt∗ ∈ arg minw∈W Rt (w), that is,b t }Tt=1 ) =RT ({w

T
T
X
1Xb t) −Rt (w
Rt (wt∗ ).
T t=1
t=1

(57)

The goal is the same as that of standard continuous covariate shift, see the problem setup in (1). The
bt (w)
key difference is that now the label of testing data is available. Therefore, the noisy feedback R
observed by the learner is an unbiased estimator with respect to the true risk Rt (w).
Tightness of Our Bound. The continuous shift with labeled feedback is essentially a non-stationary
stochastic optimization problem [39]. For general convex functions, [39] showed that any gradientbased algorithm will sufferb t }Tt=1 )] = Ω T −1/3 (VTL )1/3
E[RT ({w(58)
PTin the worst case. In above, VTL = t=2 maxw∈W |Rt (w) − Rt−1 (w)| measures the fluctuation of
the risk function. For the same performance measure, our algorithm achievese √1 + max{T −1/3 V 1/3 , T − 21 } .
b t }Tt=1 ) = O
RT ({w
(59)
T
N0
 √

1
e 1/ N0 + T −1/3 V 1/3 . As
In the non-stationary case, i.e., VT ≥ Θ(T − 2 ), our bound becomes OT√ediscussed below Proposition 1, the first term O(1/N0 ) is hard to be improved. We thus focus one −1/3 V 1/3 ).
the tightness of the second term O(TT

Although the definition of VT in our upper bound is different from the that of VTL in the lower
bound (58), the two bounds share the same dependence on the time horizon T , which provides
evidence that our bound is hard to improve. Indeed, consider the 1-dimensional case where theunderlying distribution only shifts once from D1 (x, y) to D2 (x, y) at a certain time t ∈ [T ]. In sucha case, VT = O(1) and our bound implies an O(T −1/3 ) rate for the second term of (59). On theother hand, the rate in the lower bound (58) becomes Ω(T −1/3 ). The same dependence on the time
horizon T indicates our bound is hard to improve. We leave a precise lower bound argument as the
future work.
38

D.7

On Other Choices of Bregman Divergence Function

In this section, we discuss how to apply the other two choice of the divergence function ψLS =(t − 1)2 /2 and ψKL = t log t − t in our framework. Our analysis crucially relies on two conditions:
1. the Bregman divergence function ψ satisfies the conditions in Theorem 1.
b t (θ) is an exp-concave and smooth function with respect to θ.
2. the induced loss function L
Choice of ψLS (t) = (t − 1)2 /2. Considering the Bregman divergence function ψ(t) = (t − 1)2 /2
and the hypothesis space Hθ = {x 7→ ϕ(x)⊤ θ | ∥θ∥2 ≤ S}, where ϕ : X → Rd represents abψ
specific basis function with bounded norm ∥ϕ(x)∥2 ≤ R. Then, the loss Lt (θ) as per (7) becomesb t (θ) = 1 ES [(ϕ(x)⊤ θ)2 ] − ES [ϕ(x)⊤ θ].
Lt
2 0
The above configuration recovers the unconstrained least-squares importance fitting (uLSIF)
method [28]. We can show such a choice of divergence function satisfy the condition requiredby Theorem 1 and enjoys favorable function properties.
• Conditions in Theorem 1: when Bregman divergence function is chosen as ψLS , wehave ∂ 2 ψLS (t) = 1 and ∂ 3 ψLS (t) = 0. Thus, ψLS is a 1-strongly convex function andt∂ 3 ψLS (t) = 0 for t ∈ R.
b t (θ) = ES [ϕ(x)ϕ(x)⊤ ]. Then, when the• Expconcavity and smoothness: We have ∇2 L
0
b t (θ) ≼ R2 Id ,input is upper bounded by ∥ϕ(x)∥2 ≤ R for all ϕ(x) ∈ X , we have ∇2 L
which implies the smoothness. Furthermore, if the initial data ensure ES0 [ϕ(x)ϕ(x)⊤ ] ≽b t (θ) is strongly convex and thus is exp-concave as shown by [60, Lemma 2]. We
αId , the Lnote that the regularity condition on the offline data ES0 [ϕ(x)ϕ(x)⊤ ] ≽ αId is used in theanalysis of LSIF algorithm [28, Lemma 1].
Choice of ψKL = t log t − t. Considering the Bregman divergence function ψ(t) = t log t − t andthe hypothesis Hθ = {x 7→ ϕ(x)⊤ θ | ∥θ∥2 ≤ S}, where ϕ : X → Rd+ represents a specific basisbψ
function with bounded norm ∥ϕ(x)∥2 ≤ R. Then, the loss Lt (θ) as per (7) becomesb t (θ) = ES [ϕ(x)⊤ θ)] − ES [log(ϕ(x)⊤ θ)].
L0tb t (θ) is well-defined, we further require the density function rb(x) = ϕ(x)⊤ θ > 1/β
To ensure L
for any x ∈ X and θ ∈ Θ ≜ {θ | ∥θ∥2 ≤ S} with a certain positive constant β > 0. The above
configuration recovers the UKL method [37] with the generalized linear model. We can show sucha choice of divergence function satisfy the condition required by Theorem 1 and enjoys favorable
function properties.
• Conditions in Theorem 1: When Bregman divergence function is chosen as ψKL , we have∂ 2 ψKL (t) = 1/t and ∂ 3 ψKL (t) = −1/t2 . When the input of the loss function satisfiest ≥ β, ψKL (t) is a 1/β-strongly convex function, and we can also check that t∂ 3 ψKL (t) ≤ 0
for all t ∈ R+ .
b t (θ) = ES [ϕ(x)ϕ(x)⊤ /(ϕ(x)⊤ θ)2 ]. Then,• Expconcavity and smoothness: We have ∇2 Ltsince the output of the density ratio function rb(x) = ϕ(x)⊤ θ ≥ 1/β for any x ∈ X
b t (θ) ≼and θ ∈ Θ and the boundedness of the input feature ∥ϕ(x)∥2 ≤ R, we have ∇2 L
2 2
β R Id for any θ ∈ Θ. Then, the loss function is a smooth function. Besides, under the
regularity condition ES0 [xx⊤ ] ≽ αId of the initial data and again the boundedness of Θb t (θ) ≽and ∥ϕ(x)∥2 ≤ R, the loss function is also a strongly convex function by ∇2 L
2α/(SR) Id , which indicates the expconcavity of the loss function.

E

Technical Lemmas

This section presents several useful technical lemmas used in the proof.
39

Lemma 16. For any a, b ∈ Rd , we have (a + b)(a + b)⊤ ≼ 2(aa⊤ + bb⊤ ), where for any matrix
A, B ∈ Rd×d , A ≼ B indicates B − A is a positive semi-definite matrix.

Proof. For any x ∈ Rd , x⊤ 2(aa⊤ +bb⊤ )−(a+b)(a+b)⊤ x = x⊤ aa⊤ +bb⊤ −ab⊤ −ba⊤ x= x⊤ (a − b)(a − b)⊤ x ≥ 0, which completes the proof.
Lemma 17 (Theorem 4 of [68]). Let {Zi : i ≥ 1} be a martingale difference with the filtration
F = {Fn : n ≥ 1} and suppose |Zi | ≤ R for all i ≥ 1. Then, for any δ ∈ (0, 1), r > 0, withprobability at least 1 − δ,
!
r
√
t
t
t
X
X
X
2t + 12t + 1
1
2
2.
Zi ≤ rZi +E[Zi | Fi−1 ] + log+ R log
r
δ
δ2
i=1
i=1
i=1

40