# HUBER ADDITIVE MODELS FOR NON-STATIONARY TIME SERIES ANALYSIS

**Yingjie Wang[1][,][2], Xianrui Zhong[3], Fengxiang He[2][,][∗], Hong Chen[4][,][∗]& Dacheng Tao[2]**

1College of Informatics, Huazhong Agricultural University, China
2JD Explore Academy, JD.com Inc, China
3Department of Computer Science, University of Illinoist at Urbana-Champaign, USA
4College of Science, Huazhong Agricultural University, China
yingjiewang1201@gmail.com, xzhong23@illinois.edu,
fengxiang.f.he@gmail.com, chenh@mail.hzau.edu.cn,
dacheng.tao@gmail.com

ABSTRACT

Sparse additive models have shown promising flexibility and interpretability in processing time series data. However, existing methods usually assume the time series
data to be stationary and the innovation is sampled from a Gaussian distribution.
Both assumptions are too stringent for heavy-tailed and non-stationary time series
data that frequently arise in practice, such as finance and medical fields. To address
these problems, we propose an adaptive sparse Huber additive model for robust
forecasting in both non-Gaussian data and (non)stationary data. In theory, the
generalization bounds of our estimator are established for both stationary and nonstationary time series data, which are independent of the widely used mixing conditions in learning theory of dependent observations. Moreover, the error bound for
non-stationary time series contains a discrepancy measure for the shifts of the data
distributions over time. Such a discrepancy measure can be estimated empirically
and used as a penalty in our method. Experimental results on both synthetic and
real-world benchmark datasets validate the effectiveness of the proposed method.
[The code is available at https://github.com/xianruizhong/SpHAM.](https://github.com/xianruizhong/SpHAM)

1 INTRODUCTION

Additive model has become one of the most powerful tools for time series analysis due to the
exemplary monograph (Stone, 1985; Hastie & Tibshirani, 1990) and companion software (Chambers
& Hastie, 1992). For the past two decades, the growing importance of algorithmic flexibility
and interpretability motivates the development of various additive models along with theoretical
explorations (Huang & Yang, 2004; Wang & Yang, 2007; Chu & Glymour, 2008; Song & Yang, 2010;
Yang et al., 2018; Chen et al., 2017; Liu et al., 2020; Chen et al., 2021a) and practical applications
(Dominici et al., 2002; Wang & Brown, 2011; Ravindra et al., 2019; Bussmann et al., 2020; Wang
et al., 2020). Although these aforementioned works have shown promising behaviours, the proposed
methods in these works require some stringent assumptions on the stochastic process, e.g., various
mixing conditions (Doukhan, 1994), stationary distribution and Gaussian innovation.

A number of attempts have been made to relax such stringent assumptions. For the purpose of dealing
with non-Gaussian innovation, Qiu et al. (2015) develop an elliptical vector autoregressive model
for estimating heavy-tailed stationary processes with parametric convergence analysis that reduces
the influence of non-Gaussian innovation. Moreover, under stationarity and decaying β-mixing
condition, Wong et al. (2020) derive nonasymptotic estimation error of lasso without assuming
special parametric form of the data generating process. Stationarity and various mixing conditions
are commonly adopted in many previous studies, see, e.g., (Mohri & Rostamizadeh, 2009; 2010;
Kock & Callot, 2015; Wong et al., 2020). In practice, the mixing and stationary conditions are too
stringent and not always valid (Baillie, 1996; Kuznetsov & Mohri, 2020a). To relax the mixing and
stationary conditions, Adams & Nobel (2010) prove asymptotic guarantees for stationary ergodic

_∗Corresponding authors._


-----

sequences. Agarwal & Duchi (2013) establish generalization bounds for asymptotically stationary
(mixing) processes in the case of stable on-line learning algorithms. Kuznetsov & Mohri (2014)
establish learning guarantees for fully non-stationary and mixing processes. Recently, Kuznetsov &
Mohri (2020a) provide data-dependent generalization risk bounds for non-stationary and non-mixing
stochastic processes.

However, the exploration of additive models for both robust and non-stationary time series forecasting
is still very limited. In this paper, we propose a class of sparse Huber additive models with theoretical
guarantees. Our main contributions are summarized as follows:

-  Algorithm design and theoretical guarantees: In Section 3.1, we first propose a novel sparse
Huber additive model (SpHAM) with Huber loss, sparsity-inducing ℓ2,1-norm regularizer
and additive data-dependent hypothesis space. This proposed method works for stationary
time series and can achieve robust forecasting and satisfactory inference (e.g., Granger
causal discovery) simultaneously. In theory, Section 3.2 establishes the upper bound of the
function approximation error of SpHAM by developing error decomposition technique (Wu
et al., 2006; Chen et al., 2021b) and employing sequential Rademacher complexity (Rakhlin
et al., 2010; 2015; Kuznetsov & Mohri, 2020a). With properly selected scale parameters, the
theoretical findings indicate that: a) For stationary time series, the function approximation
consistency with convergence rate O(T _[−]_ 2[1] ) can be pursued, even if the innovation is non
Gaussian distribution (see Theorem 1 and Corollary 1 for more details). Moreover, this
consistency analysis appears to be novel because no explicit data dependence assumptions
are not imposed here, e.g., various mixing conditions used in (Doukhan, 1994; Mohri &
Rostamizadeh, 2010; Zou et al., 2009; Wong et al., 2020); b) For non-stationary time series,
the function approximation error is bounded by a discrepancy measure, which characterizes
the drifts of the data distributions along with time (see Theorem 2 for more details). By
penalizing such a discrepancy measure, Section 3.3 further proposes an adaptive SpHAM
for non-stationary time series and provides its theoretical upper bound correspondingly (see
Theorem 3 for more details).

-  Optimization and empirical evaluations: The proposed SpHAM and adaptive SpHAM can
be implemented efficiently by Fast Iterative Shrinkage-Thresholding Algorithm (FISTA)
(Beck & Teboulle, 2009). Experimental results on both synthetic and real-world benchmark
CauseMe (Runge et al., 2019) validate the effectiveness of the proposed method.

**Related works: There exist some interesting studies towards sparse additive model for time series**
analysis from both theoretical and practical viewpoints; see, e.g., Chu & Glymour (2008); Song &
Yang (2010); Yang et al. (2018). Although they show promising interpretability and modeling capacity,
all of them require Gaussian innovation, stationarity and various mixing conditions. Moreover, these
assumptions are not always valid in heavy-tailed and non-stationary time series data. In contrast to
these previous additive models, we are seeking to formulate our method and investigate its asymptotic
properties without resorting to such strict assumptions.

Robust forecasting is a vibrant area of research in time series analysis (Qiu et al., 2015; Wong et al.,
2020). As one of the triumphs and milestones of robust statistics, the success stories of Huber
algorithms are mainly on robust prediction tasks with i.i.d datasets (Huber, 1964; Huber & Ronchetti,
2009; Loh, 2017; Feng & Wu, 2020), and their extensions to time series prediction are fairly sparse.
To our best knowledge, this is the first work that considers Huber additive models for non-stationary
and dependent time series data.

Our theoretical analyses are inspired from the successful usage of sequential Rademacher complexity
in Kuznetsov & Mohri (2020b;a). The generalization risk bounds are established in their works
for a general scenario of non-stationary and non-mixing stochastic processes. However, it should
be mentioned that we are chiefly concerned about the function approximation analysis, which is
essentially different from theirs and is very crucial for Huber algorithms, since the convergence
of generalization risk cannot imply the convergence of function approximation (Sun et al., 2019;
Feng & Wu, 2020; He & Tao, 2020). Moreover, the development of analysis techniques (e.g., error
decomposition and sequential Rademacher complexity) for assessing our method may shed light on
other robust models for nonstationary time series analysis.

To better highlight the novelty of our method, Table 1 summarizes the algorithmic properties of our
method and other related works, e.g., Sparse additive models for time series (TS-SpAM) (Yang et al.,


-----

|Col1|Table 1: The properties of related methods|
|---|---|
||TS-SpAM DBF R-Dantzig Ours|
|Hypothesis space|Additive Additive Kernel-based Linear Spline-based Kernel-based|
|Loss|Squared Squared Quantile-based Huber|
|Mixing condition|Yes No Yes No|
|Stationarity|Yes No Yes No|
|Robustness|No No Yes Yes|
|Sparsity|Yes No Yes Yes|


2018), Discrepancy-based Forecasting (DBF) (Kuznetsov & Mohri, 2020a) and Robust Dantzigselector-type estimator (R-Dantzig) (Qiu et al., 2015).

The remainder of this paper is organized as follows. Section 2 recalls the background of additive
models. Section 3 mainly provides our methods and theoretical guarantees. We provide empirical
evaluations in Section 4. Finally, Section 5 concludes this paper. The source code package is available
[at https://github.com/xianruizhong/SpHAM.](https://github.com/xianruizhong/SpHAM)

2 PRELIMINARY

Let {Z _[t]}t[∞]=−∞_ [be a stochastic time series with time index][ t][, where variable][ Z] _[t][ = (][X]_ _[t][, Y][ t][)][ takes]_
values in the compact input space X ⊂ R[p] and the output space Y ⊂ R. We consider a common
nonparametric model

_Y_ _[t]_ = f _[∗](X_ _[t]) + ε[t], E(ε[t]) = 0,_ (1)

where f _[∗](·) is the ground truth function, and the innovation ε[t]_ is i.i.d. across time t ∈ Z. For the
sake of simplicity, we denote ρ[t] and ρ[t] [as the jointed distribution of][ (][X] _[t][, Y][ t][)][ and the corresponding]_
_X_
marginal distribution with respect to X _[t], respectively. This setup actually covers a large number of_
scenarios commonly used in practice. For instance, the case X _[t]_ that contains p lagged values of Y _[t]_

(e.g., X _[t]_ = Y [(][t][−][p][)] _× · · · × Y_ [(][t][−][1)]) corresponds to the p-order autoregressive models. Moreover, this
case can be viewed as a vector autoregressive model in the sense that input X _[t]_ includes the historical
information of multiple variables.

Although such a nonparametric model (1) makes very few assumptions on data generation, the related
nonparametric algorithms suffer so-called “curse of dimensionality”, see Fan & Gijbels (1996) for
further discussion. An effective strategy for solving this problem is additive model. Usually, the
additive structure is obtained by decomposing the input spaceUnder the assumption that the ground truth admits an additive structure X ∈ R f[p][∗]into= _Xj ==1_ _X[f][ ∗]j_ 1[, the additive] × ... × Xp.
model can be defined by

_Y_ _[t]_ = f1[∗][(][X]1[t][) +][ · · ·][ +][ f][ ∗]p [(][X]p[t][) +][ ε][t][,] [P][p] (2)

where each component fj[∗] [:][ X][j][ →] [R][ is a smooth function. In linear time series analysis, a weak]
stationarity condition (i.e., the first two moments of time series are time invariant) is preferred (Han
et al., 2015; Qiu et al., 2015). In contrast, strict stationarity is primarily used to analyse the nonlinear
time series (Fan & Yao, 2005).
**Definition 1. A stochastic process {Z** _[t]}t[∞]=−∞_ _[is strictly stationary if][ (][Z]_ [1][, ..., Z] _[t][)][ and]_
(Z [1+][k], ..., Z _[t][+][k]) have the same joint distributions for any t ∈_ Z and k ∈ Z.

Note that, if not otherwise stated, the stationarity in this paper refers to strict stationarity. Suppose
that we are given T size time series data (x[t], y[t]) _t=1_
data-generating model (2). Under stationarity condition and zero-mean Gaussian innovation with { _}[T]_ _[∈Z]_ _[T][ which are drawn from an additive]_
finite variance, the widely used methods that learn the ground truth f _[∗]_ usually integrate squared loss
and a smoothness- or sparsity-inducing regularizer Ω(·) into a structural risk minimization scheme:


(y[t] _−_
_t=1_

X


_fj(x[t]j[))][2][ + Ω(][f]_ [)][,]
_j=1_

X


min
_f_ _∈H_


-----

wheresubspace H := Hj { could be reproducing kernel Hilbert space (Chen et al., 2017; Kandasamy & Yu, 2016;f1 + ... + fp : fj ∈Hj, j = 1, ..., p} is an additive hypothesis space. Commonly, each
Raskutti et al., 2012), the orthogonal basis inducing space (Ravikumar et al., 2009; Meier et al., 2009;
Yang et al., 2018), or the composite function space with the neural network as a typical (Agarwal
et al., 2020; Bussmann et al., 2020).

However, when facing heavy-tailed innovation, these methods may have degraded performance due
to the amplification of the squared loss to large residuals. As one commonly used statistic in robust
learning community (Peng et al., 2019; Wang et al., 2017b; Feng & Wu, 2020), Huber loss is defined
as

(f (xt) _yt)2,_ if _f_ (x[t]) _y[t]_ _< σ_
_ℓσ(f_ (x[t]) − _y[t]) =_ 2σ _f_ (x −[t]) _y[t]_ _σ[2],_ if |f (x[t]) − _y[t]|_ _σ,_ (3)
 _|_ _−_ _| −_ _|_ _−_ _| ≥_

where σ is a positive hyper-parameter. Note that in the previous studies (Huber & Ronchetti, 2009;
Loh, 2017), the hyper-parameter σ is set to be fixed according to the 95% asymptotic efficiency rule.
However, Huber regression with a fixed scale parameter may not be able to learn the ground truth
when the noise is asymmetric, as argued recently in Feng & Wu (2020); Sun et al. (2019). In this
paper, we choose the scale parameter σ by relating it to the moment condition of the noise distribution
and the sample size so that the resulting regression estimator can asymptotically converge to the
ground truth function.

3 METHOD

3.1 SPARSE HUBER ADDITIVE MODELS


In this paper, we choose reproducing kernel Hilbert space (RKHS) _Kj_ _, j = 1, .., p, to form the_
_H_
the additive hypothesis space H, where each HKj is associated with a symmetric and positive
semi-definite Mercer kernel Kj : Xj × Xj → R. An additive RKHS is defined as

_HK = {f1 + ... + fp : fj ∈HKj_ _, j = 1, ..., p}_ (4)

with kernel norm ∥f _∥K[2]_ [= inf][{∥][f] _[∥]K[2]_ 1 [+][ ...][ +][ ∥][f] _[∥]K[2]_ _p_ _[}][. By integrating the Huber loss (3), additive]_
RKHS and kernel norm inducing regularizer into a Tikhonov regularization scheme, the regularized
Huber additive model with kernel-norm can be formulated as


_ℓσ(y[t]_
_−_
_t=1_

X


_fj(x[t]j[)) +][ η]_
_j=1_

X


_τj∥fj∥K[2]_ _j_ _[}][,]_ (5)
_j=1_

X


_fˆη =_ _j=1_ _fˆη,j =_ _f_ =[P][p]jarg min=1 _[f][j]_ _[,f][j]_ _[∈H][K]j_ _{t=1_ _ℓσ(y[t]_ _−_ _j=1_ _fj(x[t]j[)) +][ η]_ _j=1_ _τj∥fj∥K[2]_ _j_ _[}][,]_ (5)

X X X X

where η is positive regularization parameter and τj is the weight for j-th kernel norm. The representer
theorem (Wahba, 1990) ensures that _f[ˆ]η can be represented as_


_fˆη =_


_αtj[η]_ _[K][j][(][x]j[t]_ _[,][ ·][)][, α]tj[η]_
_t=1_ _[∈]_ [R][.]

X


_fˆη =_


_j=1_


To offer the method with sparsity, we consider the following sparsity-inducing penalty


_αtjKj(x[t]j[,][ ·][)][}][.]_
_t=1_

X


Ω(f ) := inf{


_τj∥αj∥2 : f =_
_j=1_

X


_j=1_


Let HZ be an additive data dependent hypothesis space defined by


_αtjKj(x[t]j[,][ ·][) :][ α][tj]_

_j=1_ _t=1_ _[∈]_ [R][}][.]

X X


_HZ = {f =_

Then the SpHAM can be formulated as

_fˆ = arg min_ [1]
_f_ _∈HZ_ _{_ _T_


_ℓσ(y[t]_
_−_
_t=1_

X


_fj(x[t]j[)) +][ λ][Ω(][f]_ [)][}][.]
_j=1_

X


-----

Denote K[t]j [= (][K][j][(][x]j[1][, x][t]j[)][, ..., K][j][(][x][T]j _[, x]j[t]_ [))][′][ ∈] [R][T][,][ α][ = (][α]1[′] _[, ..., α]p[′]_ [)][′][ ∈] [R][T p][ and][ α][j] [=]
(α1j, ..., αT j)[′] _∈_ R[T], where αj[′] [here refers to the transpose of][ α][j][ for avoiding the confusion.]
The SpHAM can be represented as


_αtj[λ]_ _[K][j][(][x][t]j[,][ ·][)]_ (6)
_t=1_

X


_fˆ =_


_j=1_


with


_α[λ]_ = arg min [1]
_αj_ _∈R[T]_ _,j=1,...,p{_ _T_


_ℓσ(y[t]_
_−_
_t=1_

X


(K[t]j[)][′][α][j][) +][ λ]
_j=1_

X


_τj_ _αj_ 2 _._ (7)
_∥_ _∥_ _}_
_j=1_

X


The optimization problem (7) can be solved by Fast Iterative Shrinkage-Thresholding Algorithm
(FISTA) (Beck & Teboulle, 2009). We provide the detailed optimization procedure in Appendix E.
**Remark 1. Inspired by the studies on group additive models (Yin et al., 2012; Pan & Zhu, 2017;**
_Chen et al., 2017), the modeling capacity of our proposed method can be improved by embedding the_
_(Hierarchical) group structure. For instance, group SpHAM can be easily formulated by replacing_
_the direct decomposition {Xj}j[p]=1[,][ X][j][ ∈]_ [R][ with the subgroups decomposition][ {X][d][}][D]d=1[, where each]
_d_ _may cover more than one variable._
_X_ _⊂X_

3.2 ASYMPTOTIC THEORY ANALYSIS

In this section, we provide the theoretical analyses of SpHAM, including the following:

-  The function approximation errors of SpHAM for stationary time series data (Theorem 1;
Section 3.2.1) and non-stationary time series data (Theorem 2; Section 3.2.2);

-  An adaptive SpHAM (inspired by above theoretical findings) and its function approximation
error for non-stationary time series (See Theorem 3; Section 3.2.3)

Due to space limitation, the high-level outline and detailed proofs are provided in the Appendix A-D.

3.2.1 FUNCTION APPROXIMATION ANALYSIS FOR STATIONARY TIME SERIES

Through this paper, the marginal distribution w.r.t Y _[t]_ is assumed to be almost everywhere supported
on [−M, M ] for some M ≥ 0.
**Assumption 1. Assume that |Y** _[t]|, ∀t ∈_ Z, is bounded and there exists a constant c > 0 such that
E|Y _[t]|[1+][c]_ _< ∞, ∀t ∈_ Z.

The moment condition in Assumption 1 is rather weak in the sense that the response variable Y _[t]_
possesses infinite variance. The same condition also applies to the distributions of the innovation ε[t],
implying that heavy-tailed innovation is allowed.


**Assumption 2. Let κ = supx∈X**


_Kj(x, x) <_ _,_ _j = 1, ..., p._
_∞_ _∀_


Assumption 2 only requires that the kernel is bounded under compact X, which holds for all Mercer
kernels (e.g., Gaussian kernel) and has been used in many learning theory literatures; see, e.g., (Wu
et al., 2006; Steinwart & Christmann, 2008; Wu & Zhou, 2008).

The following definitions are needed for Assumption 3. For any j = 1, ..., p, we define a kernel
integral operator LKj _,T +1 : L2(ρ[T][ +1]j_ [)][ →] _[L][2][(][ρ][T][ +1]j_ [)][ associated with the kernel][ K][j][ by]
_X_ _X_

_LKj_ _,T +1(f_ )(x[T]j [ +1]) = _Kj(x[T]j_ [ +1], uj)f (uj)dρ[T][ +1]j [(][u][j][)][.]
ZXj _X_

Note that LKj _,T +1 is a compact and positive operator on L2(ρ[T][ +1]j_ [)][. According to Mercer theorem, we]
_X_

can find the corresponding normalized eigenpairs {(ζi[j][, ψ]i[j][)][}][i][≥][1][ such that][ {][ψ]i[j][}][i][≥][1][ is an orthonormal]
basis ofby _L2(ρ[T]X[ +1]) and ζi[j]_ _[→]_ [0][ as][ i][ →∞][. Then for given][ r >][ 0][, we defined the][ r][-th power][ L]K[r] _j_ _,T +1_
_L[r]Kj_ _,T +1[(]_ _βi[j][ψ]i[j][) =]_ _βi[j][(][ζ]i[j][)][r][ψ]i[j][.]_
Xi≥1 Xi≥1


-----

**Assumption 3. We assume that fj[∗]** [:][ X][j][ →] [R][,][ ∀][j][ = 1][, ..., p][ is a function of the form][ f][ ∗]j [=]
_L[r]Kj_ _,T +1[(][g]j[∗][)][,][ ∀][r][ ∈]_ [(0][,][ 1]2 []][ with some][ g]j[∗] _j_ [)][,][ ∀][h][ ∈] [Z][.]

_[∈]_ _[L][2][(][ρ]X[T][ +1]_

Assumption 3 is a natural extension from i.i.d setting (Assumption 1 in (Wu et al., 2006; Christmann
& Zhou, 2016; Chen & Wang, 2018)) to non-i.i.d time series. Indeed, this assumption stands in
most practical cases. For instance, if r = 0.5, the ground truth function f _[∗]_ needs to be a real-valued
function in the RKHS _K. The RKHS admits a large class of bounded and continuous real-valued_
_H_
functions that are in Hilbert space. Moreover, this assumption has been widely used in learning
theory; please refer to (Smale & Zhou, 2003; Cucker & Zhou, 2007) for more discussions.
**Theorem 1. Suppose that the process {Z** _[t]}t[∞]=−∞_ _[is stationary. Let Assumptions 1-3 be true. By]_
_taking σ = T_ _[m], η = T_ _[−]_ 4[1]r and λ = T _[−]_ 4[1]r _[−][m], we have for any 0 < δ < 1,_

_f_ _f_ _C log(1/δ)T_ [Ψ(][m,c,r][)]
_∥_ [ˆ] − _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ [e]

_with confidence at least 1 −_ _δ, where_ _C is a positive constant independently of T, λ, η, δ and σ and_

max 12 _[, m][ −]_ [1][,][ −][cm][ +][ 1]4r [+][ m][ −] [1][}][,] _if m_ 1 4r
Ψ(m, c, r) = _{−_ [e] _≤_ _−_ [1]
max{− 2[1] _[, m][ −]_ [1][,][ −][cm][ +][ 1]2r [+ 2][m][ −] [2][}][,] _if m > 1 −_ 4[1]r _[.]_

**Remark 2. In the stationary case, our bound appears to be novel for the following reasons: a) the**
_result is completely independent of the various mixing conditions which are widely used in the theory_
_analysis of non-i.i.d dependent time series (Mohri & Rostamizadeh, 2009; Yang et al., 2018; Mohri_
_& Rostamizadeh, 2010; Guo & Shi, 2011; Qiu et al., 2015); b) compared with Yang et al. (2018), the_
_innovation assumption is rather weak in the sense that the innovation possesses infinite variance and_
_thus admits a heavy-tailed distribution._
**Corollary 1. Suppose that the process {Z** _[t]}t[∞]=−∞_ _[is stationary. Let all the conditions in Theorem 1]_
_be true. We then have for any 0 < δ < 1_

_f_ _f_ _C log(1/δ)T_ [Ψ(][m,c][)]
_∥_ [ˆ] − _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ [e]

_with confidence at least 1 −_ _δ, where_

Ψ(m, c) = max{− 12 _[,][ −][cm][ +][ m][ −]_ [1]2 _[}][,]_ _if m ≤_ [1]2

max _m_ 1, _cm + +2m_ 1 _,_ _if m >_ [1]2 _[.]_

 _{_ _−_ _−_ _−_ _}_


Figure 1 summaries the convergence
rates in Corollary 1 by taking different σ and c. If the innovation is Gaussian distribution with finite variance
(i.e., Assumption 1 holds for any
_c > 1), one can arbitrarily select a_
_σ = T_ _[m]_ (0 < m < [1]2 [)][ to obtain con-]

vergence rates O(T _[−]_ 2[1] ). Moreover,

we can see that the convergence rate
will decrease as m increases. Combined with the conclusion in Lemma
1 (i.e., the equivalence relation between Huber loss based empirical
risk and MSE as σ →∞), it indicates that σ indeed plays a tradeoff role between algorithmic robustness and variance-reduction. For the
weak moment condition (e.g., 0 < Figure 1: The convergence rates under different σ and c.
_c < 1), one may get slower conver-_
gence rates (e.g., O(T [(1][−][c][)][m][−] 2[1] ) or

_O(T_ [(2][−][c][)][m][−][1])), which also coincides with our intuitive understanding that small σ may be conducive
to robust forecasting. Note that our method will not converge when σ and c are both located in the
white area in Figure 1.


-----

As a comparison, due to the non-robustness of the squared loss, the most existing convergence
rates are established under the assumption that the innovation is Gaussian distribution with finite
variance (see, e.g., Yang et al. (2018); Han et al. (2015); Wang & Yang (2007); Kock & Callot
(2015)). However, from Theorem 1 and Corollary 1, the asymptotic convergence of SpHAM can
be obtained under weaker moment condition, which verifies the robustness of our method. Recall
the learning rate O(T 2−d2+1d ) derived in Yang et al. (2018), where d is the order of smoothness of

the component function fj, j = 1, ..., p. Relatively slow convergence rate O(T 2 ) we obtained

_[−]_ [1]

indicates the sacrifice for the absence of mixing condition.

3.2.2 FUNCTION APPROXIMATION ANALYSIS FOR NON-STATIONARY TIME SERIES

In nonstationary time series setting, different Z _[t]s may follow different distributions. Thus c. For_
given weights {st}t[T]=1[, we define the estimator][ ˆ]f **[s]** = _j=1_ _f[ˆ]j[s]_ [as the minimizer of the following]
weighted objective

_T_ _p_ [P][p]

_λ[( ˆ]f_ **[s]) := min** _stℓσ(y[t]_ _fj(x[t]j[)) +][ λ][Ω(][f]_ [)][}][.] (8)
_E_ **[s]** _f_ _Z[{]_ _−_
_∈H_ _t=1_ _j=1_

X X

Next, we introduce a discrepancy measure to characterize the discrepancy between the target distribution and the distributions of observations (Kuznetsov & Mohri, 2020a).

**Definition 2. For any f ∈HZ, the discrepancy measure with respect to Huber loss is defined as**


Eℓσ(f (x[T][ +1]) − _y[T][ +1]) −_


_stEℓσ(f_ (x[t]) − _y[t])_
_t=1_

X


disc(s) := sup
_f_ _∈HZ_


The discrepancy measures the non-stationarity of the stochastic process {Z _[t]}t[∞]=−∞_ [with respect to]
both the loss function ℓσ and the hypothesis set HZ.

**Theorem 2.σ = T** 21c, λ = Let Assumptions 1-3 be true. We assume that T _[−][1]_ _and η = T_ _[−]_ 2[1], we have for any 0 < δ < 1 fj[∗] _[∈H][K]j_ _[,][ ∀][j][ = 1][, ..., p][. By taking]_

1
_∥f[ˆ][s]_ _−_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ [disc(][s][) +][ E]λ[s][( ˆ]f **[s]) +** _C1∥s∥2T_ 2c + _C2 log(1/δ)T_ _[−]_ [1]2

_with confidence at least 1 −_ _δ, where_ _C1,_ _C2 are two positive constants independently of[e]_ e _T, λ, η, δ_
_and σ._

[e] [e]

Note that there are several existing studies towards analyzing nonstationary and non-mixing time
series, see, e.g., Kuznetsov & Mohri (2020a). Different from them, the error bound we derived is
with respect to the function approximation rather than generalization risk, which is very crucial for
Huber regression problem, since the convergence of generalization risk cannot imply the convergence
of function approximation (Sun et al., 2019; Feng & Wu, 2020).

3.2.3 ADAPTIVE SPARSE HUBER ADDITIVE MODEL FOR NONSTATIONARY TIME SERIES

Theorem 2 illustrates that we shall minimize the following optimization problem for non-stationary
time series forecasting:


_fj(x[t]j[)) + disc(][s][) +][ λ][1][Ω(][f]_ [) +][ λ][2][T]
_j=1_

X


1

2c **s** 2 _,_
_∥_ _∥_ _}_


min _stℓσ(yt_
_f_ _∈HZ[{]t=1_ _−_

X


where λ1 and λ2 are two positive regularization parameters, and c is a positive constant introduced
in Assumption 1. Although the discrepancy measure disc(s) is crucial for such an optimization
problem, we cannot obtain its exact value since we do not have access to the distributions of Z _[t], t ∈_ Z.
Hence, we need to estimate the approximated discrepancy from given data. Inspired by Kuznetsov &
Mohri (2020a), one natural and necessary assumption is that there exists an underlying representation
relationship between distribution ρ[T][ +1] and distributions ρ[t], t = 1, ..., T .


-----

**Assumption 4.following term is sufficiently small: Denote by a probability set q[∗]** = {qt[∗][}]t[T]=1 _[with][ P]t[T]=1_ _[q]t[∗]_ [= 1][. We assume that the]

_T_

disc(q[∗]) := supf _∈HZ[Eℓσ(f_ (x[T][ +1]) − _y[T][ +1]) −_ _t=1_ _qtEℓσ(f_ (x[t]) − _y[t])]._

X

Note that the priori q[∗] can be any distribution, which shall be given empirically according to the
trend of the time series. For instance, in a particular scenario of the distribution ρ[T][ +1] does not change
drastically compared with the distributions ρ[t], t = 1, ..., T, Kuznetsov & Mohri (2020a) have proven
that the Assumption 4 holds if the s[∗] is an uniform distribution over last l > 0 observations, where l
is a hyper-parameter that can be tuned in practical applications.
**Theorem 3. Let Assumption 4 and the conditions in Theorem 2 be true. We have for any 0 < δ < 1**


_∥f[ˆ][s]_ _−_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ _fsup∈HZ_ _t=1(qt[∗]_ _[−]_ _[s][t][)][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][) +][ E]λ[s][( ˆ]f_ **[s])**

X


1
+C1 log(1/δ)(∥q[∗] _−_ **s∥2 + ∥s∥2T** 2c ) + _C2 log(1/δ)T_ _[−]_ 2[1]

_with confidence at least 1 −_ _δ, where[e]C1,_ _C2 are two constants independently of e_ _T, λ, η, δ and σ._

Finally, the optimization problem of adaptive SpHAM can be formulated as following two stages:

[e] [e]

**Step A: finding the weight ˆs:**


_t=1(qt[∗]_ _[−]_ _[s][t][)][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][) +][ λ][1][∥][q][∗]_ _[−]_ **[s][∥]2[2]** [+][ λ][2][T]

X


1

2c ∥s∥2[2][}] (9)


ˆs = arg min
**s** _[{]f[ sup]∈HZ_

**Step B: forecasting:**

where


_fˆ[s](x[T][ +1]) =_


_αtj[s]_ _[K][(][x]j[t]_ _[, x][T]j_ [ +1]),
_t=1_

X


_j=1_


_α[s]_ = arg min
_αj_ _∈R[T]_ _,j=1,...,p_


_sˆtℓσ(y[t]_
_−_
_t=1_

X


(K[t]j[)][′][α][j][) +][ λ]
_j=1_

X


_τj_ _αj_ 2 _._ (10)
_∥_ _∥_ _}_
_j=1_

X


The optimization problems (9) can be solved by standard gradient descent method. Similarly to
the strategy for optimization problem (7), we use Fast Iterative Shrinkage-Thresholding Algorithm
(FISTA) (Beck & Teboulle, 2009) for Step B. We provide the detailed procedure in Appendix E.

4 EXPERIMENT

This section validates the effectiveness of SpHAM and adaptive SpHAM. In all experiments, the
Gaussian kernel Kj(u, v) = exp( 2d[2] 2 ), where j = 1, ..., p and bandwidth d > 0, is employed
_−_ _[∥][u][−][v][∥][2]_

for constructing the additive data dependent hypothesis space. Due to limited space, the evaluations
on real-world data are provided in Appendix F. We consider two synthetic examples as below:

_Example A: Inspired by (Kuznetsov & Mohri, 2020a), a stationary time series is generated according_
to the non-linear additive autoregressive model:

_Y_ _[t]_ = [3]

2 [sin(] _[π]2_ _[Y][ t][−][2][)][ −]_ [sin(] _[π]2_ _[Y][ t][−][3][) + 2][ε][t][,]_


where the innovation ε[t]s are i.i.d. drawn from Gaussian distribution N (0, 0.5) and Student distribution
with degree of freedom 2, respectively.

_Example B: Inspired by (Kuznetsov & Mohri, 2020a), a time series with smooth drift is generated by_


_t_

400 [sin(][Y][ t][−][1][) + 1]2 _[ε][t][,]_


_Y_ _[t]_ =


-----

where the distributions of noise ε[t] are the same as above.

_Hyper-parameter selection and evaluation criterions: Recall that SpHAM algorithm requires three_
hyper-parameters: regularization parameter λ, bandwidth of kernel d and Huber parameter σ. We set
these parameters according to the suggestions in our Theorems. Based on the suggestion in Theorem
1
1-3, the selection of Huber parameter σ is σ = T 48 and the regularization parameter is λ = T _[−][1]._

Moreover, we set the bandwidth d = 0.5 and tune l ∈{100, 150, 200, ..., 350}. The evaluation

criterions for forecasting used here contains Average Sample Error(ASE)= _N[1]_ _Nt=1[( ˆ]f_ (x[t]) _y[t])[2]_

_−_

and True Deviation (TD)= _N[1]_ _Nt=1[( ˆ]f_ (x[t]) − _f_ _[∗](x[t]))[2], where N is the number of test samples.qP_

For each example, we generate time series withqP 4000 sample points. The samples at time t =
_{1500, 1501, ..., 1899} are used as a training set, and the samples at next time t = {1900, ..., 1999}_
are considered as the test data. The competitors include Simple Exponential Smoothing (SES),
TS-SpAM (Yang et al., 2018) and Vanilla Long Short-Term Memory (LSTM). The tuning parameters
of other competing methods such as LSTM and SES, are chosen according to their original python
packages. All the evaluations are repeated for 50 times. The average results for Examples A-B are
presented in Table 2.

From Table 2, the results on Example A verify that our SpHAM is competitive with other approaches
based on square loss under Gaussian noise, and performs better in presence of heavy-tailed t noise.
Moreover, the results on Example B shows the promising performance of our adaptive SpHAM for
non-stationary time series forecasting.


Table 2: The results on synthetic data.

Gaussian noise Student noise

Methods ASE (Std) TD (Std) ASE (Std) TD (Std)

SES 0.114(±.008) 0.098(±.005) 0.493(±.313) 0.105(±.009)
LSTM **0.063(±.005)** **0.022(±.005)** 0.538(±.342) 0.219(±.089)
TS-SpAM 0.070(±.007) 0.035(±.008) 0.505(±.332) 0.131(±.042)
SpHAM (ours) 0.076(±.005) 0.042(±.008) **0.492(±.333)** **0.095(±.006)**

SES 0.435(±.015) 0.424(±.015) 0.402(±.080) 0.366(±.017)
LSTM 0.120(±.125) 0.114(±.127) 0.247(±.125) 0.181(±.088)
TS-SpAM 0.114(±.022) 0.105(±.021) 0.173(±.125) 0.104(±.051)
SpHAM (ours) 0.115(±.028) 0.112(±.029) 0.179(±.116) 0.096(±.026)
Adaptive SpHAM (ours) **0.102(±.026)** **0.096(±.026)** **0.172(±.118)** **0.085(±.024)**


Example A

Example B


5 CONCLUSION

We propose an adaptive sparse Huber additive model by integrating Huber loss and ℓ2,1-norm
regularizer into an additive data dependent hypothesis space. We theoretically explore the asymptotic
properties of our method for both non-Gaussian and (non)stationary time series. Experimental results
on both synthetic and real-world data validate the effectiveness of the proposed method.

ACKNOWLEDGMENTS

This work is supported by the Major Science and Technology Innovation 2030 “New Generation Artificial Intelligence” key project (No. 2021ZD0111700) and the National Natural Science Foundation
of China under Grant No. 12071166. We sincerely appreciate the anonymous ICLR reviewers for
their helpful comments.

REFERENCES

Terrence M. Adams and Andrew B. Nobel. Uniform convergence of vapnik-chervonenkis classes
under ergodic sampling. The Annals of Probability, 38(4):1345–1367, 2010.

Alekh Agarwal and John C. Duchi. The generalization ability of online algorithms for dependent
data. IEEE Transactions on Information Theory, 59(1):573–587, 2013.


-----

Rishabh Agarwal, Nicholas Frosst, Xuezhou Zhang, Rich Caruana, and Geoffrey E. Hinton. Neural
additive models: Interpretable machine learning with neural nets. arXiv:2004.13912v1, 2020.

Richard T. Baillie. Long memory processes and fractional integration in econometrics. Journal of
_Econometrics, 73(1):5–59, 1996._

Amir Beck and Marc Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse
problems. SIAM Journal on Imaging Sciences, 2(1):183–202, 2009.

Bart Bussmann, Jannes Nys, and Steven Latre. Neural additive vector autoregression models for
causal discovery in time series data. arXiv:2010.09429v1, 2020.

John M. Chambers and Trevor Hastie. Statistical Models in S. Wadsworth & Brooks/Cole Advanced
Books & Software, 1992.

Hong Chen and Yulong Wang. Kernel-based sparse regression with the correntropy-induced loss.
_Applied and Computational Harmonic Analysis, 44(1):144–164, 2018._

Hong Chen, Xiaoqian Wang, Cheng Deng, and Heng Huang. Group sparse additive machine. In
_Advances in Neural Information Processing Systems (NIPS), pp. 198–208. 2017._

Hong Chen, Changying Guo, Yingjie Wang, and Huijuan Xiong. Sparse additive machine with ramp
loss. Analysis and Applications, 19(3):509 – 528, 2021a.

Hong Chen, Yingjie Wang, Feng Zheng, Cheng Deng, and Heng Huang. Sparse modal additive
model. IEEE Transactions on Neural Networks and Learning Systems, 32(6):2373 – 2387, 2021b.

Andreas Christmann and Ding-Xuan Zhou. Learning rates for the risk of kernel based quantile
regression estimators in additive models. Analysis and Applications, 14(3):449–477, 2016.

Tianjiao Chu and Clark Glymour. Search for additive nonlinear time series causal models. Journal of
_Machine Learning Research, 9:967–991, 2008._

Benoit Colson, Patrice Marcotte, and Gilles Savard. An overview of bilevel optimization. Annals of
_Operations Research, 153(1):235–256, 2007._

Felipe Cucker and Ding-Xuan Zhou. _Learning Theory: An Approximation Theory Viewpoint._
Cambridge University Press, 2007.

Francesca Dominici, Aidan McDermott, Scott L. Zeger, and Jonathan M. Samet. On the use of
generalized additive models in time-series studies of air pollution and health. American Journal of
_Epidemiology, 156(3):193–203, 2002._

Paul Doukhan. Mixing: Properties and Examples. Lecture Notes in Statistics. Springer, 1994.

Jianqing Fan and Irene Gijbels. Local Polynomial Modelling and Its Applications. Chapman and
Hall, 1996.

Jianqing Fan and Qiwei Yao. Nonlinear time series: Nonparametric and parametric methods.
Springer, 2005.

Yunlong Feng and Qiang Wu. A statistical learning assessment of huber regression.
_arXiv:2009.12755v1, 2020._

Zheng-Chu Guo and Lei Shi. Classification with non-i.i.d. sampling. Mathematical and Computer
_Modelling, 54(5):1347–1364, 2011._

Fang Han, Huanran Lu, and Han Liu. A direct estimation of high dimensional stationary vector
autoregressions. Journal of Machine Learning Research, 16:3115–3150, 2015.

Trevor J. Hastie and Robert J. Tibshirani. Generalized additive models. London: Chapman and Hall,
1990.

Fengxiang He and Dacheng Tao. Recent advances in deep learning theory. _arXiv preprint_
_arXiv:2012.10931, 2020._


-----

Jianhua Z. Huang and Lijian Yang. Identification of non-linear additive autoregressive models.
_Journal of the Royal Statistical Society. Series B (Statistical Methodology), 66(2):463–477, 2004._

Peter J. Huber. Robust estimation of a location parameter. The Annals of Mathematical Statistics, 35
(1):73–101, 1964.

Peter J. Huber and Elvezio M. Ronchetti. Robust Statistics. Wiley, 2009.

Egide Kalisa, Sulaiman Fadlallah, Mabano Amani, Lamek Nahayo, and Gabriel Habiyaremye.
Temperature and air pollution relationship during heatwaves in birmingham, uk. Sustainable Cities
_and Society, 43:111–120, 2018._

Kirthevasan Kandasamy and Yaoliang Yu. Additive approximations in high dimensional nonparametric regression via the SALSA. In International Conference on Machine Learning (ICML), pp.
69–78, 2016.

Anders Bredahl Kock and Laurent Callot. Oracle inequalities for high dimensional vector autoregressions. Journal of Econometrics, 186(2):325–344, 2015.

Vitaly Kuznetsov and Mehryar Mohri. Generalization bounds for time series prediction with nonstationary processes. In International Conference on Algorithmic Learning Theory, 2014.

Vitaly Kuznetsov and Mehryar Mohri. Discrepancy-based theory and algorithms for forecasting
non-stationary time series. Annals of Mathematics and Artificial Intelligence, 88:367–399, 2020a.

Vitaly Kuznetsov and Mehryar Mohri. Learning theory and algorithms for forecasting non-stationary
time series. In Advances in Neural Information Processing Systems (NIPS), 2020b.

Guodong Liu, Hong Chen, and Heng Huang. Sparse shrunk additive models. In International
_Conference on Machine Learning (ICML), 2020._

Jiajia Liu, Yudong Ye, Chenglong Shen, Yuming Wang, and R Erdelyi. A new tool for cme arrival
time prediction using machine learning algorithms: Cat-puma. The Astrophysical Journal, 855(2):
109–118, 2018.

Po-Ling Loh. Statistical consistency and asymptotic normality for high-dimensional robust mestimators. The Annals of Statistics, 45(2):866–896, 2017.

Lukas Meier, Sara Van De Geer, and Peter Buhlmann. High-dimensional additive modeling. The
_Annals of Statistics, 37(6B):3779–3821, 2009._

Mehryar Mohri and Afshin Rostamizadeh. Rademacher complexity bounds for non-i.i.d. processes.
In Advances in Neural Information Processing Systems (NIPS), 2009.

Mehryar Mohri and Afshin Rostamizadeh. Stability bounds for stationary φ-mixing and β-mixing
processes. In Advances in Neural Information Processing Systems (NIPS), 2010.

George Ruchathi Mwaniki, Chelsea Rosenkrance, H. Will Wallace, B. Tom Jobson, Mathew H.
Erickson, Brian K. Lamb, Rick J. Hardy, Rasa Zalakeviciute, and Timothy M. VanReken. Factors
contributing to elevated concentrations of pm2.5 during wintertime near boise, idaho. Atmospheric
_Pollution Research, 5(1):96–103, 2014._

Chao Pan and Michael Zhu. Group additive structure identification for kernel nonparametric regression. In Advances in Neural Information Processing Systems (NIPS), pp. 4907–4916. 2017.

Jiangtao Peng, Luoqing Li, and Yuan Yan Tang. Maximum likelihood estimation-based joint sparse
representation for the classification of hyperspectral remote sensing images. IEEE Transactions on
_Neural Networks and Learning Systems, 30(6):1790–1802, 2019._

Huitong Qiu, Sheng Xu, Fang Han, Han Liu, and Brian Caffo. Robust estimation of transition
matrices in high dimensional heavy-tailed vector autoregressive processes. In Proceedings of the
_32nd International Conference on Machine Learning, volume 37, pp. 1843–1851, 2015._


-----

Alexander Rakhlin, Karthik Sridharan, and Ambuj Tewari. Online learning: random averages,
combinatorial parameters, and learnability. In Advances in Neural Information Processing Systems
_(NIPS), pp. 1984–1992, 2010._

Alexander Rakhlin, Karthik Sridharan, and Ambuj Tewari. Sequential complexities and uniform
martingale laws of large numbers. Probability Theory and Related Fields, 161(1):111–153, 2015.

Garvesh Raskutti, Martin J. Wainwright, and B Yu. Minimax-optimal rates for sparse additive models
over kernel classes via convex programming. Journal of Machine Learning Research, 13(2):
389–427, 2012.

Pradeep Ravikumar, Han Liu, John Lafferty, and Larry Wasserman. SpAM: sparse additive models.
_Journal of the Royal Statistical Society: Series B, 71:1009–1030, 2009._

Khaiwal Ravindra, Preety Rattan, Suman Mor, and Ashutosh Nath Aggarwal. Generalized additive
models: Building evidence of air pollution, climate change and human health. Environment
_International, 132:104987, 2019._

Meir Ron. Nonparametric time series prediction through adaptive model selection. Machine Learning,
39(1):5–34, 2000.

Jakob Runge, Sebastian Bathiany, Erik M. Bollt, Gustau Camps-Valls, Dim Coumou, Ethan R Deyle,
Clark Glymour, Marlene Kretschmer, Miguel D. Mahecha, Jordi Munoz-Mar˜ ´ı, Egbert H. van
Nes, Jonas Peters, Rick Quax, Markus Reichstein, Marten Scheffer, Bernhard Scholkopf, Peter L.¨
Spirtes, George Sugihara, Jie Sun, Kun Zhang, and Jakob Zscheischler. Inferring causation from
time series in earth system sciences. Nature Communications, 10, 2019.

Stephen Smale and Ding-Xuan Zhou. Estimating the approximation error in learning theory. Analysis
_and Applications, 01:17–41, 2003._

Qiongxia Song and Lijian Yang. Oracally efficient spline smoothing of nonlinear additive autoregression models with simultaneous confidence band. Journal of Multivariate Analysis, 101(9):
2008–2025, 2010.

Ingo Steinwart and Andreas Christmann. Support Vector Machines. Springer Science and Business
Media, 2008.

Charles J. Stone. Additive regression and other nonparametric models. The Annals of Statistics, 13
(2):689–705, 1985.

Qiang Sun, Wen-Xin Zhou, and Jianqing Fan. Adaptive huber regression. Journal of the American
_Statistical Association, 115(529):254–265, 2019._

Pham Dinh Tao and Le Thi Hoai An. A d.c. optimization algorithm for solving the trust-region
subproblem. SIAM Journal on Optimization, 8:476–505, 1998.

Grace Wahba. Spline models for observational data. Society for Industrial and Applied Mathematics,
1990.

Li Wang and Lijian Yang. Spline-backfitted kernel smoothing of nonlinear additive autoregression
model. The Annals of Statistics, 35(6):2474 – 2503, 2007.

Xiaofeng Wang and Donald E. Brown. The spatio-temporal generalized additive model for criminal
incidents. In Proceedings of 2011 IEEE International Conference on Intelligence and Security
_Informatics, pp. 42–47, 2011._

Xiaoqian Wang, Hong Chen, Weidong Cai, Dinggang Shen, and Heng Huang. Regularized modal
regression with applications in cognitive impairment prediction. In Advances in Neural Information
_Processing Systems (NIPS), 2017a._

Yimin Wang, Jiajia Liu, Ye Jiang, and Robert Erdelyi. Cme arrival time prediction using convolutional´
neural network. The Astrophysical Journal, 881(11):15, 2019.


-----

Yingjie Wang, Hong Chen, Feng Zheng, Chen Xu, Tieliang Gong, and Yanhong Chen. Multi-task
additive models for robust estimation and automatic structure discovery. In Advances in Neural
_Information Processing Systems (NIPS), 2020._

Yulong Wang, Yuan Yan Tang, and Luoqing Li. Correntropy matching pursuit with application to
robust digit and face recognition. IEEE Transactions on Cybernetics, 47(6):1354–1366, 2017b.

Kam Chung Wong, Zifan Li, and Ambuj Tewari. Lasso guarantees for β-mixing heavy-tailed time
series. The Annals of Statistics, 48(2):1124–1142, 2020.

Qiang Wu and Ding Xuan Zhou. Learning with sample dependent hypothesis spaces. Computers and
_Mathematics with Applications, 56(11):2896–2907, 2008._

Qiang Wu, Yiming Ying, and Ding-Xuan Zhou. Learning rates of least-square regularized regression.
_Foundations of Computational Mathematics, 6(2):171–192, 2006._

Qiang Wu, Yiming Yin, and Ding-Xuan Zhou. Multi-kernel regularized classfiers. Biometrika, 23:
108–134, 2007.

Wenjun Yan, Marcus A. Worsley, Thang Pham, Alex Zettl, Carlo Carraro, and Roya Maboudian.
Effects of ambient humidity and temperature on the no2 sensing characteristics of ws2/graphene
aerogel. Applied Surface Science, 450:372–379, 2018.

Yingxiang Yang, Adams Wei Yu, Zhaoran Wang, and Tuo Zhao. Detecting nonlinear causality in
multivariate time series with sparse additive models. arXiv:1803.03919v2, 2018.

Junming Yin, Xi Chen, and Eric P Xing. Group sparse additive models. In International Conference
_on Machine Learning (ICML), pp. 871–878, 2012._

Bin Zou, Luoqing Li, and Zongben Xu. The generalization performance of erm algorithm with
strongly mixing observations. Machine Learning, 75(3):275–295, 2009.


-----

A PROOF SKETCH

This section sketches the proof focusing on the conceptual aspects. The full proofs are provided in
Appendix B-D.

|Col1|Col2|
|---|---|
|Lemma 1 Lemma 2||


|Col1|Col2|
|---|---|
|Proposition 2 Proposition 3 Proposition 4||
|Lemma 3 Lemma 4 Lemma 5 Lemma 6||


Theorem 1

(Approximation Error Bound)

**Proposition 5**


Figure 2: The important ingredients for Theorem 1.

Figure 2 summaries the important ingredients for the proof of Theorem 1. Lemma 1 describes the
properties off[ˆ] and Lemma 2 illustrates the relation between huber-based risk and MSE. Following
these two lemmas and previous works (Wang et al., 2017a; Feng & Wu, 2020), we then decompose
the function error into three important parts: sample error, hypothesis error and approximation error.
To bound the sample error, we need to employ the sequential Rademacher complexity (Kuznetsov &
Mohri, 2020a), which is the key to measuring the capacity of the data-dependent hypothesis space
for non-i.i.d data. Lemmas 3-4 give an important concentration inequality and the upper bound
of sequential Rademacher complexity, respectively. Finally, Theorem 1 is obtained by combining
Propositions 2-4. The proof of Theorems 2 is similar to the above process, except that we develop a
different error decomposition in Proposition 5. On the basis of Theorem 2, by applying an assumption
on data distribution, we derive Theorem 3.

The full proofs of Theorems 1-3 are provided in Sections C-E, respectively.

B PROOF OF THEOREM 1

We first illustrate two key properties of _f[ˆ]._

**Lemma 1. Let Assumptions 1-2 be true. From the definition of** _f[ˆ] in Eq. (6), there hold_

2κ[2]M [2]
_∥f[ˆ]∥∞_ _≤_ _M ˆf_ [:=] _λT minj=1,...,p τj_ _._


_and_


_p_ 2M [2]

_τj_ _αj[λ]_ _f_ [=]
_j=1_ _∥_ _[∥][2]_ _[≤]_ _[κ][−][2][M][ ˆ]_ _λT minj=1,...,p τj_

X


Ω( f[ˆ]) =


_where M is a positive constant such that |Y_ _[t]| ≤_ _M, ∀t ∈_ Z.

_Proof. Denote by ET[σ][(][f]_ [) :=][ P]t[T]=1 _[ℓ][σ][(][f]_ [(][x][(][t][)][)] _[−]_ _[y][(][t][)][)][ for notational convenience. From the definition]_
of _f[ˆ], we know that_
_ET[σ][( ˆ]f_ ) + λΩ( f[ˆ]) ≤ET[σ][(0)][.]
It means that


Ω( f[ˆ]) =


_j=1_ _τj∥αj[λ][∥][2]_ _[≤]_ _[λ][−][1][E]T[σ][(0)][.]_

X


-----

From the definition of Huber loss (3), we obtain


_T_ [(0) = 1]
_E_ _[σ]_ _T_


_ℓσ(y[(][t][)])_ _σ[2], if_ _y[(][t][)]_ _σ,_ _t = 1, ..., T,_
_≤_ _|_ _| ≤_ _∀_
_t=1_

X


and

_T_ [(0) = 1]
_E_ _[σ]_ _T_

Therefore,


_ℓσ(y[(][t][)]) = 2σ_ _y[(][t][)]_ _σ[2]_ 2σM _σ[2], if_ _y[(][t][)]_ _> σ,_ _t = 1, ..., T._
_|_ _| −_ _≤_ _−_ _|_ _|_ _∀_
_t=1_

X


_p_ 2M [2]

_j=1_ _τj∥αj[λ][∥][2]_ _[≤]_ _λT minj=1,...,p τj_

X

According to the property of RKHS, we conclude that


_p_ 2κ[2]M [2]

_j=1_ _∥αj[λ][∥][2]_ _[≤]_ _λT minj=1,...,p τj_

X


_∥f[ˆ]∥∞_ _≤_ _κ∥f[ˆ]∥K ≤_ _κ_


_fj_ _Kj_
_j=1_ _∥_ [ˆ] _∥[2]_ _[≤]_ _[κ][2]_

X


This completes the proof by denoting M ˆf [=] _λT min2κj[2]=1M,...,p[2]_ _τj_ [.]

Let _H[¯]Z be a data-dependent hypothesis space constrained by sparsity-inducing regularizer Ω(·), i.e.,_


_p_ 2M [2]

_τj_ _αj_ 2 _._
_j=1_ _∥_ _∥_ _≤_ _λT minj=1,..,p τj_ _}_

X


_αijKj(x[(]j[t][)][, x][j][) :]_
_i=1_

X


_H¯Z = {f_ (x) =


_j=1_


For any h ∈ Z, denote by ET[σ] +1[(][f] [) :=][ E][ℓ][σ][(][f] [(][x][T][ +1][)][−][y][T][ +1][)][ and][ E][T][ +1][(][f] [) :=][ E][(][f] [(][x][T][ +1][)][−][y][T][ +1][)][2][,]
respectively. In the following, we describe the relationship between ET[σ] +1[(][f] [)][ and][ E][T][ +1][(][f] [)][,][ ∀][f][ ∈] _H[¯]Z_ .
Note that a similar proof is given in (Feng & Wu, 2020) for the i.i.d case.
**Lemma 2. Let Assumptions 1-2 be true. For any h ∈** Z and f ∈ _H[¯]Z, there holds_

[ _T +1[(][f]_ [)][ −E]T[σ] +1[(][f][ ∗][)]][ −] [[][E][T][ +1][(][f] [)][ −E][T][ +1][(][f][ ∗][)]]
_E_ _[σ]_ _≤_ _[M]σ[c][c][,]_

_where both c and Mc = 2[3+][c](Mf ∗_ + M ˆf [+ 1)][2][E][(][|][Y][ T][ +1][|][1+][c][)][ are positive constants.]

_Proof. For any σ > max{M + M ˆf_ _[,][ 1][}][, we denote two events][ I][Y][ and][ II][Y][ as follows]_

_IY T +1 =_ _Y_ _[T][ +1]_ : _Y_ _[T][ +1]_ _σ/2_
_{_ _|_ _| ≥_ _}_

and
_IIY T +1 =_ _Y_ _[T][ +1]_ : _Y_ _[T][ +1]_ _< σ/2_ _._
_{_ _|_ _|_ _}_
Then, for any f ∈ _H[¯]Z and ∥f_ _[∗]∥∞_ = Mf ∗ _< ∞, we have_

[ET[σ] +1[(][f] [)][ −E]T[σ] +1[(][f][ ∗][)]][ −] [[][E]T +1[(][f] [)][ −E]T +1[(][f][ ∗][)]]

= h [Z]Z _ℓσ(f_ (x[T][ +1]) − _y[T][ +1])dρ[T][ +1]_ _−_ ZZ _ℓσ(f_ _[∗](x[T][ +1]) −_ _y[T][ +1])dρ[T][ +1][i]_ _−∥f −_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1])

= h [Z]X ZY _ℓσ(f_ (x[T][ +1]) − _y[T][ +1]) −_ _ℓσ(f_ _[∗](x[T][ +1]) −_ _y[T][ +1])dρ[T]Y|X[ +1][dρ][T]X[ +1]i_ _−∥f −_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1])

_≤_ ZX ZIY T +1 ∪IIY T +1 _ℓσ(f_ (x[T][ +1]) − _y[T][ +1]) −_ _ℓσ(f_ _[∗](x[T][ +1]) −_ _y[T][ +1])dρ[T]Y|X[ +1][dρ][T]X[ +1]_

+ ZX ZIY T +1 ∪IIY T +1 (f (x[T][ +1]) − _y[T][ +1])[2]_ _−_ (f _[∗](x[T][ +1]) −_ _y[T][ +1])[2]dρ[T]Y|X[ +1][dρ][T]X[ +1]_

_≤_ ZX ZIY T +1 _ℓσ(f_ (x[T][ +1]) − _y[T][ +1]) −_ _ℓσ(f_ _[∗](x[T][ +1]) −_ _y[T][ +1])dρ[T]Y|X[ +1][dρ][T]X[ +1]_

+ ZX ZIY T +1 (f (x[T][ +1]) − _y[T][ +1])[2]_ _−_ (f _[∗](x[T][ +1]) −_ _y[T][ +1])[2]dρ[T]Y|X[ +1][dρ][T]X[ +1]_ _._


-----

The last inequality is based on the fact that

_|y[T][ +1]_ _−_ _f_ (x[T][ +1])| < |y[T][ +1]| + ∥f _∥∞_ _< σ, ∀(x[T][ +1], y[T][ +1]) ∈X × IIY T +1_

and hence
_ℓσ(f_ (x[T][ +1]) _y[T][ +1]) = (f_ (x[T][ +1]) _y[T][ +1])[2]._
_−_ _−_
Recalling that the Huber loss is 2σ-Lipschitz continuous, then we get
ZX ZIY T +1 _ℓσ(f_ (x[T][ +1]) − _y[T][ +1]) −_ _ℓσ(f_ _[∗](x[T][ +1]) −_ _y[T][ +1])dρY|X[T][ +1][dρ][T]X[ +1]_


_≤_ 2σ ZX ZIY T +1 _|f_ (x) − _f_ _[∗](x)|dρY|X[T][ +1][dρ][T]X[ +1]_

2σ _f_ _f_ Prob _IY_ _._
_≤_ _∥_ _−_ _[∗]∥∞_ _{_ _}_

According to Markov’s inequality, for any constant c > 0, we further have

Prob{IY T +1 _} = Prob{|Y_ _[T][ +1]| ≥_ _σ/2} = Prob{|Y_ _[T][ +1]|[1+][c]_ _≥_ (σ/2)[1+][c]} ≤ [2][1+][c][E][(]σ[|][Y][1+][ T][c][ +1][|][1+][c][)]

Then we have


ZX ZIY T +1 _ℓσ(f_ (x[T][ +1]) − _y[T][ +1]) −_ _ℓσ(f_ _[∗](x[T][ +1]) −_ _y[T][ +1])dρY|X[T][ +1][dρ][T]X[ +1]_

_._

_≤_ [2][2+][c][∥][f][ −] _[f][ ∗][∥]σ[∞][c][E][(][|][Y][ T][ +1][|][1+][c][)]_

Moreover, the second term in the right-hand side can be bounded by
ZX ZIY T +1 (f (x[T][ +1]) − _y[T][ +1])[2]_ _−_ (f _[∗](x[T][ +1]) −_ _y[T][ +1])[2]dρY|X[T][ +1][dρ][T]X[ +1]_


_≤∥f −_ _f_ _[∗]∥∞_ ZX ZIY T +1 _|2y[T][ +1]_ _−_ _f_ (x[T][ +1]) − _f_ _[∗](x[T][ +1])|dρY|X[T][ +1][dρ]X[T][ +1]_

_≤∥f −_ _f_ _[∗]∥∞_ _IY T +1_ 2|y[T][ +1]|dρY[T][ +1] + (∥f _∥∞_ + ∥f _[∗]∥∞)Prob(IY T +1_ )
h [Z]

According to Holder inequality, we have


1

1+c
_≤_ [2][c][E][(][|][Y][ T]σ[c][ +1][|][1+][c][)]


_IY T +1_ _|Y_ _[T][ +1]|dρ[T]Y[ +1]_ _≤_ (Prob(IY T +1 ))


_c_

1+c (E|Y _[T][ +1]|[1+][c])_


Then we can deduce that

_∥f −_ _f_ _[∗]∥∞_ _IY T +1_ 2|y[T][ +1]|dρY[T][ +1] + (∥f _∥∞_ + ∥f _[∗]∥∞)Prob(IY T +1_ )
h [Z] i

+ [2][1+][c][(][∥][f] _[∥][∞]_ [+][ ∥][f][ ∗][∥][∞][)][2][E][(][|][Y][ T][ +1][|][1+][c][)]

_≤_ [(2][2+][c][ + 2][1+][c][)][∥][f][ −]σ[f][c][ ∗][∥][∞][E][(][|][Y][ T][ +1][|][1+][c][)] _σ[1+][c]_


Finally, there holds
h [Z]Z _ℓσ(f_ (x[T][ +1]) − _Y_ _[T][ +1]) −_ ZZ _ℓσ(f_ _[∗](x[T][ +1]) −_ _Y_ _[T][ +1])i_ _−∥f −_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]

with
_Mc = 2[3+][c](Mf ∗_ + M ˆf [+ 1)][2][E][(][|][Y][ T][ +1][|][1+][c][)][.]

This completes the proof.

Define a stepping-stone function with respect to the distribution ρ[T][ +1]:


≲ _[M][c]_

_σ[c]_


_fη,T[∗]_ +1 [=] _j=1_ _fη,T[∗]_ +1,j [= arg min]f _∈HK_ [E][ℓ][σ][(][f] [(][x][T][ +1][)][ −] _[Y][ T][ +1][) +][ η]_

X


_τj∥fj∥K[2]_ _j_ _[.]_
_j=1_

X


To establish the bound of function approximation for stationary time series setting, we make the
following error decomposition.


-----

**Proposition 1. Let Assumptions 1-2 be true. For any f ∈** _H[¯]Z, there holds_

_f_ _f_
_∥_ _−_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ _[E][1][ +][ E][2][ +][ E][3][ + 2]σ[M][c][,][c]_

_where_
_E1 = {ET[σ] +1[(][f]_ [)][ −E]T[σ][(][f] [)][}][ +][ {E]T[σ][(][f][ ∗]η,T +1[)][ −E]T[σ] +1[(][f][ ∗]η,T +1[)][}]


_E2 = ET +1(fη,T[∗]_ +1[)][ −E][T][ +1][(][f][ ∗][) +][ η]


_τj∥fη,T[∗]_ +1,j[∥][2]Kj
_j=1_

X


_and_


_E3 = {ET[σ][( ˆ]f_ ) + λΩ( f[ˆ]) −ET[σ][( ˆ]fη) − _η_


_τj∥f[ˆ]η,j∥K[2]_ _j_ _[}][.]_
_j=1_

X


_Proof. According to Lemma 2, for any f_ _Z, we can make following error decomposition_
_∈_ _H[¯]_

_∥f −_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]) [=][ E][T][ +1][(][f] [)][ −E][T][ +1][(][f][ ∗][)][ ≤E]T[σ] +1[(][f] [)][ −E]T[σ] +1[(][f][ ∗][) +][ M]σ[c][c]

= _{ET[σ] +1[(][f]_ [)][ −E]T[σ][(][f] [)][}][ +][ {E]T[σ][(][f][ ∗]η,T +1[)][ −E]T[σ] +1[(][f][ ∗]η,T +1[)][}] +{ET[σ][(][f] [)][ −E]T[σ][(][f][ ∗]η,T +1[)][}]
_E1_
| {z } + _T +1[(][f][ ∗]η,T +1[)][ −E]T[σ] +1[(][f][ ∗][) +][ M][c]_

_E_ _[σ]_ _σ[c]_

_p_

_≤_ _E1 + ET +1(fη,T[∗]_ +1[)][ −E]T +1[(][f][ ∗][) +][ η] _τj∥fη,T[∗]_ +1,j[∥][2]Kj +{ET[σ][(][f] [) +][ λ][Ω(][f] [)][ −E]T[σ][(][f][ ∗]η,T +1[)]

_j=1_

X

_E2_

_p_

| {z }

_−_ _η_ _τj∥fη,T[∗]_ +1,j[∥][2]Kj _[}][ + 2]σ[M][c]_ _[c]_

_j=1_

X


+ [2][M][c]

_σ[c]_


_E1 + E2 + {ET[σ][(][f]_ [) +][ λ][Ω(][f] [)][ −E]T[σ][( ˆ]fη) − _η_


_τj∥f[ˆ]η,j∥K[2]_ _j_ _[}]_
_j=1_

X


_E3_
{z


_E1 + E2 + E3 + [2][M][c]_

_σ[c][ .]_


In statistical machine learning community, we call E1, E2 and E3 sample error, approximation
error and hypothesis error, respectively. The sample error E1 describes the divergence between
the empirical risk ET[σ][(][f] [)][ and the expected risk][ E] _[σ][(][f]_ [)][. The hypothesis error][ E][2][ characterizes the]
difference between the empirical regularized risks with HK and HZ. The approximation error
measures the approximation ability of RKHS HK to H.

B.1 SAMPLE ERROR

In this section, we focus on providing the bound of sample error E1. Unfortunately, the traditional
tools for complexity analysis such as covering number (Ron, 2000; Guo & Shi, 2011), Rademacher
complexity (Mohri & Rostamizadeh, 2009)) and concentration inequalities ((Wu et al., 2007)), cannot
be applied into this non-i.i.d. setting directly. To solve this problem, we employ the sequential
Rademacher complexity developed in (Rakhlin et al., 2010; Kuznetsov & Mohri, 2020a). We define a
function-based random variable as

_ξf_ (z) = ℓσ(f (x) − _y) −_ _ℓσ(fη,T[∗]_ +1[(][x][)][ −] _[y][)][, f][ ∈]_ _H[¯]Z._

Now, we turn to establish the bound of


Eξf (z[T][ +1])
_−_ _T[1]_


_ξf_ (z[t]), _f_ _Z._
_∀_ _∈_ _H[¯]_
_t=1_

X


-----

A necessary ingredient needed for our analysis is data-dependent sequential Rademacher complexity
Rakhlin et al. (2010), which we review in the following. We adopt the following definition of a
complete binary tree: a Z-valued complete binary tree v is a sequence (v1, ..., vT ) of T mappings
_vt :_ 1 _, t_ [1, T ]. A path in the tree is γ = (γ1, ..., γT 1) 1 . To simplify the
_{±_ _}[t][−][1]_ _→Z_ _∈_ _−_ _∈{±_ _}[T][ −][1]_
notation, we will write vt(γ) instead vt(γ1, ..., γt 1), even though vt depends only on the first t 1
_−_ _−_
elements of γ. The following definition generalizes the classical notion of Rademacher complexity
to sequential setting.

**Definition 3. The sequential Rademacher complexity RT (G) of a function class G is defined by**


_RT (G) = supv_ [E][[sup]g∈G


_γtg(vt(γ))],_
_t=1_

X


_where the supremum is taken over all complete binary trees of depth T with values in Z and γ is a_
_sequence of Rademacher random variables._

Based on the definition of sequential Rademacher complexity, we have the following concentration
inequality:

**Lemma 3. Suppose that the time series {Z** _[t]}t[∞]=−∞_ _[is strictly stationary. For any][ δ >][ 0][, with]_
_probability at least 1 −_ _δ, the following inequality holds for all f ∈_ _H[¯]Z and all α > 0_


8 log [1]


_Mℓσ_
4π log T _RT (G) +_


Eξf (z[T][ +1])
_≤_ _T[1]_


_ξf_ (z[t]) +
_t=1_

X


+ 6Mℓσ


_where G := {ℓσ(f, z) + ℓσ(fη, z) : f ∈_ _H[¯]Z} and Mℓσ is a positive constant such that ℓσ(f_ (x[(][t][)]) −
_y[(][t][)]) ≤_ _Mℓσ for any f ∈_ _H[¯]Z and t ∈_ Z.

According to the property of sequential Rademacher complexity (see Proposition 14 in (Rakhlin et al.,
2015)), we have


_RT (G) = RT (ℓσ ◦_ _H[¯]Z, z) ≤_ 16σ(1 + 4

with the 2σ-Lipschitz constant of ℓσ( ).

_·_

**Lemma 4. Under Assumptions 1-2, there holds**


2 log[3][/][2](eT [2]))RT ( H[¯]Z) (11)


1
2 log[3][/][2](eT [2])]M [2]p 2 κ


_RT (G) ≤_ [32][σ][[1 + 4]


_Proof. In fact, to every kernel Kj, we can associate a feature map φ with inner product <_
_φ(xj), φ(x[′]j[)][ >][=][ K][(][x][j][, x][′]j[)][ for any][ x][j][, x][′]j_

_[∈X][j][ and][ j][ = 1][, ..., p][. Then the constrained data-]_
dependent hypothesis space _H[¯]Z can be rewritten as_


_p_ 2M [2]

_τj_ _αj_ 2 _._
_j=1_ _∥_ _∥_ _≤_ _λT minj=1,..,p τj_ _}_

X


_αijφ(x[(]j[t][)][)][φ][(][x][j][) :]_
_t=1_

X


_H¯Z = {f_ (x) =


_j=1_


-----

By direct computation, we have

1
_T ( [¯]Z)_ = sup Eσ[ sup
_R_ _H_ _{x¯t,y¯t}t[T]=1_ _α∈R[T p]_ _T_

1
= sup Eσ sup

_T_ _{x¯[(][t][)]}t[T]=1_ _α∈R[T p]_

1
= sup Eσ sup

_T_ _{x¯[(][t][)]}t[T]=1_ _α∈R[T p]_

1
= sup Eσ sup

_T_ _{x¯[(][t][)]}t[T]=1_ _α∈R[T p]_


_αijφ(x[(]j[i][)][)][φ][(¯]x[(]j[t][)][(][σ][))][ −]_ _[y][(][t][)][)]]_
_i=1_

X


_σt(_
_t=1_

X


_j=1_


_T_

_αijφ(x[(]j[i][)][)][φ][(¯]x[(]j[t][)][(][σ][))][ −]_ [1]

_T_

_i=1_

X

_T_

_αijφ(x[(]j[i][)][)][φ][(¯]x[(]j[t][)][(][σ][))]_
_i=1_

X

_T_

_αijK(x[(]j[i][)][,][ ¯]x[(]j[t][)][(][σ][))]_
_i=1_

X


sup Eσ
_y¯[(][t][)]_ _t=1_
_{_ _}[T]_


_σty[(][t][)](σ)_
_t=1_

X


_σt_
_t=1_

X

_T_

_σt_
_t=1_

X

_T_

_σt_
_t=1_

X


_j=1_

_p_

_j=1_

X

_p_

_j=1_

X


_σtKt(σ), Kt(σ) = (K(x[(1)]1_ _[,][ ¯]x[(]1[t][)][(][σ][))][, ..., K][(][x][(]p[T][ )], ¯x[(]p[t][)][(][σ][))][′][ ∈]_ [R][T p]
_t=1_

X


sup Eσ sup
_{x¯[(][t][)]}t[T]=1_ _α∈R[T p][ α][′]_


_αsup∈R[T p][ ∥][α][′][∥]_ _{x¯supt}t[T]=1_ Eσ∥


_σtKt(σ)_ 2
_∥_
_t=1_

X


2M [2]

sup Eσ
_T_ [2] _x¯[(][t][)]_ _t=1_ _∥_
_{_ _}[T]_


_σtKt(σ)_ 2
_∥_
_t=1_

X


2M [2]

sup
_T_ [2] _x¯[(][t][)]_ _t=1_
_{_ _}[T]_

2M [2]

sup
_T_ [2] _x¯[(][t][)]_ _t=1_
_{_ _}[T]_

1
2M [2]p 2 κ


Eσ[ _σtσsK[′]t[(][σ][)][K][s][(][σ][)]]_
v
u _s,t=1_
u X
t


Eσ[ **K[′]t[(][σ][)][K][t][(][σ][)]]**
v
u _t=1_
u X
t


Thus, combining above result with Eq. (11). We get the desirable result.

Combining the results in Lemma 4 with the inequality in Lemma 3, we then obtain the bound of
sample error E1 for the strictly stationary time series setting.
**Proposition 2. Under Assumptions 1-2, for any δ > 0, the following inequality holds for all f ∈** _H[¯]Z_
_and all α > 0_


8 log [1]


_Mℓσ_
4π log TR +


_E1_ + 6Mℓσ
_≤_ _√T_

_with probability at least 1 −_ _δ, where_

_R = [32][σ][[1 + 4]_

B.2 APPROXIMATION ERROR


1
2 log[3][/][2](eT [2])]M [2]p 2 κ


Since the output function fη,T[∗] +1 [is the optimal estimator in RKHS][ H][K][ for time][ T][ + 1][, the learning]
rates of the learning algorithm indeed depend on the approximation ability of the hypothesis space _K_
_H_
with respect to the optimal risk (f ) measured by the approximation error E2. For any j = 1, ..., p
_E_ _[∗]_
and h Z, we first define the kernel integral operator LKj _,T +1 : L2(ρ[T][ +1]j_ [)][ →] _[L][2][(][ρ][T][ +1]j_ [)][ associated]
_∈_ _X_ _X_
with the kernel Kj by


_LKj_ _,T +1(f_ )(x[T]j [ +1]) =


_Kj(x[T]j_ [ +1], uj)f (uj)dρ[T][ +1]j [(][u][j][)][,][ ∀][h][ ∈] [Z][.]

ZXj _X_


-----

Note that LKj _,T +1 is a compact and positive operator on L2(ρ[T][ +1]j_ [)][. According to Mercer theorem,]
_X_

we can find the corresponding normalized eigenpairs {(ζh,i[j] _[, ψ]h,i[j]_ [)][}][i][≥][1][ such that][ {][ψ]h,i[j] _[}][i][≥][1][ is an]_
orthonormal basis of L2(ρ[T]X[ +1]) and ζh,i[j]
power L[r]Kj _,T +1_ [by] _[→]_ [0][ as][ i][ →∞][. Then for given][ r >][ 0][, we defined the][ r][-th]


_L[r]Kj_ _,T +1[(]_ _βh,i[j]_ _[ψ]h,i[j]_ [) =] _βh,i[j]_ [(][ζ]h,i[j] [)][r][ψ]h,i[j] _[.]_
Xi≥1 Xi≥1

We introduce an intermediate function as follows:
_f˜η,T +1,j = (LKj_ _,T +1 + ητjI)[−][1]LKj_ _,T +1fj[∗][,][ ∀][j][ = 1][, ..., p.]_

**Lemma 5. Under Assumption 3, for the intermediate function** _f[˜]η,T +1,j defined above, there holds_

E∥f[˜]η,T +1,j(x[T]j [ +1]) − _fj[∗][(][x]j[T][ +1])∥_ + ητj∥f[˜]η,T +1,j(xj[T][ +1])∥K[2] _j_ _[≤]_ [2(][ητ][j][)][2][r][∥][L]K[−][r]j _,T +1[f][ ∗]j_ _[∥]2[2][.]_


_Proof. Under Assumption 3, for any h ∈_ Z, we know that fj[∗] [=][ L]K[r] _j_ _,h[(][g]j,h[∗]_ [)][ for some][ g]h,j[∗] [=]

_i_ 1 _[β][h,i][ψ]h,i[j]_ _j_ [)][. Then we have]
_≥_ _[∈]_ _[L][2][(][ρ]X[T][ +1]_

P _fj[∗]_ [=][ L]K[r] _j_ _,T +1[(]_ _βh,iψh,i[j]_ [) =] (ζh,i[j] [)][r][β][h,i][ψ]h,i[j] _[.]_

Xi≥1 Xi≥1

The case r = [1]2 [means each][ f][ ∗]j [lies in the RKHS][ H][K]j [. Then we have]

_f˜η,h,j_ _fj[∗]_ [= (][L][K]j _[,h]_ [+][ ητ][j][I][)][−][1][L][K]j _[,h][f][ ∗]j_ _j_ [=] _ητj(ζh,i[j]_ [)][r][β]h,i[j] _[ψ]h,i[j]_ _._
_−_ _[−]_ _[f][ ∗]_ Xi≥1 _ζi[j]_ [+][ ητ][j]

Then we have


( _ητj(ζh,i[j]_ [)][r][β]h,i[j] )[2] = (ητj)[2][r][ X]( _ητj_ )[2][−][2][r]( _ζh,i[j]_ )[2][r](βh,i[j] [)][2]
Xi≥1 _ζh,i[j]_ [+][ ητ][j] _i≥1_ _ζh,i[j]_ [+][ ητ][j] _ζh,i[j]_ [+][ ητ][j]

(ητj)[2][r][ X](βh,i[j] [)][2][ = (][ητ][j][)][2][r][∥][L]K[−][r]j _,T +1[f][ ∗]j_ 2[.]

_i≥1_ _[∥][2]_


_∥f[˜]η,T +1,j −_ _fj[∗][∥]2[2]_ =

_≤_

Similarly, we also have


(ζh,i[j] [)][1+][r] _ζh,i[j]_ _ητj_

_ητj∥f[˜]η,T +1,j∥2[2]_ = (ητj)[2][ X]i≥1( _ζh,i[j]_ [+][ ητ][j] )[2](βh,i[j] [)][2][ = (][ητ][j][)][2][r][ X]i≥1( _ζh,i[j]_ [+][ ητ][j] )[2+2][r]( _ζh,i[j]_ [+][ ητ][j] )[2][−][2][r](βh,i[j] [)][2]

(ητj)[2][r] _L[−]K[r]j_ _,T +1[f][ ∗]j_ 2[.]
_≤_ _∥_ _[∥][2]_

**Proposition 3. Under Assumptions 3, there holds**


_E2 = E(fη,T[∗]_ +1[)][ −E][(][f][ ∗][) +][ η] _j=1_ _τj∥fη,T[∗]_ +1,j[∥][2]K _[≤]_ [(][p][ + 1)]

X


(ητj)[2][r] _L[−]K[r]j_ _,T +1[f][ ∗]j_ 2[.]
_j=1_ _∥_ _[∥][2]_

X


_Proof. According to Lemma 1, we have_

_ET +1(fη,T[∗]_ +1[)][ −E][T][ +1][(][f][ ∗][) +][ η]

_≤_ _ET[σ] +1[(][f][ ∗]η,T +1[)][ −E]T[σ] +1[(][f][ ∗][) +][ η]_

_≤_ _ET[σ] +1[( ˜]fη,T +1) −ET[σ] +1[(][f][ ∗][) +][ η]_

_T +1( f[˜]η,T +1)_ _T +1(f_ ) + η
_≤_ _E_ _−E_ _[∗]_


_τj∥fη,T[∗]_ +1,j[∥][2]K
_j=1_

X

_p_

_τj∥fη,T[∗]_ +1,j[∥][2]K [+][ M]σ[c][c]
_j=1_

X

_p_

_τj_ _fη,T +1,j_ _K_ [+][ M][c]
_∥_ [˜] _∥[2]_ _σ[c]_
_j=1_

X


_p_

_τj_ _fη,T +1,j_ _K_ [+ 2][M][c]
_∥_ [˜] _∥[2]_ _σ[c]_
_j=1_

X


-----

Then we have

_T +1( f[˜]η,T +1)_ _T +1(f_ ) + η
_E_ _−E_ _[∗]_


_j=1_ _τj∥f[˜]η,T +1,j∥K[2]_ = _∥f[˜]η,T +1 −_ _f_ _[∗]∥2[2]_ [+][ η]

X


_τj_ _fη,T +1,j_ _K_
_∥_ [˜] _∥[2]_
_j=1_

X


_j=1[ f[˜]η,h,j −_ _fj[∗][]][∥]2[2]_ [+][ η]

X


_τj_ _fη,h,j_ _K_
_∥_ [˜] _∥[2]_
_j=1_

X


_j=1_ _p∥f[˜]η,T +1,j −_ _fj[∗][∥]2[2]_ [+][ η]

X


_τj_ _fη,T +1,j_ _K_
_∥_ [˜] _∥[2]_
_j=1_

X


(ητj)[2][r] _L[−]K[r]j_ _,T +1[f][ ∗]j_ 2[.]
_j=1_ _∥_ _[∥][2]_

X


(p + 1)


Combing above results, we get the desired result.

B.3 HYPOTHESIS ERROR

This section focuses on bounding the hypothesis error

_E3 = {ET[σ][( ˆ]f_ ) + λΩ( f[ˆ]) −ET[σ][( ˆ]fη) − _η_


_∥f[ˆ]η,j∥K[2]_ _[}][.]_
_j=1_

X


We first give the properties of _f[ˆ] and then use them to bridge_ _f[ˆ] and fη._

**Lemma 6. For all j = 1, ..., p, there holds**


_τj∥αj[η][∥][2][ = 1]ηT_


**_Kjt[′]_** _[α]j[η][))][2][ ≤]_ [2][σ]1
_j=1_ _ηT_ 2

X


(ℓ[′]σ[(][y][t][ −]
_t=1_

X


_Proof. Recall the represent theorem which ensures that_


_j=1_ _t=1_ _αtj[η]_ _[K][j][(][x]j[t]_ _[,][ ·][)][, α]tj[η]_ _[∈]_ [R][.]

X X


_fˆη =_


For notation simplicity, denote αj[η] [= (][α]1[η]j[, ..., α]T j[η] [)][′][ ∈] [R][T][ and][ α][η][ = ((][α]1[η][)][′][, ...,][ (][α]p[η][)][′][)][′][ ∈] [R][T p][.]
From the definition of Eq. (7), we deduce that


_T_ _p_

**_α[η]_** = arg min _ℓσ(y[t]_ **_Kjt[′]_** _[α][j][) +][ η]_
_α∈R[T p][ {][ 1]T_ _t=1_ _−_ _j=1_

X X

where Kj = _K(x[s]j[, x][i]j[)][}][n]i,s=1_
_{_ _[∈]_ [R][T][ ×][T][ . we then have]


_τj(αj[η][)][′][K][j][α][j][}][,]_
_j=1_

X


_ℓ[′]σ[(][y][t][ −]_
_t=1_

X


**_Kjt[′]_** _[α]j[η][)][K]ji[′]_ [=][ ητ][j][K][j][α]j[η][.]
_j=1_

X


It is easy to deduce that

_αj[η]_ [=]


1

_σ[(][y][1][ −]_
_ητjT_ [(][ℓ][′]


**_Kj[′]_** 1[α]j[η][)][, ..., ℓ]σ[′] [(][y][T][ −]
_j=1_

X


**_KjT[′]_** _[α]j[η][))][′][.]_
_j=1_

X


Then we obtain that, for any j = 1, ..., p,


_τj∥αj[η][∥][2][ = 1]ηT_

This completes the proof.


_p_

**_Kjt[′]_** _[α]j[η][))][2][ ≤]_ [2][σ]1
_j=1_ _ηT_ 2

X


(ℓ[′]σ[(][y][t][ −]
_t=1_

X


-----

**Proposition 4. Under Assumptions 3, there holds**


_p_

_τj_ _αj[η]_ 1
_j=1_ _∥_ _[∥][2][ ≤]_ [2]ηT[λpσ]2 _[.]_

X


_E3_ _λΩ( f[ˆ]η) = λ_
_≤_


_Proof._


_E3_ = _{ET[σ][( ˆ]f_ ) + λΩ( f[ˆ]) −ET[σ][( ˆ]fη) − _η_


_∥f[ˆ]η,j∥K[2]_ _[}]_
_j=1_

X


= _ET[σ][( ˆ]f_ ) + λΩ( f[ˆ]) −ET[σ][( ˆ]fη) − _λΩ( f[ˆ]η) −_ _η_

_λΩ( f[ˆ]η)_
_≤_

Combining the above inequality with Lemma 6, we get that


_∥f[ˆ]η,j∥K[2]_ [+][ λ][Ω( ˆ]fη)
_j=1_

X


_τj_ _αj[η]_ 1
_j=1_ _∥_ _[∥][2][ ≤]_ [2]ηT[λpσ]2 _[.]_

X


_E3_ _λΩ( f[ˆ]η) = λ_
_≤_


**Proof of Theorem 1: Combining Propositions 1-4, for Tη ≤** 1, we have with confidence 1 − _δ_

_f_ _f_ _C log(1/δ)(T_ 2 + _T_ _σ +_ _η[2][r]_ + _λη[−][1]T_ 2 σ + _σ[−][c]λ[−][2]T_ + _σ[−][c]λ[−][1]T_ ),
_∥_ _−_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ [e] _[−]_ [1] _[−][1]_ _[−]_ [1] _[−][2]_ _[−][1]_

where _C is a positive constant independently of T, λ, η, δ, σ. By taking σ = T_ _[m], η = T_ _[β]_ and
_λ = T_ _[γ], we then have with confidence 1 −_ _δ_

[e] _f_ _f_ _C log(1/δ)T_ [Ψ(][m,β,γ,c,r][)],

_∥_ _−_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ [e]

where

Ψ(m, β, γ, c, r) = max
_{−_ 2[1] _[, m][ −]_ [1][,][ 2][rβ, γ][ −] [1]2

By taking σΨ( =m, c, r21c [,][ β][ =]) =[ −]4[1]rmaxmax[and]{−{−[ γ][ =][1]221 _[,][, m][ −][ 1]2c[ −]4[1][−]r_ _[−][1][1][,][}][ −]2[,][1]c[cm][Direct computation shows that][ +][−][ 1]2r[β][+ 2][ +][ m,][m][ −][ −][cm][2][}][,][ −]_ ifif[2][γ] m m >[ −] ≤[2][,] 1[ −]1 − −[cm]44[1][1][ −]rr _[.]_ _[γ][ −]_ [1][}][.]

This completes the proof.

C PROOF OF THEOREM 2


Denote by s = {st}t[T]=1 [a probability set with][ P]t[T]=1 _[s][t][ = 1][. We define following functions associated]_
with the probability set s:

_p_

_fT[∗] +1_ [=] _j=1_ _fT[∗] +1,j_ [= arg min]f _∈HK_ [E][ℓ][σ][(][f] [(][x][T][ +1][)][ −] _[y][T][ +1][)][,]_

X


_fT[∗][,][s]_


_j=1_ _fT,j[∗][,][s]_ [= arg min]f _∈HK_

X


_stℓσ(f_ (x[t]) _y[t]),_
_−_
_t=1_

X


_stℓσ(y[t]_
_−_
_t=1_

X


_fˆη[s]_ [=] _fˆη,j[s]_ [=] arg min

_j=1_ _f_ =[P][p]j=1 _[f][j]_ _[,f][j]_ _[∈H][K]j_

X


_fj(x[t]j[)) +][ η]_
_j=1_

X


_τj∥fj∥K[2]_ _j_ _[}][,]_
_j=1_

X


and


_fˆ[s]_ = _j=1_ _fˆj[s]_ [=] _f_ =[P][p]jarg min=1 _[f][j]_ _[,f][j]_ _[∈H][K]j_ _{t=1_ _stℓσ(y[t]_ _−_ _j=1_ _fj(x[t]j[)) +][ λ][Ω(][f]_ [)][}][.]

X X X

Correspondingly, we denote by ET[σ,][s][(][f] [) =][ P]t[T]=1 _[s][t][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][)][.]_


-----

**Proposition 5. Let Assumptions 1-2 be true. For any f ∈** _H[¯]Z, there holds_

_p_

_∥f −_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1]) _[≤]_ _[E][1][ +][ E][2][ +][ E][3][ +][ E]T[σ][(][f][ ∗]T_ [) +][ η] _j=1_ _τj∥fT,j[∗]_ _[∥][2]Kj_ [+ 2]σ[M][c][,][c]

X

_where_
_E1 = ET[σ] +1[( ˆ]f_ **[s]) −ET[σ,][s][( ˆ]f** **[s])**


_E2 = ET[σ,][s][( ˆ]f_ **[s]) + λΩ( f[ˆ][s]) −ET[σ,][s][( ˆ]fη[s][)][ −]** _[η]_


_τj_ _fη,j[s]_ _[∥][2]Kj_
_∥_ [ˆ]
_j=1_

X


_and_


_E3 = ET +1(fη,T[∗]_ +1[)][ −E][T][ +1][(][f][ ∗][) +][ η]


_τj∥fη,T[∗]_ +1,j[∥][2]Kj _[.]_
_j=1_

X


_Proof. According to Lemma 2, for any f_ _Z, we can make following error decomposition_
_∈_ _H[¯]_

_∥f −_ _f_ _[∗]∥L[2]_ 2(ρ[T]X[ +1])

= _T +1(f_ ) _T +1(f_ )
_E_ _−E_ _[∗]_

_T +1[(][f]_ [)][ −E]T[σ] +1[(][f][ ∗][) +][ M][c]
_≤_ _E_ _[σ]_ _σ[c]_

_≤_ _{ET[σ] +1[(][f]_ [)][ −E]T[σ,][s][(][f] [)][}] +{ET[σ,][s][(][f] [)][ −E]T[σ,][s][(][f][ ∗]T[,][s][)][}][ +][ E]T[σ,][s][(][f][ ∗]T[,][s][) +][ E]T[σ] +1[(][f][ ∗]η,T +1[)][ −E]T[σ] +1[(][f][ ∗][)]
_E1_

_p_

| {z }

+ η _τj∥fη,T[∗]_ +1,j[∥][2]Kj _[}][ +][ M]σ[c][c]_

_j=1_

X


_E1 + {ET[σ,][s][(][f]_ [) +][ λ][Ω(][f] [)][ −E]T[σ,][s][(][f][ ∗]T[,][s][)][ −] _[η]_


_τj∥fT,j[∗][,][s][∥]K[2]_ _j_ _[}][ +][ E]T[σ,][s][(][f][ ∗]T[,][s][) +][ η]_ _τj∥fT,j[∗][,][s][∥]K[2]_ _j_
_j=1_ _j=1_

X X


_p_

+ ET[σ] +1[(][f][ ∗]η,T +1[)][ −E]T[σ] +1[(][f][ ∗][) +][ η] _τj∥fη,T[∗]_ +1,j[∥][2]Kj _[}][ +][ M]σ[c][c]_

_j=1_

X


_E1 + {ET[σ,][s][(][f]_ [) +][ λ][Ω(][f] [)][ −E]T[σ,][s][( ˆ]fη[s][)][ −] _[η]_


_τj∥f[ˆ]η,j[s]_ _[∥][2]Kj_ _[}]_
_j=1_

X


+ET[σ,][s][(][f][ ∗]T[,][s][) +][ η] _τj∥fT,j[∗][,][s][∥]K[2]_ _j_

_j=1_

X


_E2_
{z


+ ET[σ] +1[(][f][ ∗]η,T +1[)][ −E]T[σ] +1[(][f][ ∗][) +][ η] _τj∥fη,T[∗]_ +1,j[∥][2]Kj _[}][ +][ M]σ[c][c]_

_j=1_

X


_p_

_τj∥fT,j[∗][,][s][∥]K[2]_ _j_ [+ 2]σ[M][c] _[c]_
_j=1_

X


_E1 + E2 + ET +1(fη,T[∗]_ +1[)][ −E]T +1[(][f][ ∗][) +][ η]


_τj∥fη,T[∗]_ +1,j[∥][2]Kj
_j=1_

X


+ET[σ,][s][(][f][ ∗]T[,][s][) +][ η]


_E3_

_p_

{z

_τj∥fT,j[∗][,][s][∥]K[2]_ _j_ [+ 2]σ[M][c] _[c]_
_j=1_

X


_E1 + E2 + E3 + ET[σ,][s][(][f][ ∗]T[,][s][) +][ η]_


_τj∥fT,j[∗][,][s][∥]K[2]_ _j_ [+ 2]σ[M][c][ .][c]
_j=1_

X


_E1 + E2 + E3 + minf_ _Z[{E]T[σ,][s][(][f]_ [) +][ λ][Ω(][f] [)][}][ +][ η]
_∈H_


Following the generalization analysis in (Kuznetsov & Mohri, 2020a), we need to discrepancy
measure to measure the discrepancy of the non-stationarity of the stochastic process {Zt}t[∞]=−∞ [with]
respect to both the loss function ℓσ and the hypothesis set _H[¯]Z._


-----

**Definition 4. The discrepancy describes the discrepancy between target distribution and the distribu-**
_tion of the sample. For any f_ _Z, the discrepancy measure with respect to Huber loss is defined_
_∈_ _H[¯]_
_as_


Eℓσ(f (x[T][ +1]) − _y[T][ +1]) −_


Estℓσ(f (x[t]) − _y[t])_
_t=1_

X


disc(s) := sup
_f_ _∈HZ_


**Lemma 7. For any δ > 0, with probability at least 1 −** _δ, the following inequality holds for all_
_f ∈_ _H[¯]Z and all α > 0_


8 log [1]

_δ [,]_


Eℓσ(f (x[T][ +1])−y[T][ +1]) ≤


_stℓσ(f_ (x[t]) _y[t])+disc(s)+_ **s** 2+6Mℓσ
_−_ _∥_ _∥_
_t=1_

X


4π log T _RT (G)+Mℓσ_ _∥s∥2_


_where G := {ℓσ(f, z) : f ∈_ _H[¯]Z} and Mℓσ is a positive constant such that ℓσ(f_ (x[t]) − _y[t]) ≤_ _Mℓσ_
_for any f ∈_ _H[¯]Z and t ∈_ Z.

**Proof of Theorem 2 The proofs of bounding errors E1, E2 and E3 proceeds similarly to the proof of**
Theorem 1 and are omitted for brevity.

According to Lemma 3, we have


8 log [1]

_δ [,]_


_E1 ≤_ disc(s) + ∥s∥2 + 6Mℓσ

where
_R =_ **s** 2[32σ[1 + 4
_∥_ _∥_

Moreover, from Proposition 4, we can get


4π log T + Mℓσ **s** 2
_R_ _∥_ _∥_

1
2 log[3][/][2](eT [2])]M [2]p 2 κ].


(ητj)[2][r] _L[−]K[r]j_ _,T +1[f][ ∗]j_ 2[.]
_j=1_ _∥_ _[∥][2]_

X


_E3_ (p + 1)
_≤_

Similarly, to bound the hypothesis error


_E2 = {ET[σ,][s][( ˆ]f_ **[s]) + λΩ( f[ˆ][s]) −ET[σ,][s][( ˆ]fη[s][)][ −]** _[η]_


_∥f[ˆ]η,j[s]_ _[∥]K[2]_ _[}][,]_
_j=1_

X


we first give the properties of _f[ˆ] and then use them to bridge_ _f[ˆ] and fη._

**Lemma 8. For all j = 1, ..., p, there holds**

_E2_ _._
_≤_ [2][λpσ]η[∥][s][∥][2]

_Proof. Recall the represent theorem which ensures that_


_j=1_ _t=1_ _αtj[η,][s][K][j][(][x]j[t]_ _[,][ ·][)][, α]tj[η,][s]_ _[∈]_ [R][.]

X X


_fˆη[s]_ [=]


For notation simplicity, denote αj[η,][s] = (α1[η,]j[s][, ..., α]T j[η,][s][)][′][ ∈] [R][T][ and][ α][η,][s][ = ((][α]1[η,][s][)][′][, ...,][ (][α]p[η,][s][)][′][)][′][ ∈]
R[T p]. We then deduce that


_stℓσ(y[t]_
_−_
_t=1_

X


**_α[η,][s]_** = arg min
_α_ R[T p][ {]
_∈_


_τj(αj[η][)][′][K][j][α][j][}][,]_
_j=1_

X


**_Kjt[′]_** _[α][j][) +][ η]_
_j=1_

X


where Kj = _K(x[s]j[, x][i]j[)][}][n]i,s=1_
_{_ _[∈]_ [R][T][ ×][T][ . We get]


_stℓ[′]σ[(][y][t][ −]_ **_Kjt[′]_** _[α]j[η,][s][)][K]ji[′]_ [=][ ητ][j][K][j][α]j[η,][s][.]
_t=1_ _j=1_

X X


-----

It is easy to deduce that

_αj[η,][s]_ =


1

(s1ℓ[′]σ[(][y][1][ −]
_ητj_


**_Kj[′]_** 1[α]j[η][)][, ..., s][T][ ℓ]σ[′] [(][y][T][ −]
_j=1_

X


**_KjT[′]_** _[α]j[η][))][′][.]_
_j=1_

X


Then we obtain that, for any j = 1, ..., p,


_τj_ _αj[η]_
_∥_ _[∥][2][ = 1]η_ vu

u
t

This completes the proof.

Under Assumptions 1-2, there holds


**_Kjt[′]_** _[α]j[η][))][2][ ≤]_ [2][σ][∥][s][∥][2]

_η_

_j=1_

X


_s[2]t_ [(][ℓ][′]σ[(][y][t][ −]
_t=1_

X


_E2_ = _{ET[σ][(][f]_ [) +][ λ][Ω( ˆ]f **[s]) −ET[σ][( ˆ]fη[s][)][ −]** _[η]_


_∥f[ˆ]η,j[s]_ _[∥][2]K[}]_
_j=1_

X


_ET[σ][( ˆ]f_ **[s]) + λΩ( f[ˆ][s]) −ET[σ][( ˆ]fη[s][)][ −]** _[λ][Ω( ˆ]fη[s][)][ −]_ _[η]_


_∥f[ˆ]η,j[s]_ _[∥]K[2]_ [+][ λ][Ω( ˆ]fη[s][)]
_j=1_

X


_λΩ( f[ˆ]η[s][) =][ λ]_

2λpσ **s** 2
_∥_ _∥_ _._


_τj_ _αj[η,][s]_
_j=1_ _∥_ _[∥][2]_

X


By combining the above results, we have with confidence 1 − _δ_

_f_ **[s]** _f_
_∥_ [ˆ] _−_ _[∗]∥L[2]_ 2(ρ[T]X[ +1])
_≤_ disc(s) + minf _Z[{E]T[σ,][s][(][f]_ [) +][ λ][Ω(][f] [)][}]
_∈H_

+C log(1/δ)(∥s∥2 + ∥s∥2σ + η[2][r] + λη[−][1]∥s∥2σ + σ[−][c]λ[−][2]T _[−][2]_ + σ[−][c]λ[−][1]T _[−][1]_ + η),

1
where _C is a positive constant independently of T, λ, η, δ, σ and p. By taking σ = T_ 2c, λ = T _[−][1]_

and η = T _[−][e][1]2, we completes the proof._

[e]

D PROOF OF THEOREM 3

According to the definition of disc(T + 1), we have


Eℓσ(f (x[T][ +1]) − _y[T][ +1]) −_


_stEℓσ(f_ (x[t]) − _y[t])_
_t=1_

X


disc(s) = sup
_f_ _∈HZ_

_≤_ _fsup∈HZ_


_qt[∗][E][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][)][ −]_
_t=1_

X


_stEℓσ(f_ (x[t]) − _y[t])] + disc(q[∗])._
_t=1_

X


Furthermore, it is easy to deduce that


(qt[∗]
_t=1_ _[−]_ _[s][t][)][E][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][)][ −]_ _f[sup]∈HZ_

X


(qt[∗]
_t=1_ _[−]_ _[s][t][)][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][)]_

X


sup
_f_ _∈HZ_

sup
_f_ _∈HZ_


(qt[∗]
_t=1_ _[−]_ _[s][t][)[][E][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][)][ −]_ _[ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][)]][.]_

X


**Proof of Theorem 3 According to the proof of Theorem 1 in Kuznetsov & Mohri (2020a), we have**


8 log [1]


(qt[∗]
_t=1_ _[−][s][t][)[][E][ℓ][σ][(][f]_ [(][x][t][)][−][y][t][)][−][ℓ][σ][(][f] [(][x][t][)][−][y][t][)]][ ≤]

X


4π log T _RT (G)+Mℓσ_ _∥q[∗]−s∥2_


sup
_f_ _∈HZ_


+6Mℓσ


-----

By combining above results, we have with confidence 1 − _δ_

_f_ **[s]** _f_
_∥_ [ˆ] _−_ _[∗]∥L[2]_ 2(ρ[T]X[ +1])

_T_

_≤_ disc(q[∗]) + supf _∈HZ_ _t=1(qt[∗]_ _[−]_ _[s][t][)][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][) + min]f_ _∈HZ[{E]T[σ,][s][(][f]_ [) +][ λ][Ω(][f] [)][}]

X

+C log(1/δ1 )(∥q[∗] _−_ **s∥2 + ∥s∥2 + ∥s∥2σ + η[2][r]** + λη[−][1]∥s∥2σ + σ[−][c]λ[−][2]T _[−][2]_ + σ[−][c]λ[−][1]T _[−][1]_ + η).

By taking σ = T 2c, λ = T _[−][1]_ and η = T _[−]_ 2[1], we complete the proof.

[e]


**Algorithm 1: Optimization procedure for adaptive SpHAM**

**Input: Data {(x[t], y[t])}t[T]=1[, Max-Iter][ Z][ ∈]** [Z][, Mercer kernel][ K][j][, j][ = 1][, ..., p][ with bandwidth][ d][,]
Weights τl, l = 1, ..., p, q[∗].
**Initialization: Lipschitz constant L, s[0].**
**Step A: Computing weights ˆs:**
**for z = 1, ..., Z do**

1. Compute A[z][−][1] via DC-programming (or gradient descent method);
2. Update s[z] via (12).


**Output: ˆs = s[Z].**
**Step B: Computing** _f[ˆ][s]:_
**for z = 1, ..., Z do**

1): Compute α[z] = pL(β[z]) via (14);

2): mz+1 = 1+[√]1+42 _m[2]z_ ;

3): β[z][+1] = α[z] + _[m]m[z]z[−]+1[1]_ [(][α][z][ −] _[α][z][−][1][)][.]_

**Output: α[ˆ]s = αZ;**
**Prediction function:** _f[ˆ][s]_ = _j=1_ _Tt=1_ _[α][ˆ]tjs_ _[K][j][(][x][t]j[,][ ·][)][;]_
**Variable selection:** _j :_ _αj[s][∥][2][ ≥]_ _[v, j][ = 1][, ..., p][}][.]_
_{_ _∥_ P

[P][P]

E EXPERIMENT OPTIMIZATION


The optimization problem (10) reduces to problem 7 when taking s[∗] = _T[1]_ [I][T][ . Recall the optimization]

problem in Step A:


_T_

(qt[∗] 2 [+][ λ][2] 2[T]
_t=1_ _[−]_ _[s][t][)][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][) +][ λ]2[1]_ 2

_[∥][q][∗]_ _[−]_ **[s][∥][2]** _[∥][s][∥][2]_

X


1

2c }.


ˆs = arg min
**s** _[{]f[ sup]∈HZ_


The above optimization problem can be equivalently rewritten as a common type of bilevel optimization problem (Colson et al., 2007), i.e,

Outer problem: mins _Tt=1[(][q]t[∗]_ _[−]_ _[s][t][)][ℓ][σ][(][f][s][(][x][t][)][ −]_ _[y][t][) +][ λ][1][∥][q][∗]_ _[−]_ **[s][∥]2[2]** [+][ λ][2][T][ 1][/][2][c][∥][s][∥]2[2][,]

_T_
Inner problem: fs = arg maxP _f_ **z** _t=1[(][q]t[∗]_
_∈H_

_[−]_ _[s][t][)][ℓ][σ][(][f]_ [(][x][t][)][ −] _[y][t][)][,]_

where the outer (min) problem parameterized byP **s, is nested within the inner (max) problem. The**
outer problem can be solved by standard gradient descent, where in each step, we need to optimize
the inner (max) problem with last updated s.

We denote by k the iteration time. Then for k + 1-th iteration, we have the following gradient update
rule

1
**s[k][+1]** = s[k] _γ(A[k]_ _λ1(q[∗]_ **s[k]) + λ2T** 2c s[k]), (12)
_−_ _−_ _−_

where γ is learning rate,


_p_

_′_
_αtj[s][k]_ _[K][j][(][x]j[t]_ _[, x][T]j_ [)][ −] _[y][T][ )]_ R[T]
_∈_
_j=1_

X 


_A[k]_ = _ℓσ(_
_−_



_αtj[s][k]_ _[K][j][(][x]j[t]_ _[, x][1]j_ [)][ −] _[y][1][)][, ..., ℓ][σ][(]_
_j=1_

X


_t=1_


_t=1_


-----

and
_α[s][k]_ = (α11[s][k] _[, ..., α]T[s][k] 1[, ..., α]1[s][k]p[, ..., α]T p[s][k]_ [)][′][ ∈] [R][T p]

is obtained by the following weighted optimization problem


_α[s][k]_ = arg max


_t=1(qt[∗]_ _[−]_ _[s]t[k][)][ℓ][σ][(]m=1_

X X


_αmjKj(x[m]j_ _[, x]j[t]_ [)][ −] _[y][t][)][.]_
_j=1_

X


Note that the inner problem is subjected to f ∈HK, i.e., ∥α∥2 = _j=1_ _[τ][j][∥][α][j][∥][2][ ≤]_ _λT min2jM=1[2],...,p τj_ [.]

A widely-used method for solving this constrained inner problem is DC programming (Tao & An,
1998). For simplicity, we here transform this inner problem into a regularized problem with[P][p] _ℓ2-norm_
regularizer, and solve this regularized problem by standard gradient method. After obtaining the
solution ˆs, we turn to solve the following weighted optimization problem in Step B:

_T_ _p_ _p_

_α[ˆ]s =_ _αj_ _∈arg minR[T]_ _,j=1,...,p{t=1_ _sˆtℓσ(y[t]_ _−_ _j=1(K[t]j[)][′][α][j][) +][ λ]_ _j=1_ _τj||αj||2}._ (13)

X X X

We can see that this problem contains non-smooth function λ _j=1_ _[τ][j][||][α][j][||][2][ which makes the]_
standard gradient descent method inapplicable. To conquer this challenge, we leverage fast iterative
shrinkage-thresholding algorithm (FISTA)(Beck & Teboulle, 2009). Our optimization problem

[P][p]
becomes


_T_

_sˆt_ _ℓσ(y[t]_ (K _[t])[′]α[k][−][1])_
_∇_ _−_ _||[2][o]_
_t=1_

X


_τj_ _αj_ 2 + _[L]_
_||_ _||_ 2 _[||][α][ −]_ [(][α][k][−][1][ −] _L[1]_


_α[k]_ = pL(α[k][−][1]) := arg min


_,_ (14)
1≤j≤p

_T_

_sˆt_ _αj_ _ℓσ(y[t]_ (K _[t])[′]α[k][−][1]))_
_∇_ _−_
_t=1_

X


1 −


_vj_


_vj_ 2
_||_ _||_


where

and


_vj = αj[k][−][1]_ _−_ _L[1]_

_αj_ _ℓσ(y[t]_ (K _[t])[′]α[k][−][1])) =_
_∇_ _−_


2((K _[t])[′]α −_ _y[t])(Kj[t][)][′]_ _|(K_ _[t])[′]α −_ _y[t]| < σ_

_αj_ _ℓσ(y[t]_ (K _[t])[′]α[k][−][1])) =_ 2σ (K _[t])[′]α_ _y[t]_ _σ_
_∇_ _−_  _−_ _≥_

 2σ (K _[t])[′]α_ _y[t]_ _σ_

_−_ _−_ _≤−_

and L = max _T[1]_ _Tt=1_  _T_ _Tt=1_

Proj(s) as a projection of[∥][(][K] s and[t][)][′][K] v >[t][∥] [is the lipschitz constant] 0 as the threshold value. To focus on the weights of samples[ 1] _[∇][ℓ][σ][(][y][t][ −]_ [(][K] _[t][)][′][α][k][−][1][)][.][ Denote]_
which are really useful for forecasting, we consider ProjP (s)P as a box projection such that st >
_qt[∗][,][ ∀][t][ = 1][, ..., p][). Finally, the optimization procedure for adaptive SpHAM can be summarized in]_
Algorithm 1. Note that if we only run Step B with ˆst = _T[1]_ _[,][ ∀][t][ = 1][, ..., T]_ [, we can further obtain the]

optimization procedure for SpHAM.
**Remark 3. The computational complexity of Algorithm 1 depends on the optimization strategy for**
_DC programming, the training size T_ _, dimension p and the iteration times Z. We denote by the_
_computational complexity of DC programming O(DC(T, p)). Then, the computational complexity of_
_Step A is O(ZDC(T, p) + ZT_ ) and the computational complexity of Step B is O(ZT [3]p[3]). Thus, the
_total computational complexity of Algorithm 1 is O(ZO(DC) + ZT_ [3]p[3]). For large scale data, we
_can further speed up Algorithm 1 by random Fourier features technique [Rahimi and Recht 2007],_
_which is leaved for future work. This has been carefully discussed in the revised manuscript._


F ADDITIONAL EXPERIMENT

F.1 EVALUATION ON BENCHMARK DATA

We test our algorithm on nonlinear dataset from CauseMe. The hyper-parameter selection is the same
as the one in Section 4. The results in Table 3 verify the effectiveness of our method.


-----

Table 3: The ASE on CauseMe data (p refers to the dimension of features).

Methods (p = 3, T = 300) (p = 5, T = 300)

LSTM 0.7681 **0.9230**
TS−SpAM 0.7782 0.9485
SpHAM **0.7548** 0.9484

F.2 EXPERIMENTS ON AIR QUALITY DATASET

[We use the Air Quality dataset obtained from UCI Machine Learning Repository (https://](https://archive.ics.uci.edu/ml/datasets/Air+quality)
[archive.ics.uci.edu/ml/datasets/Air+quality) to test our model’s ability to detect](https://archive.ics.uci.edu/ml/datasets/Air+quality)
the Granger causality. This dataset includes 9358 hourly air quality data in an Italian city, collected
from March 2004 to February 2005. The details of the dataset can be obtained on the UCI website.
The Granger causal network we detect is shown in Figure 3. Given the network, we observe that
the temperature (T) influences ozone (O3) and nitrogen dioxide (NO2) which is validated in Kalisa
et al. (2018). Mwaniki et al. (2014) confirm the relationship between relative humidity (RH) and
nitrogen dioxide (NO2). Yan et al. (2018) verify the relationship between humidity and nitrogen
dioxide (NO2).

F.3 EXPERIMENTS ON CORONAL MASS EJECTIONS DATASET

Coronal Mass Ejections (CMEs) are the most violent eruptions in the Solar System. Despite
machine learning approaches have been applied to these tasks recently Wang et al. (2019); Liu
et al. (2018), there is no any work for interpretable prediction with Granger causal network. CMEs
[data are provided in The Richardson and Cane List (http://www.srl.caltech.edu/ACE/](http://www.srl.caltech.edu/ACE/ASC/DATA/level3/ icmetable2.htm)
[ASC/DATA/level3/icmetable2.htm). From this link, we collect 152 ICMEs observations](http://www.srl.caltech.edu/ACE/ASC/DATA/level3/ icmetable2.htm)
[from 1996 to 2016. The features of CMEs are provided in SOHO LASCO CME Catalog (https:](https://cdaw.gsfc. nasa.gov/CME_list/)
[//cdaw.gsfc.nasa.gov/CME_list/). In-situ solar wind parameters can be downloaded](https://cdaw.gsfc. nasa.gov/CME_list/)
[from OMNIWeb Plus (https://omniweb.gsfc.nasa.gov/). A total of 9 features are chosen](https://omniweb.gsfc.nasa.gov/)
as input, including: (1) Central PA (CPA), (2) Angular Width, (3) three approximated speeds ( Linear
Speed, 2nd-order Speed at final height,and 2nd-order Speed at 20 Rs), (4) Mass, (5) Kinetic Energy,
(6) MPA and (7)CMEs arrival time. Figure 4 shows that Granger causal network when the output
is CMEs Arrival time. Some interesting findings are concluded from this Granger causal network.
For instance, Speed and Mass, as the significant variables causing CMEs arrival time, have been also
screened out in Liu et al. (2018). Morover, the CMEs Angular Width does not cause the CMEs arrival
time forecasting, while Liu et al. (2018) state that they have a significant correlation. This indicates
that the CMEs arrival time is affected by the Angular Width at current time, but not by the historical
Angular Width.


-----

Figure 3: The Granger causal network on Air Quality dataset.


-----

Figure 4: The Granger causal network on CMEs dataset.


-----

