Huber Additive Models for Non-stationary Time Series Analysis
Authors: Yingjie Wang
Xianrui Zhong
Fengxiang He
Hong Chen
Dacheng Tao
Venue: ICLR 2022 Poster (2022)
Keywords: Sparse additive models
variable selection
Huber
non-stationary
robust forecasting

Abstract: Sparse additive models have shown promising ﬂexibility 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 ﬁnance and medical ﬁelds. 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.

Published as a conference paper at ICLR 2022
HUBER ADDITIVE MODELS FOR NON-STATIONARY TIME SERIES ANALYSIS
Yingjie Wang Xianrui Zhong", Fengxiang He Hong Chen & Dacheng Tao? 'College of Informatics, Huazhong Agricultural University, China 2ID) Explore Academy, JD.com Inc, China
Department of Computer Science, University of Illinoist Urbana- Champaaan,,,,,n USA *College of Science, Huazhong Agricultural University, China yingrieeangg1201 gmail 11. com, xzhong23@1l1ing fengxi ang. f. he@gmail 1. com, chenh@ma i 11 hzau. edu. cn, da dacheng tao i com
edu,,
ABSTRACT
Sparse additive models have shown promising flexibility and interpretability in pro- cessing 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
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, propose an adaptive sparse Huber additive model for robust forecasting in both non- non-Gaussian data and (non)sttaiinss generalization bounds of our estimator are established for both stationary and non- stationary time series data, which are independent of the widely used mixing condi- tions in learning theory of dependent observations. Moreover, the error bound for non non- stat-ttttonayyy 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 penalty in our method. Experimental results on both synthetic and real-world benchmark datasets validate the effectiveness of the proposed method. com/ wwwpiiiiaiiiiiiaa www.rruuizoogggioo SPHMM
data. In theory, the
The code is available at
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 WWng & Yang [0077] Chu & Glymour 2008 Song Yang 2010
Yang et al.) 2018 Chen et al. [2017] Liu al. |020 Chen et al. 222 la) and practical applications
(Dominici et al. 2002] Wang & Brown 2011 I, Ravindra et al.) 2019 Bussmann et al. 2020; Wang
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 al. 2015)) develop an elliptical vector autoregressive model for estimating heavy-tailed stationary processes with parametric convergence analysis that reduces the influence of non- GaussGausian innovation. Moreover, under stationarity and decaying (-mmxxig condition, WWong et al. (2020) derive nonasympootic 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.

---

Published as a conference paper at ICLR 2022
sequences. Agarwal & Duchi
(2033) 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 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, first propose novel sparse Huber additive model (SpHAM) with Huber loss, sparsity-inducing 42, -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) simultancously. In theory. Section 3.2 establishes the upper bound of the function approximation error of SpHAM by developing error decomposition technique (Wu et
2006 Chen et al. 2021b) and employing sequential Rademacher complexity (Rakhlin
et al. 2010 2013 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 b) can be pursued, even if the innovation is non- Gaussian distribution (see Theorem | and Corollary 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 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.
Shrinkage-Threshoddhe
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 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.
Related works:
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 I summarizes the algorithmic properties of our method and other related works, e.g., Sparse additive models for time series (TS-SpAM) (Yang et al.
2

---

Published as a conference paper at ICLR 2022
Table 1: The properties of related methods
TS-SpAM
DBF
R-Dantzig
Ours
Additive
Additive
Hypothesis space
Kerne--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-base based Forecasting (DBF) Kuznetsov & MMohri 2020a) and Robust Dantzig-
selector-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 Section 4. Finally, Section concludes this paper. The source code package available at httts:://sat http::://tttttt com www.iiiaaiiiiiiao stteuuiiihnng//pppea
2 PRELIMINARY
Let {- be a stochastic time series with time index t. where variable Z = takes values in the compact input space X' c RRP and the output space y We consider common nonparametric model
=x= B(e)) =0
where (*()) is the ground truth function, and the innovation e is i.i.d. across time / E Z. For the sake of simplicity, we denote p² and p'x as the jointed distribution of (X*, Y*) and the corresponding marginal distribution with respect to X respectively. This setup actually covers a large number of scenarios commonly used in practice. For instance, the case X that contains p lagged values of y (e.g., X* = (t-p) x ... x Y ((-1)) corresponds to the p-order autoregressive models. Moreover, this case can be viewed as vector autoregressive model the sense that input X includes the historical information of multiple variables.
Although such nonparametric model (I) 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 space X' E R into X' = X1 x X Xp. Under the assumption that the ground truth admits an additive structure f = > fj. the additive model can be defined by
-(=)
+5'
where each component f3 X₂ + 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) 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 ( is strictly stationary if (Z) ..... Z²) and (2l+k Zttk) have the same joint distributions for any E Z and Z.
Note that, if not otherwise stated, the stationarity in this paper refers to strict stationarity. Suppose that we are given size time series data (()",")) z which are drawn from an additive data-generating model (2). Under stationarity condition and zero-mean Gaussian innovation with finite variance, the widely used methods that learn the ground truth f usually integrate squared loss and smoothness- -aassity-nduciig regularizer 2(-) into structural risk minimization scheme:
P
min FEH
Σ
(g'
(())) j=1
3

---

Published as a conference paper at ICLR 2022
where 7/ := {fr + ....
+ +
Hy.J 1...,,pp is an additive hypothesis space. Commonly, each
subspace 7 could be reproducing kernel Hilbert space (Chen et al. 2017 Kandasamy & Yu 2016
Raskutti et al. (2012). the orthogonal basis inducing space (Ravikumar et al. 2009; Meter et al. 2009
Yang et al. 2018), or the composite function space with the neural network as a typical (Agarwal
al. 2020 Bussmann et al. 2020).
However, when facing heavy-tailed innovation, these methods may have degraded performance due the amplification of the squared loss 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
(((())) - v²) = [ff(") - y²)² f(x') - yO1 -?? if If(x²) y" 2 a,
if If(x²) 3/1 <o
where a is positive hyper- papre-.armmeter. Note that in the previous studies (Huber & Ronchetti (009)] Loh, 2017), the hyper-parameter a is set to be fixed according to the 95% asymptotic efficiency rule. However, Huber regression with 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, choose the scale parameter by relating 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.
METHOD
3.1 SPARSE HUBER ADDITIVE MODELS
In this paper. we choose reproducing kernel Hilbert space (RKHS) HK,.j = L.... to form the the additive hypothesis space H. where each Hx, associated with symmetric and positive semi-definite Mercer kernel Kj Xj X3 - R. An additive RKHS is defined as
HK = (/- + + Sp E HK,jj= = 1,...
with kernel norm || = inf (1/kk,, + |1111 J. 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
f, where 17 is positive regularization parameter and Tj is the weight for j-th kernel norm. The representer theorem (Wahba 1990) ensures that f, can be represented as
arg min fin
lo - ED
Σ j=1
(5)
f E J.RK t=1
T
f
a o
K,
og ER.
j=
To offer the method with sparsity, consider the following sparsity-inducing penalty
S(
inf
Σ
ary K, K
j=1
j=1
Let Hz be an additive data dependent hypothesis space defined by
T ar anjKy(z), auj E R).
Hz = (f= =
Σ
1=
Then the SpHAM can be formulated
1 Σ
P lo -
- arg min(
+ (Q(1))
T
JEHz
4

---

Published as a conference paper at ICLR 2022
Denote Kj (Ky KGlf,,, R a/))...,,''' RTP and aj (anj..... OrTy'' E R², where of here refers to the transpose of aj for avoiding the confusion. The SpHAM can be represented as
ai (z))))
(6)
j=
with
T
a
arg min 0y f
Ca
EiKYas) + 3 1=1
T [1
(7)
The optimization problem 77 can be solved by Fast Iterative Bhrinkage-Thresholee
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 al. 2012 Pan & Zhu 2017 Chen et al. [017)) 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 {KzYi=ı E R with the subgroups decomposition {Xa) d= where each Na X may cover more than one variable.
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) !
Section 3.2.1) and non-stationary time series data (Theorem 2 Section 3e.22on3.2.22):
An adaptive SpHAM (inspired by above theoretical findings) and its function approximation error for non-stationary time series (See Theorem Section 3.2.3)
Due to space limitation, the high-level outline and detailed proofs are provided the Appendix A
3.2.1 FUNCTION APPROXIMATION ANALYSIS FOR STATIONARY TIME SERIES
Through this paper, the marginal distribution w.r.t Y is assumed to be almost everywhere supported on M, .) for some M 2 O.
Assumption Assume that |Y Vt E Z, is bounded and there exists constant c >> such that Z.
B| < oo, Vt
The moment condition in Assumption is rather weak in the sense that the response variable y possesses infinite variance. The same condition also applies to the distributions of the innovation st. implying that heavy-tailed innovation is allowed.
Assumption 2. Let supzex / (x, x) < oc, Vj = 1, ......
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 al. 2006 Steinwart & Christmann 2008 Wu & Zhou 2008).
The following definitions are needed for Assumption For any j = I,....P. p, we define kernel T+1 > T=2(2%)t) T+11 associated with the kernel K2 by integral operator I.K,TT+I T22
/ K2 uz) (1, dp (uz).
LK,T (1)cf)))
x,
T-1 ). According to Mercer theorem, we Note that Lx, is compact and positive operator on L2(p can find the corresponding normalized eigenpairs {(C) v )22 such that { J22 is an orthonormal basis of (2(T++) T+1) and ç + 0 as + cc. Then for given r > 0, we defined the r-th power L Lk, by
T
1>1
5

---

Published as a conference paper at ICLR 2022
Assumption 3. We assume that fj X3 R. 1....P is function of the form fj Lkk, LKk,T+109)) Vr (0, bll with some oz E L2(ex) wh &.
Assumption 3 is natural extension from i.idd setting (Assumption in Wu et al. 2006 Christmann & "Zhou, 2016, Chen Wang 2018) to non-i.i. time series. Indeed, this assumption stands in most practical cases. For instance, if 0.5, the ground truth function f* needs to be real-valuee function the RKHS Hx. The RKHS admits a large class of bounded and continuous real-valued 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 )= is stationary. Let Assumptions 1-3 be true. By taking o T ,7 T d. and X = T m,, we have for any 0 < 5 < 1,
Il 1- L s (Clo((1/))TMMmnn
with confidence at least - 8, where Ö is a positive constant independently of X, â and a and
max 5. m 1, -cm + 1 + m - - if m S 1 1
((m,crr)=
max m -m + 2m 2}. if > 1 4r
In the stationary case, bound appears to be novel for the following reasons: a) the result 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 QQu et 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.
Remark 2.
Corollary 1. Suppose that the process [2= is stationary. Let all the conditions in Theorem be true. We then have for any 0 < <<1
Ilf - 12llatty 1-
C (log(L18))
with confidence least 8, where
[maxf -cm 1 m max m -m + +2m
(m,
1}.
if m
2
Figure T summaries the convergence rates in Corollary by taking differ- ent a and c. If the innovation Gaus- sian distribution with finite variance holds for any 1), one arbitrarily select
2
00)
(i.e., Assumption
0(r))
o = T (0 < m < ( to obtain con- vergence rates O(T- $). Moreover, we can see that the convergence rate will decrease as increases. Com- bined with the conclusion Lemma 1 (i.e., the equivalence relation be- tween Huber loss based empirical risk and MSE as a oo), it in- dicates that indeed plays trade- off role between algorithmic robust- ness and variance-reduucion. For the weak moment condition (e-g., <
O cmn
Ofgli-ii
001 1)
value for r
Figure 1: The convergence rates under different and
1), may get slower conver- gence rates (e.g., O(T(I-c)m- 4) or
(2 )), which also coincides with our intuitive understanding that small a may be conducive
to robust forecasting, Note that our method will not converge when o and c are both located in the white area in Figure]

---

Published as a conference paper at ICLR 2022
As 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 al. (2018); Han al. (2013); Wang & Yang (007): Kock & Callot (2015)). However, from Theeoem and Corollary 1 the asymptotic convergence SpHAM can be obtained under weaker moment condition, which verifies the robustness of our method. Recall =2d
learning rate O(T derived Yang et al. (2018), where the order of smoothness of
the component function fj.j = 1, ,....P p. Relatively slow convergence rate OTT²)) we obtained indicates the sacrifice for the absence of mixing condition.
3.2.2 FUNCTION APPROXIMATION ANALYSIS FOR NON- NON-STATIONARY TIME SERIES
In nonstationary time series setting, different Z's may follow different distributions. Thus For given weights {st [ we define the estimator f E}=1 ff the minimizer of the following weighted objective
P l (v'
{I(ft) = min
+(()))+000)
JEHz
Next, we introduce discrepancy measure to characterize the discrepancy between the target distribu- tion and the distributions of observations (Kuznetsov & Mohri 2020a).
Definition 2. For any f E Hz, the discrepancy measure with respect to Huber loss is defined as
disc(s): = sup fEHz
Bl, Blo TT+ Be -f((t
The discrepancy measures the non- non-stationarity of the stochastic process { ft=-oo with respect to both the loss function l, and the hypothesis Hz.
Theorem 2. Let Assumptions 1-3 be true. We assume that fj E HK, Vj I,.... By taking
== T A
T-1 and 17 = T -, we have for any 0 <
|lf fr 1' ( I disc(s) + =E(())) + E 2| + log(1/887
with confidence at least 1 8, where Ü₁ C2 are two positive constants independently Tf X, n, and a.
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 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 illustrates that we shall minimize the following optimization problem for non-stationary time series forecasting:
T E
min JEHz
sclo
Σn(a)) disc(s) + x1 + X |||2)) 3=1
where A and X2 are two positive regularization parameters, and e is positive constant introduced in Assumption [ 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', E 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 pT+1 and distributions p¹, = 1,...??
7

---

Published as a conference paper at ICLR 2022
Assumption 4. Denote by a probability set q" {qi]e=) with El= o8 1. We assume that the following term is sufficiently small:
disc(q") = sup [Bb JEHz
y Σ Bl -/)]l
Note that the priori q' 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 does not change
drastically compared with the distributions ,², 1..,, T. Kuznetsov & Mohri (2020a) have proven that the Assumption 4 holds if the s* is an uniform distribution over last > observations, where / a apper-parameter that can be tuned in practical applications.
Theorem 3. Let Assumption F and the conditions in Theorem
2 be true. We huve for any 0 < <
11 - MI (
Σ (g (o)s))(((((()))) - v²) + Ex((²)
sup JEH2
Er log with confidence at least - 8, where G₁ C2 are two constants independently of X, n, 8 and a.
log((1/))) s + |
Finally, the optimization problem of adaptive SpHAM can be formulated as following two stages:
Step A: finding the weight S:
$ = arg m sup JEHz t=1
Σ(E - s,)co(((())) Ailla" - -s + To |
Step B: forecasting:
p T
(TT11)
a
ax +111,
j=1t=1
where
a
EyER" ;I Ealll' Σ(K))a)) - A Σ 7311051123-
arg min
(10)
The optimization problems can be solved by standard gradient descent method. Similarly to the strategy for optimization problem (7). we use Fast Iterative Shrinkage-Thresholle Thresholding Algorithm (FISTA) (Beck & Tebbulle 2009) for Step B. We provide the detailed procedure in Appendix
EXPERIMENT
This section validates the effeetiveness of SpHAM and adaptive SpHAM. In all experiments, the -
Gaussian kernel (,, exp(- 2d for constructing the additive data dependent hypothesis space. Due limited space, the evaluations on real-world data are provided in Appendix
where j 1, I,... and bandwidth d > 0, is employed
We consider two synthetic examples as below:
Example A: Inspired by (Kuznetsov & Mohri 2020a). stationary time series is generated according to the non-linear additive autoregressive model:
yt sin " y1-2) sin( yt 2
where the innovation e's are i.i.d. drawn from Gaussian distribution NO, 0.5) and Student distribution with degree of freedom 2, respectively.
Example B: Inspired by (Kuznetsov & Mohri 2020a). time series with smooth drift is generated by
yt_ 400
sin y 1) + 2 =
8

---

Published as a conference paper at ICLR 2022
where the distributions of noise e are the same above.
Shyper-pramette selection and evaluation criterions: Recall that SpHAM algorithm requires three pyper-parameters: regularization parameter X, bandwidth of kernel d and Huber parameter We set these parameters according to the suggestions our Theorems. Based on the suggestion in Thcorem 3 the selection of Huber parameter a is o = T and the regularization parameter is À = 7-1, Moreover, we set the bandwidth d 0.5 and tune / {100, 150, 200, ..., 350). The evaluation
criterions for forecasting used here contains Average Sample Error(ASE)= 7े E - y²)²
and True Deviation (TD)= N V E (f(x') ft (x*))? where N is the number of test samples.
For each example, we generate time series with 4000 sample points. The samples at time t {1500, 1501, ]899} are used as training set, and the samples next time f {1900,.... .... 1999) considered as the test data. The competitors include Simple Exponential Smoothing (SES). TS- TS-SpAM (Yang et al. 2018) and Vanilla Long Short-Tene 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
From T 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 & 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(2.008)
0.098(E.005)
0.493(E.313)
0.105(E.009)
LSTM
0.06((+.005)
0.022(=.005)
0.538(+.312)
0.219(+.089)
Example A
TS-SpAM
0.070(+.0007)
0.035(#.008)
0.505(+.332)
0.131(+.042)
SpHAM (ours)
0,076(+.005)
0.04((#.008)
0.492(+.333)
0.095(+.006)
SES
0.435(+.0115)
0.42((+.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.17((:.125)
0.104(:..05)
Example
SpHAM (ours)
0.115(+.028)
0.112(#.029)
0.179(+.116)
0.096(+.026)
Adantive SpHAM (ours)
0.102(++.02)
0.096(+0026)
0.172(+.118)
0.085(+.024)
5 CONCLUSION
We propose an adaptive sparse Huber additive model by integrating Huber loss and l 1 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 Ar- tificial Intelligence" key project (No. 20217D01 11700) 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): 2010.
Alekh Agarwal and John C. Duchi. The generalization ability of online algorithms for dependent data. IEEE Transactions Information Theory, (91))55335587 2013.

---

Published as a conference paper at ICLR 2022
Rishabh Agarwal, Nicholas Frosst, Xuezhou Zhang, Rich Caruana, and Geoffrey E. Hinton. Neural additive models: Interpretable machine learning with neural nets. arXiv:2004.1391211. 2020.
Richard T. Baillie, Long memory processes and fractional integration in econometrics. Journal of Econometrics, 73(1);5-59, 1996.
Amir Beck and Mare Teboulle. A fast iterative thrikkige-thresholding algorithm for linear inverse problems. SIAM Journal on Imaging Sciences, 2(1):183-200) 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- 5228, 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) (:3311))232222333 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-25) 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. wwXiv:2009.12255v1. 2020.
Zheng-Chu Guo and Lei Shi. Classification with non-i.idd sampling. Mathematical and Computer 54(5):1347-1366 5
2011.
Fang Han, Huanran Lu, and Han Liu. direct estimation of high dimensional stationary vector autoregressions. Journal of Machine Learning Research, 16:3115-3150 2015.
Trevor J. Hastic 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:20/2.10931, 2020.
10

---

Published as a conference paper at ICLR 2022
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-47) 477, 2004.
Peter J. Huber. Robust estimation of location parameter. The Annals of Mathematical Statistics, 35 (1):73-101, 1964.
Peter J. Huber and Elvezio M. Ronchetti. Robust Statisfics. 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 high dimensional nonpara- metric regression via the SALSA. In International Conference Machine Learning (ICML), pp- 69-78, 2016.
Anders Bredahl Kock and Laurent Callot. Oracle inequalities for high dimensional vector autoregres- sions. Journal of Econometrics, 186(2):325-38) -444,, 2015.
Vitaly Kuznetsov and Mehryar Mohri. Generalization bounds for time series prediction with non- stationary processes. In International Conference on Algoritumic Learning Theory, 2014.
Vitaly Kuznetsov and Mehryar Mohri. Discrepancy- caareppance-baeedd -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 Learming (ICML), 2020.
Jiajia Liu, Yudong Ye, Chenglong Shen, Yuming Wang, and Erdelyi. A new tool for eme 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 mestima- tors. 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 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 d-mixing and 8-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 regres- sion. 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 Neural Networks and Learning Systems, 3096): 2019.
Huitong Qiu, Sheng Xu, Fang Han, Han Liu, and Brian Caffo. Robust estimation of transition matrices in high dimensional heavy-tailes vector autoregressive processes. In Proceedings of the 32nd International Conference Machine Learning, volume 37, pp. 1843-1851, 2015.
11

---

Published as a conference paper at ICLR 2022
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 (111::111153,,, 2015.
Garvesh Raskutti, Martin J. Wainwright, and Yu. Minimax-optimal rates for sparse additive models 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-1033, 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, Sehastian Bathiany, Erik M. Bollt, Gustau Camps- Valls, Dim Coumou, Ethan R Deyle, Clark Glymour, Marlene Kretschmer, Miguel D. Mahecha, Jordi Muñoz-Marí, Eghert H. van Nes, Jonas Peters, Rick Quax, Markus Reichstein, Marten Scheffer, Bernhard Schölkopf, 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 autore- gression models with simultaneous confidence band. Journal of Multivariate Analysis, 101(9): 2008--025, 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, 1 115(529):254-205 2019.
Pham Dinh Tao and Le Thi Hoai An. 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. Jajja Liu, Ye Jiang, and Robert Erdélyi. Cme arrival time prediction using convolutional neural network. The Astrophysical Journal, 881 (11):15, 2019.
12

---

Published as a conference paper at ICLR 2022
Yingjie Wang. Hong Chen, Feng Zheng. Chen Xu, Tieliang Gong, and Yanhong Chen. Multi-task additive models for robust estimation and automatic structure discovery. 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. JEEE Transactions on Cybernetics, 47(6):1354-1366, 2017b.
Kam Chung Wong, Zifan Li, and Ambuj Tewari. Lasso guarantees for B-mixing heavy-tailed time series. The Annals of Statistics, 48(2):1124-11444 2020.
Qiang Wu and Ding Xuan Zhou. Learning with sample dependent hypothesis spaces. Computers and Mathematics with Applications, 56(11):8896-2977, 2008.
Qiang Wu, Yiming Ying, and Ding-Xuan Zhou. Learning rates of least-square regularized regression. Foundations of Computational Mathematics, 6(2):171-19, 2006.
Qiang Wu, Yiming Yin, and Ding-Xuan Zhou. Multi-kerne 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 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. 3xXi::80330391111 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())2275-295, 2009.
13

---

Published as a conference paper at ICLR 2022
A
PROOF SKETCH
This section sketches the proof focusing on the conceptual aspects. The full proofs are provided in Appendix B D
Theorem (Approximation Error Bound)
Proposition 5
Proposition 2
Proposition (Hypothesis Error)
Proposition 4
(Error Decomposition)
(Sample Error)
(Approximation Error)
Lemma
Lemma
Lemma 3
Lemma
Lemmass
Lemma
Figure 2: The important ingredients for Theorem
Figure summaries the important ingredients for the proof of Theorem Lemma describes the properties of and Lemma illustrates the relation between huber-based risk and MSE. Following these two lemmas and previous works (Wang et al. 2017a Feng & Wul 2020), we then decompose the function error into three important parts: sample eтror, hypothesis error and approximation error. To bound the sample error, we need to employ the sequential Rademacher complexity (Kuznetsov Mohri (020a), which is the key to measuring the capacity of the data-dependent hypothesis space for non-i.i.d data. Lemmas 3 give an important concentration inequality and the upper bound of sequential Rademacher complexity, respectively. Finally, Theorem is obtained by combining Propositions 2-4. The proof of Theorems
2 similar to the above process, except that we develop a different error decomposition in Proposition On the basis of Theorem by applying assumption on data distribution, we derive Theorem
3
The full proofs of Theorems
provided in Sections C-E, respectively.
B PROOF OF THEOREM
We first illustrate two key properties of
Lemma 1. Let Assumptions 1-2 be true. From the definition of f in Eq. 66. there hold
2x²M
11 S M
X7
and
V Σ
2M
2())
Fllا 1l2 S Sx 2M
XI minn........p
j=l
where M is a positive constant such that |4 M.
Proof. Denote by 87(F) Ei ()) for notational convenience. From the definition
of f. we know that
Eq(1) + + $
means that
8
2 - 5
14

---

Published as a conference paper at ICLR 2022
From the definition of Huber loss
obtain
T
(7)
so? if 1 St Vt I,...
and
T Σ
()) la (y)
20/ 0² < 20 - -o?, if 1 VI = =--,.... I.....
1=1
Therefore,
p
2M
Σ Trlla311 S
XI m'iijj_l....p Tj
According to the property of RKHS, we conclude that
2=2 Σ Elll Sk² Σ laglla X7 minj mi... Tj
l S skllflk S K
This completes the proof by denoting M
2x
min:-1......77
Let Hz be data-dependent hypothesis space constrained by sparsity regularizer 22()). i.e.
T 1z= = - ony (5 mz)) Σ Tj layll2
2M²
XT Xmminj=1...pp TJ }.
J=1
j=1
For any h E 2, denote by (l respectively. In the following, describe the relationship between (7f and (f), vf Hz. Note that a similar proof is given in (Feng Wu 2020) for the i.i.d case.
and (T) = B( B[ T+212
Lemma 2. Let Assumptions 1-2 be true. For any Z. and f Hzz there holds
M
(0 ET
[17
ge
where both and M. 23t My 1)) are positive constants.
Proof.
For any a > max(M M + Mf. 1}, we denote two events Iy and IIy as follows
lytn (Y |T+ 2 0/2)
and
1ly++1 = (T |YT < o/22.. Then, for any f E Fl and ||0 = M1- < oc, we have
|Cf+1(0) - -7+1((") - (E - Er+((()
ITʃ - 1d - 7
|I 1'1 (
(t I1 - (o T+1 11 dpx T+1 - 1| - 1 M
T+1
Vel Spprrsultt (1))=))) - ) (o(5 (x x+1) - T+1 dd T+1 dp
+ [ I 1yT- (f(xT))) - Tt1p - ( T1T T+ dp T+1
$ 1 Jay c(((x741) +4 - "-(1" (1tt)) y T+1 d2 +1
+| 1 1YT-1
(f(ot+))- T2222 ( T
15

---

Published as a conference paper at ICLR 2022
The last inequality is based on the fact that
T-1 f(x < 1 T+ + llfllo <0, v(I , [+]) E & 1lyr-1
and hence
Lo (f(x , (f(x*t)) y ,
Recalling that the Huber loss is 20- 2i-Lipschitz continuous, then we get
[+ U IYTHI (a( - - - ( -d
$20| 20 1x ly7+1 1f(x) - (x)) T+1
$ -f I
According to Markov's inequality, for any constant c > O we further have
2
ProbfIyr+a+ Prob[[y?Tt] 2 a/2) Prob Prob(l (0/2)
altc
Then we have
U (((((r+t)) - T+ - (f - y+ d T+1 2f -f
o= Moreover, the second term in the Light-hand side can be bounded by
1/ Jy (f(xIt1) - +2+ (x y d T+1
dp
$ - -11--51 1 [ 1y741
12y T+1 - [(x - ( (xx [d T+1
S llf - 1 /
2/ 21y T+1 + (11fl +
lyT+1
According to Holder inequality, we have
2² PFB(YT+114+
J "T+""dpp"
THr (B)))))
Then deduce that
|l - 1- 1 lyr+
2/1/ 1dp T+1 + (11/1>> + [1/"1l1>)
(22+c + 2
2 115 17
gc
oltc
Finally, there holds
11/ la (1
(6( - - M
Me $
ge
with
M. 23 (Mp- + Mj + 1)
This completes the proof.
Define stepping- stone function with respect to the distribution TT!
M Σ Eh arg min Bl, (f. -Y jERx
- + Ta all
To establish the bound of function approximation for stationary time series setting, we make the following error decomposition.
16

---

Published as a conference paper at ICLR 2022
Proposition 1. Let Assumptions 1-2 be true. For any f Hz, there holds
2
|l f 1' L21p3 - Ei + Eg + Ez +
ac
where
E (8741 - - + (EF - E7+ (1)++11)
Es E -0r+((f+)) + Σ Tylf |lk,
and
Es (E7(f) ( - n
Σ Ty flkg
Proof. According to Lemma 2, for any / E Hz, we can make following error decomposition
Ilf - -lloaee ) (2+1lf) -++1(f) - (000 Me oe
-1 -
E1
+8++((144)) (7+1 + Ma go
< Ei + - - + Σ T li, +187(1) AN(F) (f(
8g
f1 1i,,t 2Mc
" Σ
P
2M. re
< E B ( - "TIlff lix, }
3=1
2M. 0=
< E1 + F2 + Fag +
In statistical machine learning community, call Er, E2 and E3 sample error, approximation error and hypothesis error, respectively. The sample error E₁ describes the divergence between the empirical risk (f) and the expected risk (f). The hypothesis error Eg 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 providing the bound of sample error E,. Unfortumately, the traditional tools for complexity analysis such as covering number (Ron) 2000 Guo Shi 201 1). Rademacher complexity (Mohri & Rostamizadeh, 2209)) and concentration inequalities ((Wu et al. (007)). 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 function- random variable as
Cy(=) (f(x)) v) lo (1) v). Hz.
Now, we turn to establish the bound of
T Σ
EE, (2 (xT+1)
Er(:)) Vf
T
17

---

Published as a conference paper at ICLR 2022
A necessary ingredient needed for analysis is data-dependent sequential Rademacher complexity Rakhlin et al. (2010), which we review in the following. We adopt the following definition
complete binary tree: Z-valued complete binary tree is sequence (U1,,.., ..., UT) of T mappings (+ + 2, E [1, T A path in the tree is 7 (71..., -r-1) E [=1]7- To simplify the notation, will write uy (*)) instead (",.... "t-1). even though depends only on the first
elements of y. The following definition generalizes the classical notion of Rademacher complexity sequential setting.
Definition 3. The sequential Rademacher complexity RT(G) of a function class G is defined by
Rr(g)= = sup R Suu 1 geg T
(g)))
where the supremum is taken over all complete binary trees of depth T with values in Z and Y is sequence of Rademacher random variables.
Based on the definition sequential Rademacher complexity, we have the following concentration inequality:
Lemma 3. Suppose that the time series = is strictly stationary. For any 8 > 0, with probability at least - 8, the following inequality holds for all f E Hz and all > 0
T +
+ +M,, v 4T log (Rr(R) + Mr, v log 1
1 VT
<
Ef(²")
T
where 9 = 2l) + Lo 2) Hz2 and Me, is positive constant such that to (f y(t)) < Me, for any Hz and Z.
(fm
According to the property of sequential Rademacher complexity (see Proposition 14 in Rakhlin et al. 2015
we have
(r) Rr(C, H2, I 16o
(11)
with the 2o- constant of l (-).
Lemma 4. Under Assumptions 1-2, there holds
RT( < 32o + 4 220g (cT2)( 2
T
Proof. fact, to every kernel Ky, we can associate a feature map with inner product < ((x;), o(*)) >= x(zj,z)) for any zj 2/ E X2 and j 1,...,.. Then the constrained data- dependent hypothesis space Hz can be rewritten as
P -==(f(a) = ΣΣ (() j=
2M2 Σ Tallalla XT min; Tj
18

---

Published as a conference paper at ICLR 2022
By direct computation, have
T
T
(xz( = sup
E sup
0
(oll))")(2 y/1)]
T
(21.8:)T
P
[=]
J
1 sup
Eg sup (=W)
P
ayol2 )0(0)) (o))
(
sup E. Σ (o)
T
DER
1-11
j=
1 sup
Eg sup "R P
a )
T
(BWN
T
1
K( (0 (o))
=
sup
Eg sup HR
aiy
T
{B00 It-1
1=1
1=1
1
T (o), K2(o) = (x)
= 1
sup
Eg
sup OERtY
( E ERT
T
[=1
1
sup TERT
sup |a''l
Eg
T
(E. 1
1=1
2M
<
sup (2(1TT
EG
(o)l(
T
1t=1
t=1
2M
sup (MONTT
Eg
oi 0,o, [(o,K(()),((t))
T
2M
T Eg
sup (3(1TT
(((o)
T
1=1
2M 222
T
Thus, combining above result with Eq. ITT. 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 E, for the strictly stationary time series setting.
Proposition 2. Under Assumptions 1-2, for any O, the following inequality holds for all f E Hz and all a >
Mr
8log
1 + +M,, V 4 log TR +
E S
VT
VT
with probability at least 1 6, where
R 32a + 4V2log 3/2 (
T
B.2 APPROXIMATION ERROR
Since the output function fi,4 is the optimal estimator in RKHS HR for time T + 1, the learning rates of the learning algorithm indeed depend on the approximation ability of the hypothesis space H. with respect to the optimal risk E((*) measured by the approximation error Eg. For any j I..., and h E Z. we first define the kernel integral operator LK, Le(p - Le(px) T+1 associated with the kernel K, by
LKgT+I (k,0+t()(( Jz8 (u)) (1;) wh EZ
19

---

Published as a conference paper at ICLR 2022
(1 ). According to Mercer theorem, we can find the corresponding normalized eigempairs {(S v 22, such that (v%,1221 is an orthonormal basis of L2(xx*) and Sh. + 0 as i + co. Then for given r > O, we defined the r-th power Lx T++ by
Note that L1,T+11 is compact and positive operator on L2(Px;
L K5,T+1 +1, C (h. Chanl
121
121
We introduce an intermediate function as follows:
fu.++11 = (Lx,,T41 + nT, LKDD"L,,,,+il Vj I.....P.
Lemma 5. Under Assumption 3. for the intermediate function fo.T+1, defined above, there holds
(=] (1
+ (1, (xT+1) lll -
)a1 "''")
Proof. Under Assumption 3, for any h E Z, we know that fz = T'K, n (03,,) for some Dh
Eii1 She E L2(ex;;)) ). Then we have
1j LK ( E
$21
The case = means each fÿ lies in the RKHS HK,- Then we have
7T (r
Inns - f3 = (LKgn + (1 LKGNfj - f3 =
C +1Ts
121
Then we have
nT (T,((h "h (1TT;)21
Sha,t
-
TTj
)2-2r(
)ar (8hh))))
Sl.i + +1T
Sh + 17Tj
Ch. t15
121
121
(nTTy)
(3) 2 (17T7)
T+1 TT,,,
Similarly, we also have
(Ch 1+r
(1T1)) )
)²(3),)) (1T7))2
)222r(
TTg
j2- 2r (8)
Sh.s !5द
Sh.5 T/TJ
Ch. + 1TT
12
(1TT))2T1
K,1 112.
Proposition 3. Under Assumptions 3 there holds
Eg = -
- E(f)) + n Tillf 11 S (o+ ++11 1) >
(nTT,
|| 2
5=1
Proof. According to Lemmal I we have
p
(r+1(TTTt1) - - + Σ
Tyll T
1=1
M ge
< (f+1((iT+41) - c0 (T+1( +
T 11
1=1
Ma
< - +1)
- -77
ge
j=1
2M.
s Er41(fo,411 -
+
20

---

Published as a conference paper at ICLR 2022
Then we have
(f+1(far+11) - Er+1(f") +
Σ
Tflff 1l = - 5 + Σ
==1
==1
p 1Eini 171122 +
Ty 11 11
j=1
Σ E -3 + Σ j=1
Ty T
==1
(p + 1) Σ (1Tr)) 2r | +1 f2
j=1
Combing above results, we get the desired result.
B.3 HYPOTHESIS ERROR
This section focuses on bounding the hypothesis error
Es (ET(f) + xR(f) (fn) Σ lik).
j=1
We first give the properties of f and then use them to bridge f and Sn.
Lemma 6. For all j = 1, I. P, there holds
1
Σ Cl Ekzaa))33 S
20
n7T
j=1
Proof. Recall the represent theorem which ensures that
Jn
aly K,(,),)),, al ER.
For notation simplicity, denote a² = (a a E RT and a² ((a?)),.... e R?2
From the definition of Eq. 7 we deduce that
p
2 arg min
(v'
EK
+1 Σ Tr(aa)KKajj) KjaJh, j=1
oERTY
where K, = (K(2)
E RR
TXT
we then have
T
Σ
- Ka,
1=1
is easy to deduce that
o3 NT,T
Σ
K PoV Σ
K Kor
1=1
Then we obtain that, for any j I,....p
T
3110ح112=
1
Σ (y'5"' - EkSay))
20 nf
nT
1=1
5=1
This completes the proof.
21

---

Published as a conference paper at ICLR 2022
Proposition 4. Under Assumptions
3 there holds
2Apo
Es S (l(fa)= = > 2
7/T
Proof.
P (87()) + (2() - -8((T)) -n -nΣ
Es
1mm
ei(++ + x2( 8q((m) 1M(f)) Σ
1inilll + M1(()))
AQ(fx) Combining the above inequality with Lemma 6 get that
2xpo T/
Es S An(s))
Ty llla 122 -
Proof of Theorem 1: Combining Propositions 1-4, for T' S we have with confidence 1 1-6 || - M S (Clos(/))((++++>>>
+ Am -T - 2 ------
where Ö is a positive constant independently of T, ,, By taking Tm, 11 = T and T7, we then have with confidence 1 5
|| - f 12 Lo
where
(1,,B,7,,,r) max{ 5. 1,,,,,21 1,2r,,
3 m, -c - 27 - 2, -m - - 11..
By taking =-=-1 -- and ?
Direct computation shows that
2c
maxx
21 -1},
1 1
if < 1
(m,, c,r) m
2. m
1,- 2F + 2m 2}, if > 4
cm
This completes the proof.
C PROOF OF THEOREM? 2
Denote bys (st)T_ probability set with Σ 8t 1. We define following functions associated with the probability set s:
P
fr
FT+1 = arg min
Blg (f(x - -T++),
JEHx
j=l
f
JTJ
arg min JEHK
Σ
s sclo ) - -3²),
j=l
t=1
F
P
arg min f
sy'o
Σ fz(z))) + Σ
T alllili,
1=1
f Dp=
t=1
5=1
j=1
and
p
T
B
f0
arg min E J,H,
Σ S1
-E((E))+000))
j=1
3m1
Correspondingly, denote by ET*(() (f) =t= Et= (s -yy)-
22

---

Published as a conference paper at ICLR 2022
Proposition 5. Let Assumptions 1-2 be true. For any f Hz, there holds
i1fi, llkk, + 2M
|lf f 111 S Ei Eg E3 + 87 (ft) +
where
Ej = =87+1 - ET- (f+t)
p
Eg = X1( 87- (f+)) - Σ
Ty ll llks,
j=l
and
P E3 ET+1 - (T+1(f)) Σ
Proof. According to Lemma 2, for any Hzz make following error decomposition
|l 1 11
(f+i(f) ET+I(f")
Me
(f) E7+1(5") +
0c
x (E)+1() - E77(())
(1) - ++7)))))))))))
-
E1
n Σ Ts
Me
oc
< E + (1) + x2() 07"0"" Σ G Ik) + 8 (1²) +
Ts 1 173
Me
+ - Ef++ +
T [1,,}
oe
j=l
$ E 177 (1) + xR(S) (
[k, + + v
1=1
E2
(T+1(G - -7++1(") + Σ
T 7
Me
1=1
2M
E. E2 (r- - + +
Ty 3ll + (f7*) +
7 Gl
ac
R8
<
E, + E2 Es ( + n Σ Tyl 315 +
2M
oc
2Mc oc
< E, + E2 + Eg +
min
{ (f) + x( + Σ T 12 +
SER2
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 {Ze 1.0 -oo with respect both the loss function f, and the hypothesis set Hz.
23

---

Published as a conference paper at ICLR 2022
Definition 4. The discrepancy describes the discrepancy between target distribution and the distribu- tion of the sample. For any f E H2, the discrepancy measure with respect to Huber loss is defined as
T
disc(s) = sup SEHz
[Ba - -
-
- v')..
t=1
Lemma 7. For any 8 > O. with probability at least 1 - 6, the following inequality holds for all fe E Hz and all a > 0
T
Bl. (f +1 ) 1-
o (f('))-)) (ssssssssss V (o 1
8log 8"
1=1
where G = (f,(f,2) : E Hzz and Me, is a positive constant such that &o (f(x)) - y²) S Mr. for any E Hz and t E Z.
Proof of Theorem 2 The proofs of bounding errors E,, Eg and Eg proceeds similarly to the proof of Theorem and are omitted for brevity.
According to Lemma we have
1 E S disc(s) + ||s||2 + 6M,, v 4n log T' + Mr, ||s/|11// 8 log 5
where
R = |||| [32o[1 + 4V222122//gg/ 2/2 (eT))
Moreover, from Proposition
4 we can get
E3 S (p 1)
(177s) 1K, 12.
Similarly, to bound the hypothesis error
Eg = 1 + -2(ff) Eq - (f/) - Σ 12c
1=1
we first give the properties of f and then them to bridge f and fo- Lemma 8. For all j 1,,,,,,,,,, .. there holds
2/po||s|12
E2
Proof. Recall the represent theorem which ensures that
T
ofj ER.
j=1t=l
For notation simplicity, denote RTp, We then deduce that
a7 n
(a y...
n,5y Rº and a. = (","))))."{{""))))
arg min
Σ soloy' Σ
(( + Σ (aaa)' Kya5},
TERTY
j=1
where K, = (K | TXT, We get
p Σ Sv Σ
K
K2 7. ) 17, k0,".
j=1
24

---

Published as a conference paper at ICLR 2022
It is easy to deduce that
p
P Σ
(s Σ
K Koraa))
Then we obtain that, for any j = I....P. P.
2011s112
Ty 1lla3 2l
Σ (e
Σ
j=l
This completes the proof.
Under Assumptions there holds
P
E2 + xN( - -
-7 >
11
11kt
5=1
= E7(8) + 1M( E7(};) (((
- > ==1
1 + 21((f)))
< xn(ft) = 12 2
21por|ss12
By combining the above results, have with confidence 1-6 8
11 1' M1
disc(s) + min JEHz
(((S)}
+0 (1g))))\[]ss|] + ||s||0 + -j >1) - + - + 11);
where C is positive constant independently of T, A, n, 6, and p. By taking T X T-1
and 1) T 4, we completes the proof.
D PROOF OF THEOREM
According to the definition of disc(T (T + 1), we have
T EC, (f(x - -y Σ
disc(s)
sup fERu
syELo sEE((f((z) - -y'}
1=1
T
q. Blo -²² - Σ
sy El (f(x²) - -v-1 + disc(q").
sup JERz
Furthermore, it is easy deduce that
T
T
sup FEHz t=1
Σ (g (f)E)
- -
sup
Σ
- (f)ec( - -²)
zEHz tm1
sup Σ st) - -/¹ - - - JEHz t=1
Proof of Theorem According the proof of Theorem in Kuznetsov & Mohri (2020a), have
T
sup JEHz
(gi -x (
+6Mr 4T (or 1q'-sl/2 8l0g
Blog
VT
25

---

Published as a conference paper at ICLR 2022
By combining above results, we have with confidence - 8
f - - Li ')
T
< disc(q") + sup
Σ
(97 - (1)o - y¹) + min { (f) + x
SEH2 t=1
JeHz
+Clog(1/6)(1l" - 2 + ||||| + |||||20 + n An + 2 2 -- + n).
By taking = T X = T and rl T we complete the proof.
Algorithm 1: Optimizationn procedure for adaptive SpHAM
Input: Data ((**, yy )))))== Max-Iter Z Z. Mercer kernel Kj,j 1,.... P with bandwidth d,
Weights n. == =1,...,p, q.
Initialization: Lipschitz constant L, s..
Step A: Computing weights 8:
for 2 1,... Z do
1. Compute A=-1 via DC-programming (or gradient descent method);
Update s* via 12
Output: $ s².
Step B: Computing f::
for z = I,.... Z do
1): Compute a² = (3*) via 14
1+ 11m
2): MM+1
2
3): g5+1 = ax
m mi+t1
(a - 1=
Output: a. oz;
Prediction function: fm =
-Ep=1 E Bi Ka(s,
Variable selection: (i: :b:10e112 v,j I..,pppp p}.
E EXPERIMENT OPTIMIZATION
The optimization problem 10 reduces to problem problem in Step A:
when taking s* = + Recall the optimization
ò = arg mine sup JEHz
E(n - ())c(((())) 3¹) |1a s112 + 2 1
The above optimization problem can be equivalently rewritten as a common type bilevel optimiza- tion problem (Colson et al. 2007), i.e.
Outer problem: ming Σi (gt (f((m)) ve) + A1llq" s|13 //2
Inner problem: fa = arg san, ΣT (gf -t) st)oo(((ar) yy.)
where the outer (min) problem parameterized by is nested within the inner (max) problem. The outer problem can be solved by standard gradient descent, where each step, we need to optimize the inner (max) problem with last updated
We denote by the iteration time. Then for k + 1-th iteration, we have the following gradient update rule
s k+1 = = - - -\i(q*
-s 22
(12)
where y is learning rate,
P Σ
4
-
ang K2 K3 x)) - ..., lo
diy K K, -
-y') R
t=1
1=1
26

---

Published as a conference paper at ICLR 2022
and
a= = (afi..... of: ofp Tpp) R2
is obtained by the following weighted optimization problem
arg maxx
Σ (i st) E1
Omy (2 24) v²).
m=1j=1
Note that the inner problem is subjected to f E Hk. i.e., ||0||2 - Σラー Tyllاl S XI min Pi 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 & -orm regularizer, and solve this regularized problem by standard gradient method. After obtaining the solution 6, we turn to solve the following weighted optimization problem Step B:
2M
arg min RyER" J=.,..P
Σ(K))a)) f=1
T [||05l12a
(13)
We can see that this problem contains non-smooth function E ع|ريyاl2 which makes the standard gradient descent method inapplicable. conquer this challenge, we leverage fast iterative algorithm (BISTA)(Beck & Teboulle 2009). Our optimization problem
becomes
a = ( = arg min {i Σo
Bllalll + 2
L lla -aak-1
Σ t=1
T
St V -
L
a
(11 |15l112
B2
(14)
15iSP
where
8i a
Σ
8 Vay lol - -KKY)(t-))
T.
and
(2(KK)a ,²)(K
1K'Ya - 88 <0
Vaylo 6o (K*)aa
20
(K*)a - 2220 20
-20
(K*)a V -0
and L = max + Ex "K is the lipschitz constant + ΣT Vloly' -KKY'* Denote Proj(s) projection of and v > 0 the threshold value. To focus the weights of samples which are really useful for forecasting, we consider Proj(s) )) as a box projection such that 8t > qt, Vt 1 L..... p). Finally, the optimization procedure for adaptive SpHAM can be summarized in Algorithm Note that if we only run Step B with = = I', we can further obtain the optimization procedure for SpHAM.
Remark 3. The computational complexity of Algorithm I 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 D(DC(T,p)) Then, the computational complexity of Step is O OZDCC(T, .p) + ZT) and the computational complexity of Step B is O pZT'pp). Thus, the total computational complexity of Algorithm I is O(ZO(DC) pT For large scale data, can further speed up Algorithm 1 by random Fourier features technique RRahmm and Recht 2007], which is leaved for future work. This has been carefully discussed in the revised manascript.
F ADDITIONAL EXPERIMENT
F.1 EVALUATION ON BENCHMARK DATA
We test our algorithm on nonlinear dataset from CauseMe. The syper-parameted selection is the same the one in Section& 4. The results in Table 3 verify the effectiveness of our method.
27

---

Published as a conference paper at ICLR 2022
Table 3: The ASE CauseMe data (p refers the dimension of features).
Methods
(p = 3, = 300)
(p = 5, 300)
LSTM
0.7681
0.9230
TS 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: // archive ics uci wwweeeu/mm/..... /ataasttsstssss Air wiir+gge (ullttuuaityyy to test model's ability to detect the Granger causality. This dataset includes 9358 hourly air quality data an Italian city, collected from March 2004 February 2005. The details of the dataset can be obtained 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 (Oa) and nitrogen dioxide (NO2) which is validated in Kalisa iet al.][[201 (2018). Mwaniki et (2014) confirm the relationship between relative humidity (RH) and nitrogen dioxide (NO2). Yan et (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 al. (019); Liu
et al. 2018), there is no any work for interpretable prediction with Granger causal network. CMES data provided in The Richardson and Cane List thttp://wee
srl ca tech.. edu/ACE
ASC//DAA lvee 1 33/cmettablee.... htm). From this link, we collect 152 ICMES observations from 1996 to 2016. The features of CMEs provided in SOHO LASCO CME Catalog (https) c/caww www....c...oooooooo nasa gow/ lsst). In-situ solar wind parameters can be downloaded
from OMNIWeb Plus (https:/ http:///mmwwwbbbbbm wwfc. waffc.nssa........m.... gov/). A total 9 features are chosen as imput, 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 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 al. (2018) state that they have 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.
28

---

Published as a conference paper at ICLR 2022
C6H6
NOx
T
RH
AH
03
NO2
NMHC
CO
Figure 3: The Granger causal network on Air Quality dataset.
29

---

Published as a conference paper at ICLR 2022
SND Speed Final
Kinetic Energy
SND Speed 20RS
MPA
CMEs Arrival Time
Mass
Linear Speed
CPA
ANG Width
Figure The Granger causal network on CMEs dataset.
30


Reviews:
[{"review_id": "jhdAW1bSFHM", "reviewer": "Reviewer_3C4P", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Overall, I thought this was a solid paper.  The SpHAM model is not that novel, but the analysis of the model appeared novel, correct and of interest to those working with real world data (e.g. financial data) that is non i.i.d., non-stationary and heavy tailed for which \"practical\" theoretical findings are not that common.\n\nPros:\n-I think this was a well written paper. The development was clean and logical and the authors nicely illustrated their contributions.\n-I think Theorem 1/Corollary 1 is a tidy result.  Expressing the rate of convergence via a moment condition, c, and the huber threshold param provides some nice intuition for situations where the method will/will not be effective.\n-The experiments nicely demonstrated the utility of the SpHAM approach for heavy tailed, stationary data and the adaptive SpHAM for non-stationary data.\n\nCons:\n-How restrictive is assumption 3 in practice? Since we don't know the true f*, it that condition testable/verifiable?\n-How do we know that a q being a uniform distribution over the most recent observations satisfies Assumption 4?  Seems like if l is to large its not likely to hold, but if l is too small then we can't capture any of the non-stationarity.\n-Can you provide a bit more color on how Step A/eq (9) is solved? \n-Can you comment on the scalability/space, time complexity of algorithm 1?  Could the method be scaled to longer time series and do you envision an incremental/online version of the procedure? \n-How does performance change if t/400 is changed in experiment B?  i.e. as the drift to noise scale ratio is changed.\n", "summary_of_the_paper": "The authors propose, theoretically analyze and empirically evaluate a sparse additive huber additive model (SpHAM), which is an additive model where each component function in the model is a kernel function that is fitted by minimizing a huber huber loss function with a functional sparsity inducing penalty term.  Statistical properties when estimating this model on stationary data with non gaussian noise and also on non-stationary data are developed. Finally, the model and fitting routines are demonstrated empirically on 2 synthetic and 1 real world time series.    ", "summary_of_the_review": "Overall, I think this was a good paper with some nice theoretical developments and an easily understandable and interpretable model.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "ja6x_T9d9fY", "reviewer": "Reviewer_KeVD", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "- Assumption 1 on boundedness of $|Y^t|$ is very strong and may not hold even for simple distributions such as Gaussian. The authors are encouraged to relax this assumption, otherwise the results of the manuscript will be limited and not general enough.\n\n- Authors need to either verify or give references to why certain kernels would satisfy Assumption 2.\n\n- The authors claim in the abstract that `\"... we propose an adaptive sparse Huber additive model for robust\nforecasting and inference (e.g., Granger causal discovery) …\". However, (a) there is no results about inference in the manuscript which would quantify the uncertainty of estimation; (b) The Granger causal network is only mentioned in the appendix for a real data set while it is not discussed how this network is estimated using the proposed algorithm and whether the estimated network makes sense scientifically. The authors need to add results related to inference, or at least modify the abstract accordingly.\n\n- The upper bound in Theorem 2 is unusual. It is appropriate for the authors to interpret the results of this theorem. Also, are the results sharp?\n\n\n- There are several hyper-parameters in the proposed algorithm, and as mentioned at top of page 9, the authors did not attempt to find data-driven methods to select such parameters. In is not clear whether the proposed methodology would still work well under different simulation settings using the same fixed hyper-parameters. A better treatment of hyper-parameter selection is needed for the method to become practical.\n\n- Experiments in Section 4.1 are repeated only 10 times which is quite low. These studies need to be replicated more, say 100 times.\n\n- The code for reproducing empirical results is not provided.\n\n\n\n\n- The sentence above Definition 2 in page 7 needs rewriting.\n\n- Statement of Theorem 3 needs rewriting: \"Let Assumption 4 the conditions in Theorem 2 be true\".\n\n\n- Typo at the beginning of section 4.2 on page 9: \"We test our algorithm on nonlinear Nonlinear-VAR dataset\".\n\n", "summary_of_the_paper": "In this manuscript, the authors proposed an adaptive sparse Huber additive model and provide an algorithm to handle non-parametric estimation of the underlying unknown function. Theoretical guarantees are provided under certain conditions and some empirical comparison with competing methods are provided.", "summary_of_the_review": "- A more comprehensive discussion on Assumptions are needed.\n- More interpretations about the theoretical results need to be provided.\n- Some claims in the abstract need to be modified.\n- A better treatment of hyper-parameter selection is required.\n- Some typos need to be fixed.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "J9TgncUEMX-", "reviewer": "Reviewer_es4A", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The novelty of this ms is mainly theoretical. It introduces a non-Gaussian noise to the additive model can come up with the error bounds. The experiments are not comprehensive. It only compares with two other auto regressive methods. I will be interested to see its performance in non-autoregressive setting, i.e. simply sum several non-linear functions and Gaussian / non-Gaussian noise and see if the ground truth could be recovered. Also, the real data experiment (F.1) result is only a figure without any interpretation. It will be nice to see if this method provides new insights compared with previous methods.", "summary_of_the_paper": "This ms tackles a timer series problem where each observation is the sum of p hidden functions and a noise (termed as innovation in ms). Traditionally the noise is Gaussian, but this ms makes the noise into heavy tailed distribution such as t distribution. The L2 loss function of the original problem is then replaced by so called Huber loss, which introduces a cutoff to the deviations and is less affected by outliers than L2 loss. Based on this new loss function, this ms sets up the sparse Huber additive model and derive the bounds of the learned hidden functions, under stationary processes. In case of non-stationary processes, the ms introduces different weights to each time points, which can be used to calculate a discrepancy that reflects the non-stationarity. With different weights to different time points, the ms come up with new error bounds. Inference algorithms are provided in the end.", "summary_of_the_review": "My general feeling is that there is theoretical contribution but the experiments could be improved.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "C-nx-fDVbE7", "reviewer": "Reviewer_4Rek", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths: \n1. The proposed time series model is both robust and interpretable. \n2. Some novel performance lower bounds for the proposed model. \n\nWeaknesses:\n1. Although the simple linear additive model is well interpretable, its modeling capacity is quite weak on the other hand. As a remedy not destroying your optimization framework, you may consider adding some interacting terms, like Y2*Y3. \n2. The proposed robust time series prediction model with the three theorems constitutes the main contribution of the paper. But to me, the meanings of the lower bound (namely each part of the right-hand side) are hard to interpret. Are they tight bounds in general? Besides, are $\\hat{f}$ solutions of the corresponding optimization problems? \n3. It is well known that robust Huber criterion function can handle non-Gaussian innovations/errors, but it remains unclear in the paper how is it used in the proposed model to handle non-i.i.d. data, which is regarded as one major contribution if I were right.\n4. Could the authors comment on the impact of the kernel selection? I see in the experiments only a Gaussian kernel (which is a stationary kernel) is demonstrated. Don't you need a non-stationary kernel instead?  \n5. Experimental validation is somewhat insufficient. The selected examples (A&B) well match the linear additive model assumption. It would be interesting to see some model mismatching examples and the performance thereof. Besides, in both examples, the innovations/errors are i.i.d., right? \n6. As benchmarks for comparison, only a simple AR model is adopted. The authors should compare the proposed model with more state-of-the-art models, such as LSTM, Gaussian process (GP) regression, and even some latest attention-based models. Actually, your model is sort of similar to a GP model when replacing Huber loss back to L2 loss. \n7. Code is not made publicly available. \n8. Many typos and grammatical errors. \n   ", "summary_of_the_paper": "The authors proposed a robust Huber additive model for non-stationary time series prediction. They combine the idea of robust Huber regression (against non-Gaussian innovations) with a linear additive model formulated with the aid of the representer theorem in RKHS; sparsity is imposed on the weights, and the overall weights optimization problem can be solved efficiently through the classic FISTA algorithm due to the special problem structure. The major contributions lie in a robust (due to the Huber loss function) and interpretable (due to the use of linear additive model) time series prediction model as well as some theoretical supports on the performance. ", "summary_of_the_review": "To conclude, I find the paper interesting and easy to follow, the technical contributions are fine but obviously need further refinement. The literature review is insufficient, and some SOTA methods, such as LSTM, GPR, etc., are missing. The section of experimental results can be improved. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes an adaptive sparse Huber additive model for for forecasting non-stationary time series. The prior work has considered similar models for Gaussian innovations which is overly restrictive for a variety of applications such as finance. The results are supported both by theory and experiments. The results are novel and are of interest to ICLR and machine learning communities in general.
Decision: Accept (Poster)