Sparse Deep Learning for Time Series Data: Theory and Applications
Authors: Mingxuan Zhang
Yan Sun
Faming Liang
Venue: NeurIPS 2023 poster (2023)
Keywords: Sparse Deep Learning
Uncertainty Quantification
Model Compression
Variable Selection
Dependent Data

Abstract: Sparse deep learning has become a popular technique for improving the performance of deep neural networks in areas such as uncertainty quantification, variable selection, and large-scale network compression. However, most existing research has focused on problems where the observations are independent and identically distributed (i.i.d.), and there has been little work on the problems where the observations are dependent, such as time series data and sequential data in natural language processing. This paper aims to address this gap by studying the theory for sparse deep learning with dependent data. We show that sparse recurrent neural networks (RNNs) can be consistently estimated, and their predictions are asymptotically normally distributed under appropriate assumptions, enabling the prediction uncertainty to be correctly quantified. Our numerical results show that sparse deep learning outperforms state-of-the-art methods, such as conformal predictions, in prediction uncertainty quantification for time series data. Furthermore, our results indicate that the proposed method can consistently identify the autoregressive order for time series data and outperform existing methods in large-scale model compression. Our proposed method has important practical implications in fields such as finance, healthcare, and energy, where both accurate point estimates and prediction uncertainty quantification are of concern.

Sparse Deep Learning for Time Series Data: Theory and Applications
Mingxuan Zhang Department of Statistics Purdue University
zhan36920purdue- edu
Yan Sun
Faming Liang
Department of Biostatistics, Epidemiology, and Informatics University of Pennsylvania
Department of Statistics
fmliang@purdue. edu
yan. www.sennmmddiicinnee.m www.eiin.nnn..nnnn.mm..mmmmmmmne upenn. edu
Abstract
Sparse deep learning has become a popular technique for improving the performance of deep neural networks in areas such as uncertainty quantifi cation, variable selection, and large-scale network compression. However, most existing research has focused on problems where the observations are independent and identically distributed (i.i.d.,, and there has been little work on the problems where the observations are dependent, such as time scries data and soquential data in natural language processing. This paper aims to address this gap by studying the theory for sparse deep learning with dependent dala. We show that sparse recurrent neural nelworks (RNNs) can be consistently estimated, and their predictions are asymptotically nor- mally distributed under appropriate assumpt ions, enabling the prediction uncertainty to be correctly quantified. Our numerical results show that sparse deep learning outperforms state-of-the-art methods, such as confor- mal predictions, in prediction uncertainty quantification for time series data. Furthermore, our results indicate that the proposed method can consistently identify the autorcgressive order for time series data and outperform ex- isting mothods in large-scale model compression. Our proposed method has important practical implications in fields such as finance, healthcare, and energy, where both accurate point estimates and prediction uncertainty quantification of concern.
Introduction
Over the past decade, deep learning has experienced unparalleled triumphs across multitude of domains, such time scries forccasting [
g 4 5.
natural language processing
and computer vision [ 0 However, challenges like generalization and miscalibration [0]]
persist, posing potential risks critical applications like medical diagnosis and aulonomous vehicles.
In order enhance the performance of deep neural networks (DNNs), significant research efforts have been dedicated to exploring optimization methods and the loss surface of the DNNs, see, c.g., [1 12
15 16. which have aimed to expodite and direct the con- vergence of DNNs towards regions that exhibit strong generalization capabilitics. While these investigations are valuable, offectively addressing both the challenges of generalization
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

and miscalibration require additional and perhaps more essential aspects: consistent esti- mation of the underlying input-output mapping and complete knowledge of the asymptotic distribution of predictions. As highly effective method that addresses both challenges, sparse deep learning has been extensively studied, see, e.g., T7 it is important to note that all the studies have been conducted under the assumption of independently and identically distributed (i.i.d.) data. However, in practice, we frequently encounter situations where the data exhibits dependence, such as time series data.
18
19
20
21 Nevertheless,
The primary objective of this paper to address this gap by establishing a theorctical foundation for sparse deep learning with time series data. Specifically, we lay the foundation within the Bayesian framework. For RNNs, by letting their parameters be subject to a mixture Gaussian prior, we establish posterior consistency, structure selection consistency, imput- output mapping estimation consistency, and asymptotic normality of predicted values. We validate our theory through numerical experiments on both synthetic and real-world datasets. Our approach outperforms existing state-of-the-art methods in uncertainty quantification and model compression, highlighting its potential for practical applications where both accurate point prediction and prediction uncertainty quantification are of concern.
Additionally, we would elaborate the contribution of this paper in a broader context of statistical modeling. As discussed in [22, two distinct cultures exist for statistical modeling: the 'data modeling culture' and the 'algorithmiic modeling culture'. The former focuses on simple generative models that explain the data, potentially lacking a consistent estimate of the true data-generating mechanism due to the model's inherent simplicity. The latter, on the other hand, aims to find models that can predict the data regardless of complexity. Our proposed method occupies middle ground between these two cultures. It seeks to identify a parsimonious model within the realm of complex models, while also ensuring a consistent estimation of the true Data-generating mechanism. From this perspective, this work and related ones, e.g., [7 [ 19 20 221, represent a hybridization of the 'algorithmic modeling culture' and the 'data modeling culture', which holds the potential to expedile advancements in modern data science.
2 Related Works
Sparse deep learning. Thcoretical investigations have been conducted the approxi- mation power of sparse DNNs across different classes of functions 23 24 Recently 20 has made notable progress by integrating sparse DNNs into the framework of statistical modeling, which offers a fundamentally distinct neural network approximation theory. Unlike traditional theories that lack data involvement and allow connection weights to assume values in an unbounded space to achieve arbitrarily small approximation errors with small networks 24. their theory 200 links network approximation error, network size, and weight bounds to the training sample size. They show that a sparse DNN of size O(n/ log(n)) can effectively approximate various types of functions, such as affine and piecewise smooth functions, as n + 00, where n denotes the training sample size. Additionally, sparse DNNs exhibit several advantagcous thcorctical guarantoes, such as improved interpretability, enabling the consistent identification of relevant variables for high-dimensional nonlinear systems. Building upon this foundation, [2] establishes the asymptotic normality of connoction weights and predictions, enabling valid statistical inference for predicting uncertainties. This work extends the sparse deep learning theory of 20 21 from the case i.i.d data to the case of time series data.
Uncertainty quantification. Conformal Prediction (CP) has emerged as a prominent technique for generating prediction intervals, particularly for black-box models like neural networks. A key advantage of CP is its capability to provide valid prediction intervals for any data distribution, even with finite samples, provided the data meets the condition of exchangeability 26 27. While i.i.d. data casily satisfics this condition, dependent data, such as time series, often docsn't. Researchers have extended CP to handle time scries data by relying on properties like strong mixing and ergodicity 28 29 In a recent work, B0 introduced a random swapping mechanism to address potentially non-exchangeable data, allowing conformal prediction to be applied on top of a model trained with weighted samples. The main focus of this approach was to provide theoretical basis for the
2

---

differences observed in the coverage rate of the proposed method. Another recent study by BT took deep dive into the Adaptive Conformal Inference (ACI) 32 leading to the development of the Aggregation Adaptive Conformal Inference (AgACT) method. In situations where a dataset contains a group of similar and independent time series, treating each time series as separate observation, applying a CP method becomes straightforward [3). For a comprehensive tutorial on CP methods, one can refer to [3]. Beyond CP,
other approaches for addressing uncertainty quantification in time series datasets include multi-horizon probabilistic forecasting 35 methods based on dropout 30 and recursive Bayesian approaches 37
3 Sparse Deep Learning for Time Series Data: Theory
Let Dn = -V1,, (y1,....,,b Bn denote time series sequence, where yi E R. Let (2, F, PO) be the probability space of Dn, and lct 0g = sup[[ (y, E A, yk+s E - P" P E A) (Uk+s E
B)] A, B F.j E ] be the k-th order a-mixing coefficient.
Assumption 3.1. The time series D, is (strictly) stationary and a-mixing with an expo- nentially decaying mixing coefficient and follows an autoregressive model of order I
y 1² 11--1 u.) + 7/6
(1)
where μ² is a non-linear function, Vi- 1--1:i-1 (y: Wi:-...... Ii-,), u, contains optional exogenous variables, and Ty i.d.d ~ o20, with x² being assumed to be a constant.
Remark 3.2. Similar assumptions are commonly adopted establish asymptotic properties of stochastic processes [8 39, 40
42 29). For example, the asymptotic normality of the maximum likelihood estimator (MLE) can be established under the assumption that the time series is strictly stationary and ergodic, provided that the model size is fixed [38]] A posterior contraction rate of the autoregressive (AR) model can be obtained by assuming it is a-mixing with Ex= Cg
41
1-2/s < oc for some 8 > 2 which is implied by an exponentially decaying mixing coefficient 12] For stochastic processes that are strictly stationary and f-mixxng, results such as uniform laws of large numbers and convergence rates of the empirical processes
39 20 can also be obtained.
Remark 3.3. Extending the results of 20 21 to the case that the dataset includes a set of i.i.d. time scries, with each time series regarded as an individual observation, is straightforward, this is because all observations are independent and have the same distribution.
3.1 Posterior Consistency
Both the MLP and RNN can be used to approximate 1²" as defined in 11, and for simplicity, we do not explicitly denote the exogenous variables u, unless it is necessary. For the MLP, we can formulate as a regression problem, where the input is æ; = Wi-1u-R: for some I < Ri < n, and the corresponding output is yi, then the dataset Dn can be expressed as {( 1 Detailed settings and results for the MLP are given in Appendix B.3 In what follows, we will focus on the RNN, which serves as an extension of the previous studies. For the RNN, we can rewrite the training dataset as Dn = [M h for some I R < Mi < n, i.e., split the entire sequence into a set of shorter sequences, where R₁ denotes an upper bound for the exact AR order I, and Mi denotes the length of these shorter sequences (see Figure[] We assume R₃ is known but not since, in practice, it is unlikely that we know the exact order I.
For simplicity of notations, we do not distinguish between weights and biases of the RNN. In this paper, the presence of the subscript n in the notation of a variable indicates its potential to increase with the sample size n. To define an RNN with Hn - 1 hidden layers, for E {1, 21, Hn); we let w and Ln denote, respectively, the nonlinear activation function and the number of hidden neurons at layer h. We set LH,, = and Lo Pa: where Pn. denotes a generic input dimension. Because of the existence of hidden states from the past, the input x; can contain only Wi-1 3i-1:i-r for some > 1. Let wh e R La-1 and h E (Rlax x denote the weight matrices at layer h. With these notations, the output of the step i of an RNN model can be expressed
p( (z/_ Hr -n- ,B)- w Hh Hn 11. (.... ...... . /wo' x; + w² 21 2/111
(2)
3

---

R-
Figure A multi-layer RNN with an input window size of k. We restrict the use of the RNN's outputs until the hidden states have accumulated sufficient quantity of past information to ensure accurate predictions.
where x = [o"["=x]] + v" =ا=, denotes the hidden state of layer h at step with 2 = 0:
and B the collection of all weights, consisting Kn = - (Ln x Ln- 1n-1)+++ - (L2) (1 Ln elements. To represent the structure for a sparse RNN, we introduce an indicator variable for each weight in B. Let 7 = (7g E {0, 13 j 1,...: Kn); which specifies the structure of sparse RNN. To include information concise, we redenote μ(z; f2t-11111110 1- only on (3 and up to -i
the network structure and keep the notation
H., -0-------00 (3) by
H.,-1 depends 3.7), as f2 f2,-1/m-1 -
Posterior consistency is an essential concept in Bayesian statistics, which forms the basis of Bayesian inference. While posterior consistency generally holds for low-dimensional problems, establishing it becomes challenging in high-dimensional scenarios. In such cases, the dimensionality often surpasses the sample size, and if the prior is not appropriately elicited, prior information can overpower the data information, leading to posterior inconsistency.
Following ie.,
20
21
let cach connection weight be subject to mixture Gaussian prior,
~ N00, + (1-1,,11))) (1,, jE (1,2....,KKK K,),,
by integrating out the structure information 7: where An (0, the mixture proportion, ogn is lypically set to a very small number, while of is relatively large. Visualizations of the mixture Gaussian priors for different Xn, ob and of,n are given in the Appendix E
We assume n can be well approximated by sparse RNN given enough past information, and refer to this sparse RNN as the true RNN model. To be more specific, we define the true RNN model
arg min (8 (8,-)== (B.7)E0... |ama My 1.-))) "" ) bg(0) San
hy|
(4)
where Gn == g(Cu, C1, e, Pn. L1, Le...., LH. denotes the space of all valid networks that satisfy the Assumption 3.4 for the given values of Hny Pn; and L4 's, and wa some sequence converging to 0 as n + xo. For any given RNN (B,7) 7), ) the error (n*(.) p(: (,(,33,7)) 7)1 can be decomposed as the approximation error (.((.) (a(;,8", x²)l and the estimation error |n(=/9", - (3,)). The former is bounded by wa, and the order of the latter will be given in Theorem 3.9 For the sparse RNN, we make the following assumptions:
Assumption 3.4. The true sparse RNN model (3",>") satisfies the following conditions:
+s, The network structure satisfies: T, log n + log L sn log < Cyn'-s, where 0 e 1 is small constant, Tn |r*1 denotes the connectivity of x", L max Ly denotes the maximum hidden state dimension, and sn denotes the input dimension of y
The network weights are polynomially bounded: || En, where E, n2 for some constant C₁ > 0.
Remark 3.5. Assumption 3.4 identical assumption A.2 of 20 which limits the con- nectivity of the true RNN model to be of o(n!-t) for some 0 1. Then, as implied by

---

an RNN of size O [ogn)) has been large enough for modeling a time series sequence length n. Refer 20 for discussions on the universal approximation ability of the neural netw under this assumption; the universal approximation ability can st.ill hold for many classes of functions, such as affine function, piecewise smooth function, and bounded a- a-Hölder smooth function.
Lemma 3.10
Remark 3.6. The existence of the sparse RNN model stems from Lemma 4.1 of [0 which, through the trick of independent block sequence construction, shows that many properties of the i.i.d processes can be extended to mixing processes. While the lemma was proven for the case of B-mixing, the author did mention her doubts about its applicability to a-mixing. Therefore, at least for the sequences of 3-mixing, which implies a-mixing, the non-empty of the sparse RNNN set in 4 can be guaranteed for many classes of functions as mentioned in Remark 3.5
Assumption 3.7. The activation function sa is bounded for h (e.g., sigmoid and tanh), and is Lipschitz continuous with a Lipschitz constant of for 2 s h S Hn (e.g., ReLU, sigmoid and tanh).
except that we require the activation function Remark 3.8. Assumption 3.7 mirrors 20 21,
for the first layer to be bounded. This extra assumption can be viewed as a replacement for the boundedness assumption for the input variables of a conventional DNN.
Let d(p1.P2) denote the integrated Hellinger distance between two conditional densitics Pı (v/æz and pz(y/z). Let x(-(D,) be the posterior probability of an event. Theorem 3.9 establishes posterior consistency for the RNN model with the mixture Gaussian prior
Theorem 3.9. Suppose Assumptions 3.1 3.4. and 3.7 hold. If the misture Gaussian prior satisfies the conditions Xn
O(1/[M, Kn 2M, (L.N)""(Lp)))" for some con- I 01, S no for some constant a, and 0on I then there ea- ists an an error sequence en = 0( 0(=))+ such that lim en = 0 and lim ne, 0o, and the posterior distribution satisfies
stant T T> O, Hn, logn + log L.11/2
min(1/[M(" VnK,n / H., )2M.)
1/[M/ Vnkkn (
π add 2 Cca, Dn,) Op- (c (c=
(5)
for sufficiently large n and C' > 0, where 22 = Hn logn log &+ pnl/n, Pp- denotes the underbying true data distribution, and pæ denotes the data distribution reconstructed by the Bayesian RNN based on ils posterior samples.
3.2 Uncertainty Quantification with Sparse RNNs
As mentioned previously, posterior consistency forms the basis for Bayesian inference with the RNN model. Based on Thcorem 3.9 we further establish structure selection consistency and asymptotic normality of connection weights and predictions for the sparse RNN. In particular, the asymptotic normality of predictions enables the prediction intervals with correct coverage rates to be constructed.
Structure Sclection Consistency It is known that the neural network often suffers from a non-identifiability issuc due to the symmetry of its structure. For instance, one can permute certain hidden nodes or simultaneously change the signs or scales of certain weights while keeping the output of the neural network invariant. To address this issue, we follow 20 to define an equivalent class of RNNs, denoted by θ which is set of RNNs such that any possible RNN for the problem can be represented by one and only one RNN in θ via appropriate weight transformations. Let v(8, 7) θ denote an operator that maps any RNN to O To serve the purpose of structure selection in the space θ we consider the marginal posterior inclusion probability (MIPP) approach. Formally, for each connection weight i = 1,..., Kn, we define its MIPP as qi = En 4 (i (.) where ())(B.7) is the indicator of i in v(8, 7). The MIPP approach selects the connections whose MIPPs exceed threshold q. Let 9q = {i: qa 4, 1..... Kn} denote an estimator K.)....). Let A(cx) = {8 d(ps. P < Gn) and
-
of 7. = {i Gip("', 1, p(Gn)) = maxi SA(e.) Ey lemp((.)) - equ("" xx) (o (718, (.) (.) )dg, which measures
5

---

the structure difference on A(e) for the true RNN from those sampled from the posterior. Then we have the following Lemma: Lemma 3.10. If the condilions of Theorem 3.9 are salisfied and p(en) - as n 00, then
(i) mixisisk, f1a: - -slu((8 3 as 00;
(ti) (sure screening) P(y. Cga) -1 as + 0o, for any prespecified @ E (0, 1);
(bii) (consistency) P(y. f..s) S as 00;
where 4 denotes convergence in probability.
Remark 3.11. This lemma implies that we can filter out irrelevant variables and simplify the RNN structure when appropriate. Please refer to Scction 5.2 for numerical illustration.
Asymptotic Normality of Connection Weights and Predictions The following two theorems establish the asymptotic normality of ((3) and predictions, where x(3) denotes a transformation of 3 which is invariant with respect to ((.88,7) while minimizing |0( 1l
We follow the same definition of asymptotic normality as in distribution for the function g(B) is asymptotically normal with center g and variance G if, for d3 the bounded Lipschitz metric for weak convergence, and on the mapping On vn(g(B) - g²), it holds, n oc, that
19
The posterior
o N(O, ((0, +0,
(6)
in P"-probability, which we also denote (.. o on ~ N)O,
The detailed assumptions and sctups for the following two thcorems are given in Appendix
C For simplicity, we let Mj = Rg + 1, and let 4, (8) и-H Li=n++1 log(Ps 11
denote the averaged log-likelihood function. Tet (3) denote the Hessian matrix of (B)),
and let (B) and h-d (B) denote the (i, element of H,) and H₂ respectively.
Theorem 3.12. Assume the conditions of Lemma 3.10 hold with p(en) = of Kg
and
additional assumptions hold given in Appendir C then क. ((vn(((()) - 3)) ~ (O, in P -"opbobabilliyy otherwise.
00, where V = (eig), and Vej E(h'' if i, E r and 0
Theorem 3.13 establishes asymptotic normality of the sparse RNN prediction, which implies prediction consistency and forms the theoretical basis for prediction uncertainty quantification as well.
Theorem 3.13. Assume the conditions of Theorem 3.12 and additional assumptions " D,)) ~ NOO, where Σ = H pl, and H = Fl- - (B'))(("))) is the Fisher information matria.
hold given in Appendir C Then T T[Vn(u((.,)) 3) ul,,
4 Computation
In the preceding section, we established a theoretical foundation for sparse deep learning with time series data under the Bayesian framework. Building on [20], is straightforward to show that Bayesian computation can be simplified by invoking the Laplace approximation theorem the maximum posteriori (MAP) estimator. This essentially transforms the proposed Bayesian method into regularization method by interpreting the log-prior density function as penalty for the log- log-liicl function in RNN training. Conscquently. we can train the regularized RNN model using an optimization algorithm, such as SGD or Adam. To address the local-trap issue ialentially suffered by these methods, we train the regularized RNN using prior annealing algorithm [11] as described in Algorithm [ For trained RNN, we sparsify its structure by truncating the weights less than a threshold to zero and further refine the nonzero weights for attaining the MAP estimator. For algorithmic

---

specifics, refer to Appendix D Below, we outline the steps for constructing prediction intervals for one-step-ahead forecasts, where 1 9) is of one dimension, and B and fy represent the estimators of the network parameters and structure, respectively, as obtained by Algorithm
Estimate o² by 6² R4
R++1 11 alli--iii- Ry: B)11
For test point yk-1.k 1ik- R. estimate E in Theorem 3.13 by
(V ,,, -V
The corresponding (1 - a3% prediction interval is given by
1: + =a/2 (n Rτ) 62
where there are n R observations used in training, and 2a/2 denotes the upper
a/2-quantile of the standard Gaussian distribution.
For construction of multi-horizon prediction intervals, see Appendix P
Numerical Experiments
5.1 Uncertainty Quantification
As mentioned in Section 3 we will consider two types of time series datasets: the first type comprises single time scries, and the second type consists of set of time series. We will compare the performance of our method against the state-of-the-as (CP) methods for both types of datasets. We set 0.1 (the error rate) for all uncertainty quantification experiments in the paper, and so the nominal coverage level of the prediction intervals is 90%.
Conformal Prediction
5.1.1 A Single Time Series: French Electricity Spot Prices
We perform onc-step-ahead forecasts on the French electricity spot prices data from 2016 to 2020, which consists of 35,064 observations. A detailed description and visualization of this time series are given in Appendix G.I. Our goal is to predict the 24 hourly prices of the following day, given the prices up until the end of the current day. As the hourly prices exhibit distinct patterns, we fit one model per hour as in the CP baseline BI. We follow the data splitting strategy used in [T. where the first three years 2016 2019 data are used as the (initial) training set, and the prediction is made for the last year 2019 2020.
For all the methods considered, we use the same underlying neural network model: an MLP with one hidden layer of size 100 and the sigmoid activation function. More details on the training process are provided
Appendix G.II For the statc-of-the-art CP methods, EnbPI-V2 29. NEX-CP 30, ACT [2] and AgACI [T. we conduct experiments in an online fashion, where the model is trained using a sliding window of the previous three years of data (refer to Figure 4 in the Appendix). Specifically, after constructing the prediction interval for each time step in the prediction period, we add the ground truth value the training set and then rein the model with the updated training set. For ACT, we conduct experiments with various values of y and present the one that yields the best performance. In the case of AgACI, we adopt the same aggregation approach used in BII namely, the Bernstein Online Aggregation (BOA) method 44 with a gradient trick. We also report the performance of ACI with y 0 as reference. For NEX-CP, use the same weights as those employed in their time-scrics experiments. For EnbPI-V2, we tunc the number of bootstrap models and select the one that offers the best performance.
Since this time series exhibits no or minor distribution shift, our method PA trained an oflline fashion, where the model is fixed for using only the observations between 2016 and 2019, and the observations the prediction range are used only for final evaluation. That is, our method uses less data information in training compared to the baseline methods.

---

40
35
30
Average length
30
Median length
20
20
15-
5-
10:
90.0
90.5
91.0
91.5
92.0
92.5
95.0
90.0
90.3
91.0
92.5
92.0
92.5
95.0
Coverage
PA offline
online
onlina
y
EntPI
NEX-CP
Comparisons of our method and basclines on hourly spot electricity prices in France. Left: Average prediction length with vertical error bars indicating the standard deviation of the prediction interval lengths. Right: Median prediction length with vertical error bars indicating the interquartile range of the prediction interval lengths. The precise values for these metrics are provided in Table 5 of Appendix G.I
Figure 2:
The results are presented in Figure 2 which includes empirical coverage (with methods that positioned closer 90.0 being more effective), poodian/average prediction interval length, and corresponding interquartile range/standard deviation. As expected, our method is able to train and calibrate the model by using only the initial training set, i.c., the data for 2016-2019, and successfully produces faithful prediction intervals. In contrast, all CP methods produce wider prediction intervals than ours and higher coverage rates than the nominal level 90%. In addition, ACI is sensitive to the choice of 31
5.1.2 A Set of Time Series
We conduct experiments using three publicly available real-world datasets: Medical Informa- tion Mart for Intensive Care (MIMICIIIT). eleetroencophallgrappe
(EEG) data, and daily COVID-19 case numbers within the United Kingdom's 's local authority districts (COVID-19) 83 A concise overview of these datascts is presented in Table 1 Our method, denoted as
PA-RNN (since the underlying prediction model is an LSTM 45 ), compared to three benchmark methods: CF-RNN [3], MQ-RNN 35 and DP-RNN 36 where LSTM is used
for all these methods. To ensure a fair comparison, we adhere to the same model struc- ture, hyperparameters, and data processing steps specified in 88 Detailed information
regarding the three datasets and training procedures can found in Appendix G.2
The numerical results are summarized in Tables 2 Note that the baseline results for EEG and COVID-19 are directly taken from the original paper [3 We reproduce the bascline results for MIMIC-III, as the specific subset used by 33 is not publicly available. Table 2 indicates that our method consistently outperforms the basclines. In particular, Our method consistently generates shorter prediction intervals compared to the conformal baseline CF-RNN while maintaining the same or even better coverage rale as CF-RNN. Both the MQ-RNN and DP-RNN methods fail to generate prediction intervals that accurately maintain a faithful coverage rate.
where "Train size" and "Test size" indicate the numbers of training and testing sequences, respectively, and "Length" represents the length of the input sequence.
Table 1: A brief description of the datasets used in Section 5.1.2
Dataset
Train size
Test size
Length
Prediction horizon
MIMIC-III
46
2692
300
[3,47]
EEG 47
19200
19200
40
10
COVID-19
48
300
80
100
50
8

---

Table Uncertainty quantification results by different methods for the MIMIC-III, EEG, and COVID-19 data, where "Coverage" represents the joint coverage rate averaged over five different random seeds, the prediction interval length is averaged across prediction horizons and five random seeds, and the numbers in the parentheses indicate the standard deviations (across five random seeds) of the respective means.
MIMIC-III
EEG
COVID-19
Model
Coverage
PI length
Coverage
PI length
Coverage
PI length
PA-RNN
90.8%(0.7%
94.69%(0.1%% 2.35(0.26)
51.02(10.50)
90.5%(1.44%
120.93
CF-RNN
(2.0%(0.3%)
(0.4%)) 3.01(0.17)
61.86(18.02)
89.7%(..44%
733.95(582.52))
MQ-RNN
85.3%(0.2%)
2.64(0.11)
48.0%(1.8%%
21.39(2.36)
15.0%(2.6%)
136.56(63.32)
DP-RNN
(0.3%)
0.12(0.01)
(0.3%)
(.39().74)
(.0%())0%)
61.18(32.37)
5.2
Autoregressive Order Selection
In this section, we evaluate the performance of our method in selecting the autoregressive order for two synthetic autoregressive processes. The first is the non-linear autoregrossive (NLAR) process
49 50,
5I:
= -..17 + 0.85 + 0.14 2 0.31y:-3 + 0.08 + 12.80G (M--1a-7) (. 2.44C2(v:-1a-)) + ne
where 77: ((, 1) represents i.i.d. Gaussian random noises, and the functions G1 and G2 are defined as:
Gi(y--a-77
exp(-0.16(0099))--1--- + 6.68)))
G2(y)-1a-7) (1 exp 1.17 10 $00.83y, 0.53y. 0.18y; 0.38) }
The second is the exponential autoregressive process 52
y: = (0.8- 1.1 exp(-500 9i-1 + Th:
where, again, Y ~ NO,,) 1) denotes i.i.d. Gaussian random noises.
For both synthetic processes, we generate five datasets. Fach dataset consists of training, validation, and test sequences. The training sequence has 10000 samples, while the validation and test sequences each contain 1000 samples. For training, employ single-layer RNN with a hidden layer width of 1000. Further details on the experimental setting can be found in Appendix G.3
For the NLAR process, we consider two different window sizes for RNN modeling: 15 (with input as 15-15 and 1 (with input as 1). Notably, the NLAR process has an order of 7. In the case where the window size is 15, the input information sullices for RNN modeling, rendering the past. information conveyed by the hidden states redundant. However, when the window size 1, this past information becomes indispensable for the RNN. In contrast, the exponential autoregressive process has an order of 1. For all window sizes we explored, namely 1, 7, 10, 15, the input information is always sufficient for RNN modeling.
We evaluate the predictive performance using mean square prediction error (MSPE) and mean square fitting error (MSFE). The model selection performance is assessed by two metrics: the false selection rate (FSR) and the negative selection rate (NSR). We define
FSR Σ, [s,/S| and NSR = Σ, |S/š,1
where S denotes the set of true variables, and >9,181
E,=1$,1
8, represents the set of selected variables for dataset j. Furthermore, we provide the final number of nonzero connections for hidden states and the estimated autoregressive orders. The numerical results for the NLAR process are presented in Table 3 while the numcrical results for the exponential autoregressive process given in Table
Our results are promising. Specifically, when the window size equal to or exceeds the true autoregressive order, all connections associated with the hidden states are pruned, effectively converting the RNN into an MLP. Conversely, the window size is smaller than the true autoregressive order, significant number of connections from the hidden states

---

are retained. Impressively, our method accurately identifies the autoregressive order noteworthy achievement considering the inherent dependencies in time series data. Although our method produces nonzero FSR for the NLAR process, it is quite reasonable considering the relatively short time sequence and the complexity of the functions G1 and G2.
Table 3: Numerical results for the NLAR process: numbers in the parentheses are standard deviations of the respective means, '_" indicates not applicable, & means lower is better, and "Hhideen link" denotes the number of nonzero connections from the hidden states. All results are obtained from five independent runs.
Model Window size FSR 4
NSR 4
AR order #hidden link
MSPE 4
MSFE +
PA-RNN PA-RNN
0
(27
1.037 1.056(0.001) 1.017(0.008)
(0.000.))
0.23
0
7.4(0.25)
0(0)
1.020(0.010)
Table 4: Numerical results for the exponetial autoregressive process.
Model
Window size
FSR +
NSR +
AR order
#hidden link
MSPE 4
MSFE
PA-RNN
0
1(0)
0(0)
1.004(0.000
1.003 (0.)00
PA-RNN
0
1(0)
0(0)
1.006(0.005)
0.999(0.004)
PA-RNN
0
1(0)
0(0)
1.000(0.004)
(.007(0.005)
PA-RNN
0
1(0)
0(0)
1.006(0.005)
1.000(0.003)
PA-RNN
0
1(0)
0(0)
(.002(0.004)
1.002 (0.06)
PA-RNN
0
1(0)
0(0)
1.001
1.002(0.007)
5.3 RNN Model Compression
We have also applied our method to RNN model compression, achieving state-of-the-art results. Please refer to Section G.4 in the Appendix for details.
6
Discussion
This paper has established the theoretical groundwork for sparse deep learning with time series data, including posterior consistency, structure selection consistency, and asymptotic normality of predictions. Our empirical studics indicate that sparse deep lcarning can outperform current cutting-edge approaches, such as conformal predictions, in prediction uncertainty quantification. More specifically, compared to conformal methods, our method maintains the same coverage rate, if not. betler, while generaling significantly shorter prediction intervals. Furthermore, our method effectively determines the outoregressoon order for time series data and surpasses state-of-ther the-art techniques in large-scale model compression. The theory developed in this paper has included LSTM [5 as a special case, and some numerical examples have been conducted with LSTM; see Section G of the Appendix for the detail. Furthermore, there is room for refining the thcoretical study under varying mixing assumptions for time scries data, which could broaden applications of the proposed method. Also, the efficacy of the proposed mothod can potentially be further improved with the elicitation of different prior distributions.
In summarr, this paper represents a siguificant advancement in statistical inference for deep RNNs, which, through sparsing, has successfully integrated the RNNs into the framework of statistical modeling. The superiority of our method over the conformal methods shows further the criticality of consistently approximating the underlying mechanism of the data generation process in uncertainty quantification.
In terms of limitations of the proposed method, one potential concern pertains to the calculation of the inverse of the Fisher information matrix. For largo-scale problems, the sparsified model could still retain large number of non-zero parameters. In such instances, the computational feasibility of calculating the Hessian matrix might become compromised. Nonetheless, alternative avenue exists in the form of the Bayesian approach, which circumvents the matrix inversion challenge. A concise description of this Bayesian strategy can be found in 21
10

---

References
[]] Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu. Modeling long-am and short-term temporal patterns with deep neural networks. In The 41st international ACM SIGIR conference on research & development in information relrieval, pages 104, 2018.
[2] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. Deepar: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting, 36(3):1181-13) 1191, 2020.
3] Syama Sundar Rangapuram, Matthias W Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski. Deep state space models for time series forecasting. Advances in neural information processing systems, 31, 2018.
[4] Cristian Challu, Kin G Olivares, Boris N Oreshkin, Federico Garza, Max Mergenthaler, and Artur Dubrawski. N-hits: Neural hierarchical interpolation for time series forocasting. arXiv preprint ar wrXtx:2201.12880 2022.
5 Hansika Hewamalage, Christoph Bergmeir, and Kasun Bandara. Recurrent neural nolw for time series forecasting: Current status and future directions. International Journal of Forecasting, 37(1):388 :388-447, 2021.
6) Gábor Melis, Tomáš Kočiský, and Phil Blunsom. Mogrifier Istm. ar preprinl 1xx::199.001 2019.
[7] Luciano Floridi and Massimo Chiriatti. Gpt-3: Its nature, scope, limits, and conse- quences. Minds and Machines, 30:681-699 694, 2020.
8) Yuki Tatsunami and Masato Taki. Sequencer: Deep 1stm for image classification. arXiv preprint ar Xic:22000 01977, 2022.
(9] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo. Swin transformer: Hicrarchical vision transformer using shifted windows. In Proceedings of the IEEE/CVF international conference on computer vision, pages 10012 10022, 2021.
[10] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger. On calibration of modern neural networks. In International conference on machine learning, pages 1321 1330. PMLR, 2017.
[11] Quynh Nguyen and Matthias Hein. The loss surface of deep and wide neural networks. In International conference on machine learning, pages 2603-2612 PMLR, 2017.
[12] Marco Gori and Alberto Tesi. On the problem of local minima in backpropagation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(1):76-88) 7666,,,666666 1992.
[13] Zeyuan Allen Allen-Zhu, Yuanzhi Li, and Zhao Song. A convergence theory for deep learning via wwer-parameterizationn In International Conference on Machine Learning, pages 242-25 252.. PMLR, 2019.
[14] Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai. Gradient descent finds global minima of deep neural networks. In International conference on machine learning, pages 1675 1685. PMLR, 2019.
[15] Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu. Gradient descent optimizes over-parameterized deep relu networks. Machine learning, 109(3):467 492, 2020.
[16] Difan Zou and Quanquan Gu. An improved analysis of training over-parameterized deep neural networks. Advonces in neural information processing systems, 32, 2019.
[17] Faming Liang, Qizhai Li, and Lei Zhou. Bayesian neural networks for selection of drug sensitive genes. Journal of the American Statistical Association, 113(523):9555-99) 972, 2018.
[18] Nicholas G Polson and Veronika Ročková. Posterior concentration for sparse deep learning. Advances in Newral Information Processing Systemss 31, 2018.
[19] Yuexi Wang and Veronika Rocková. Uncertainty quantification for sparse deep learning. In International Conference on Artificial Intelligence and Statistics, pages 298 308. PMLR, 2020.
11

---

[20] Yan Sun, Qifan Song, and Faming Liang, Consistent sparse deep learning: Theory and
computation. Journal of the American Statistical Association, pages 1 1-15, 2021.
[21] Yan Sun, Wenjun Xiong, and Faming Liang. Sparse deep learning: A new framework immune to local traps and miscalibration. Advances in Neural Information Processing Systems, 34:22301-2253
2021.
[22] L. Breiman. Statistical modeling: The two cultures (with comments and a rejoinder by the author). Statistical Science, 16:199 2019,2211 2001.
[23] Johannes Schmidt-Hieber. Nonparametric regression using deep neural networks with relu activation function. The Annuls of Slatislics, 48(4):1916 1921, 2020.
[24] Helmut Boleskei, Philipp Grohs, Gitta Kutyniok, and Philipp Petersen. Optimal approx- imation with sparsely connected deep neural networks. SIAM Journal on Mathematics
of Data Science, 1(1):8 45, 2019.
[25] Vitaly Maiorov and Allan Pinkus. Lower bounds for approximation by mlp neural networks. Newrocomputing. 25(1-3):81-91, 1999.
[26] Jing Lei, Max G'Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman. Distribution-free predictive inference for regression. Journal of the American Statistical Association, (236553):094411111 2018.
[27] Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman. Inductive confidence machines for regression. In Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19-23, 2002 Proceedings 13, pages 345-356. Springer, 2002.
[28] Victor Chernozhukov, Kaspar Wüthrich, and Zhu Yinchu. Exact and robust conformal inference methods for predictive machine learning with dependent data. In Conference On learning theory, pages 732-749. PMLR, 2018.
[29] Chen Xu and Yao Xie. Conformal prediction interval for dynamic time-series. In International Conference on Machine Learning, pages 11559 11569. PMLR, 2021.
[30] Rina Foygel Barber, Emmanuel J Candes, Aaditya Ramdas, and Ryan J Tibshirani. Conformal prediction beyond exchangeability. The Annals of Statistics, (1(2):816 845, 2023.
[31] Margaux Zaffran, Olivier Féron, Yannig Goude, Julie Josse, and Aymeric Diculeveut. Adaptive conformal predictions for time scries. In International Conference on Machine Learning, pages 25834 25866. PMLR, 2022.
[32] Isaac Gibbs and Emmanuel Candes. Adaptive conformal inference under distribution shift. Advances in Newral Information Processing Systems, 34:1660-1672, 2021.
[33] Kamile Stankeviciute, Ahmed M Alaa, and Mihaela van der Schaar. Conformal time- series forecasting. Advances in Neural Information Processing Systems, 34:6216-6222 2021.
[34] Anastasios N Angelopoulos and Stephen Bates. A gentle introduction conformal pre- diction and distribution-free uncertainty quantification. ar preprint arXio:2107.07511 2021.
[35] Ruofeng Wen, Kari Torkkola, Balakrishnan Narayanaswamy, and Dhruv Madeka. A multi-horizon quantile recurrent forecaster. Xiv preprint arXiv: 1711. 11053, 2017.
[36] Yarin Gal and Zoubin Ghahramani. Dropout as bayesian approximation: Representing model uncertainty in deep learning. In international conference on machine learning, pages 1050 1059. PMLR, 2016.
[37] Derrick T Mirikitani and Nikolay Nikolaev. Recursive bayesian recurrent neural networks for time-scrics modeling. IEEE Transactions on Neural Networks, 21(2):262-274, 2009.
[38] Shiqing Ling and Michael McAleer, A general asymptotic theory for time-series models. Statistica Neerlandica, 64(1)):7-111 2010.
[39] Bin Yu. Rates of convergence for empirical processes of stationary mixing sequences. The Annals of Probebility, pages 94 116, 1994.
[40] Ron Meir. Performance bounds for nonlinear time series prediction. Tn Proceedings of the tenth annual conference on computational learning theory, pages 122 129, 1997.
12

---

[41] Cosma Shalizi and Aryeh Kontorovich. Predictive pac learning and process decomposi- tions. Advances in neural information processing systems, 26, 2013.
[42] Subhashis Ghosal and Aad van der Vaart. Convergence rales of posterior distributions for non-i.i.d. observations. Annals of Slatistics, 35:192 223, 2007.
[43] Ismael Castillo and Judith Rousseau. A bernstein von mises theorem for smooth functionals in semiparametric models. Annals of Slalislics, 43:2353 2383, 2013.
[44] Olivier Wintenberger. Optimal learning with bernstein online aggregation. Machine Learning, 106:119-1444 141, 2017.
[45] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Newral Comput., 9(8):1735-1788) nov 1997.
[46] Alistair EW Johnson, Tom J Pollard, Lu Shen, Li-wei H Lehman, Mengling Feng, Mohammad Ghassemi, Benjamin Moody, Peter Szolovits, Leo Authony Celi, and Roger G Mark. Mimic-iii, a freely accessible critical care database. Scientifie data, 3(1):1-9, 2016.
[47] Catherine Blake. Uci repository of machine learning databases. http://wwnee ics. uci. edu/" mlearn/Maa /LRRpoostoory. html, 1998.
[48] GOV UK et al. Coronavirus (covid-19) in the uk, 2020.
[49] Souhaib Ben Taieb and Amir F Atiya. A bias and variance analysis for multistep-aheut time series forecasting. IEEE tronsactions on neural networks and leorning systems, 27(1):62-76, 2015.
[50] Marcelo C Medeiros, Timo Teräsvirta, and Gianluigi Rech. Building neural network models for time series: a statistical approach. Journal of Forecasling, 25(1):49 75, 2006.
[51] Timo Terasvirla and Auders Kock. Forecasting with nonlinear time series models. Econometries: Mulliple Equation Models eJournal, 2010.
[52] Bjørn Auestad and Dag Tjøstheim. Identification of nonlinear time series: First order characterization and order determination. Biometrika, ((()))6696887, 1990.
[53] Wenxin Jiang. Bayesian variable selection for high dimensional generalized linear models: convergence rates of the fitted densities. The Annals of Statistics, 144) 2007.
[54] Subhashis Ghosal and Aad Van der Vaart. Fundamentals of nonparametric Bayesian inference, volume 44. Cambridge Universily Press, 2017.
[55] Andre M. Zubkov and Aleksandr A Serov. A complete proof of universal inequalities for the distribution function of the binomial law. Theory of Probabilily & Its Applications, 57()):539 544, 2013.
[56] Pentti Saikkonen. Dependent versions of a central limit theorem for the squared length of sample mean. Statistics & probability letters, 1253):185-194, 1995.
[57] Stephen Portnoy. Asymptotic behavior of likelihood methods for exponential families when the number of parameters tends infinity. The Annals of Statistics, pages 356-366,
1988.
[58] Yi-An Ma, Tianqi Chen, and Emily Fox. A complete recipe for stochastic gradient meme. Advances in newral information processing systems, 28, 2015.
(59] Grigoris A Dourbois and Pandelis N Biskas. European market coupling algorithm incorporating clearing conditions of block and complex orders. In 2015 IEEE Findhoven PowerTech, pages 1-6. IEEE, 2015.
[60] Zhen Lin, Shubhendu Trivedi, and Jimeng Sun. Conformal prediction with temporal quantile adjustments. arXiv preprint arXim:205 099400 2022.
[61] Davis Blalock, Jose Javier Gonzalez Ortiz, Jonathan Frankle, and John Guttag, What is the state of neural network pruning? Proceedings of machine learning and systems, 2:129-146 2020.
[62] Torsten Hoefler, Dan Alistarh, Tal Ben-Nun, Nikoli Dryden, and Alexandra Peste. Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks. J. Mach. Learn. Res., 124 2021.
13

---

[63] Wojciech Zaremba, Ilya Sutskever, and Oriol Vinyals. Recurrent neural network regu- larization. arXiv preprint orXiv:1409. 2329, 2014.
[64] Nadezhda Chirkova, Ekaterina Lobacheva, and Dmitry Vetrov. Bayesian compression for natural language processing. aiv preprinl arXiv:1810. 10927, 2018.
[65] Maxim Kodryan, Artem Grachev, Dmitry Ignatov, and Dmitry Vetrov. Efficient language modeling with automatic relevance determination in recurrent neural networks. In Proceedings of the 4th Workshop on Representation Learning for NLP (RepL4NLP.
2019), pages 40-48, 2019.
[66] Ekaterina Lobacheva, Nadezhda Chirkova, and Dmitry Vetrov. Bayesian sparsification of gated recurrent neural networks. arXiv preprint
2018
[67] Artem M Grachev, Dmitry I Tgnatov, and Andrey V Savchenko. Compression of recurrent neural networks for efficient language modeling. Applied Soft Computing, 79:354 362, 2019.
[68] Michael Zhu and Suyog Gupta. To prune, not to prune: exploring the efficacy of pruning for model compression. arxiv preprint Xiv: 110000188877 2017.
14

---

Appendix for "Sparse Deep Learning for Time Series Data: Theory and Applications"
A Mathematical Facts of Sparse RNNs
For sparse RNN model with H, recurrent layers and single output layer, several mathematical facts can be established. Let's denote the number of nodes in each layer as L,y..., LH.- Additionally, let rw, represent the number of non-zero connection weights for w', and let To denote the number of non-zero connection weights for v', where w, and v; denote the connection weights at two layers of sparse RNN.
Furthermore, we define (,, xx.1) as the output value of the j-th neuron of the i-th layer at time t. This output depends the parameter vector 8 and the input sequence Lt.1 = (att,.... ..t,...,221 . 11), where 24 represents the generic input to the sparse RNN at time step t.
Lemma A.1. Under Assumption 3.7 if sparse RNN has at most non-zero connection Hn
weights (i.e., EE Tw,
Bn Li=1 rr,) and || Ens then the summation of the absolute outputs of the ith layer at time is bounded by
(1, 1
II rwn
Twn
II Tuve sk=1
(Fr) Isis H, - 1,
k=1
L O'Hn. (3, x(:1)1 s t
Hn
Hn
Hn-1
t-1
II
II
Tiu
(En) (HHn
(En)
Proof. For simplicity, we rewrite (B, a11) as Oh when appropriate. The lemma is the result from the following facts:
For = and 1 sis Hr (Temma from
20
L
Lj
j=1
For 2 2 and i = (recursive relationship)
09 S Entwy + Ento,
1
For 12 2 and 2 2ii H, (recursive relationship)
L Σ
h--1
L Σ
10tyl & F Tws,
Σ
0 1
+ Enr..
0--1
1=1
1=1
For t2 2 and i = Hn (recursive relationship)
LR.
L
Σ
1O'H..s I S EnTwae
10
1H.-1sh-
We can verify this lemma by plugging the conclusion into all recursive relationships. We show the steps to verify the case when 2 2 and 2 < i < H, 1, since other cases are trivial
15

---

to verify.
La-1
I Σ 1 + EnIo
Oggl S Entws
1
1=1
<t-1
II Tux
(En)
(
Tws
k=1
+ (t 1)²
Tv, (tn)(tt
Tark
Tek
k=1
<1-1
Twn
II Tix
(En)*
+ 1)¹
II Twn
II ruv k=1
(En)
k=1
St
II Twn
(E,)
Fux
k==1
where the last inequality is due to the fact that
gi-1
(t - 1)a
< + 1-
+
=1.
t
t
Lemma A.2. Under Assumplion 3.7 consider two RNNs, B and B, where the former one a sparse network satisfying || Tn and || Ens and it's network structure vector
is denoted by 7. Assume that if [8 B,1 < 61 for all i E and B, -B,1 < 82 for all & 7 then
max Iu(EE:1, (3) u( B)1 #i1
Hr.
Hn-1
t-1
HD 81 H,
II Twe
(En + (,),(1+1))2
II
For
k=1
Hn bx)
Hr
HH
(-1
1HH022(pp
[1 + Twn (En + tdi]]
II
Ld2L + Fop. (on + 81)]
k=1
Proof. Define B such that B, B, for all y and for all i y. Let Oi denote O 0, (3, am.). Then, from the following facts that
For / and 1 S i H, (Lemma S4 from 20
L
o 0fg1 S (E, + 61)
rine
For 2 2 and i i= (recursive relationship)
Σ -1,11 $ rw, 81 + re, 81 Σ 1 (En + 81)
Σ
1 0 |.
j=1
For 2 2 and 2 sis H
(recursive relationship)
L Σ
Li-1 O u.l + Tw, Tw, + 61)
Li-1
1 - 0.51 rw, 81
Σ
1 - O.
j=l
k=
k=1
L
L.
Tey 61
1 + Tw, (En + 81)
1 -
k=1
16

---

For t 2 and i = Hn (recursive relationship)
LH Σ 2
L
Ö O'n..
O'Hni S Twnn 61
Σ
10
Hn -1,8
1=1
Lun-1
+ Twua
(En + ò()
Ö Hr 11, - O2 -1,6-
We have
For isis Hn
Σ
II riee
Il T'v
(En + 81) (ittle))-
For Hn
|p( ()) -1p(x111, /))
2
in..s O'n...g
Hn
Hn -1
f-1
S tn Hn
II
(En - (ta
Tark
Tek
We can verify the above conclusion by plugging into all recursive relationships. We show the steps to verify the case when t 2 2 and 2 < i < Hn since other cases are trivial to
17

---

verify.
1--1
L--1
Ö O,,11 S rw. 81
10 11 + +rw, (En + +di)
00 1 1-1 11-1,
==1
k=1
1
+ Te, δ1
1
1 + +Te (En + +61)
2, - -
k=1
k=1
6-1
S 6 1-
II Tunk
II k=1
(En + &1)
Tix
i-1
1-1
+ 18 - 1)
Tink
II
Tix
(En + 0
x=1
i II
t-2
+ 11 - 1)a
Tank
Fon
Tn; (En + :81)
k=1
k=1
1-2
+ +1 - 1)tt1
(Fin + +51) 51)
R=1
k=1
t-1
S 6,1-11
II
(Fn + 51) (1+1)0-2
Tive
A==1
(-1
+ (1 - 1)
II
II A=11
(En + 56 (14
Tine
W=1
i
1-1
+.1 - 1)x
II
(En + +5z)
Fwn
Fen
A=1
1-1
+ 5t - 1)it1
IIrua Tick
II km
(En + (6411))8
Tix
t-1
$ 18
II Twn
II Tee
(En + 114
where the last inequality is due to the fact that for and 2 2, it is easy see that
1-1 + (i
1)0' + (1 - 1)a + i(l 1)tt1
i-11
i+1
+
ittt1
it
it t
t
1
(8-10 1)2
1
3
t+
2t
+12
- 21 +
2(t
1)3
+
212
2t3
2t
31
+
< 1.
2+3
Now let Öi denotes O1,5, (B, mnx), then from the facts that
For 1 < T and sis Hn
L
t-1
Σ 1
5t
Twn
II
Tex
(En (En+
isis 11, 1stsT,
k=
-=
Lua
Hn
Hn-
8-1
Σ
1 ÖHnul S H4
Tws
T's
(En, + 81) 11)+++
1stsT.
k=1
A=1
18

---

This can be easily derived by Lemma A.1 and the fact that B
En + 81.
For = 1 and 1 i Hn (Lemma S4 from 20
10 - - s 82 Ia Σ
Σ
11 + rw (Pn + +5i)l
1=1
k=1
k=
For 12 2 and i (recursive relationship)
La
In Σ
In + + (En + +61)
Σ
10 - Oi S L1L1nn + 82h1
1
-Oi,h
Σ
- 01,11
1=1
k=1
k=1
L
+ ااعة
1
k=1
For t2 2 and 2 2 H
(rocursive relationship)
L
11-1
Σ 0h 01,1 s 82
Σ
1 0111, rw;
(Fr Fun
+ +6)))
Σ
Ó
-O1-1K
k=1
L1-1
8g I4
Σ
10;
821.1 Σ
1 -
k=1
Li
+ Tv. E + 81)
1
- 1| + 8ZLر
Σ
i
k=
k=
For 2 2 and i Hn (recursive relationship)
Lun
LHa
i'n... - Öirm.sl 52LH
1 -1,-1,1
0 11,
Lin
Twna
(En +
Ö H..k O'Hm-1.h
Lun-1
+ 82LH
H, ÖHH,-1,1
k=1
We have
For sis H
4.
Σ
1 وو,ن s & Pn Pnl + Σ Tok
[
182T + rw (Fa (Fn + 51)]
1=
k=1
II
182Lx Tun (En + 813]
For = Hn
Lu
|p( (B)
ul --1=
Σ
Ö O
O O'Hnol
F=1
Hn.
Hn
HH -1
6-1
< Hn Prtto +
Tak
II
T2 + rw, (Fn + 51)]
To + row (Fn + 51)]
We can verify the above conclusion in similar approach. The proof is completed by summation of the bound for |u( Iu(0,1, 3)1 and [u ()) u(
19

---

B Proofs on Posterior Consistency: A Single Time Series
To establish posterior consistency for DNNs with i.i.d. data, 20 utilized Proposition 1 from [3.. This lemma provides three sufficient conditions for proving posterior consistency for general statistical models with i.id. data, along with a posterior contraction rate. In this paper, we aim to establish posterior consistency for DNNs with stochastic processes, specifically time series, that are strictly stationary and a-mixing with an exponentially decaying mixing coefficient 54 42.
Consider time series Dn (21, 22....,,2) 2n) defined probability space F,, P*), which satisfies the assumptions outlined in Assumption 1. For simplicity, we assume that the initial values 21-1,--.;20 Zo are fixed and given.
Let Px, denote a set of probability densities, let IP denote the complement of P,, and let 62 denote sequence of positive numbers. Let NN (N,,,, PPn) be the minimum number of Hellinger balls of radius en that are needed to cover Pn, i.e., (en, PP.) is the minimum of all number k's such that there exist sets Sj = {p: d(p.p)) S end.j I...,k, with P, C S holding, where
d(p.9)
denotes the integrated Hellinger distance p(2[z) and g(2|=), where a E TR' contains the history up to time steps of 2, and v(x) is the probability density function of the marginal distribution of æ. Note that is invariant with respect to time index i due to the strictly stationary assumption.
42 54
between the two conditional densities
For Da, denote the corresponding true conditional density by p². Define T(-) as the prior density, and (T as the posterior. Define #(c) = π rld(p.p") > el Dn] for each e > 0. Assume the conditions:
(a) log N(e,, Pn) S ne, for all sufficiently large n. -bbbn for some b > 0 and all sufficiently large n.
(b) xPP)) Se bnch
(c) Let (f
then <(Ipp - -plll S blcn) 2 e e- for some 84 > 0, 2, and all sufficiently large
Lemma B.1. Under the conditions (a), (b) and (c), given sufficiently large lim, 0, and lim n 00, we have for some large C > 0,
(CC,)) = (d(p.p") 2 Cc, Dr) = Op (c (n..) =(Cc.)=
Proof. This Lemma can be proved with similar arguments used in Section 9.5.3, Theorem 8.19, and Theorem 8.29 of 54
B.1 Postcrior Consistency with a General Shrinkage Prior
Let B denote the vector of all connection weights of RNN. To prove Theorem 3.9 first consider a general shrinkage prior that all entries of B are subject to an independent contimuous prior Tb. iee, ((3) = I Ta6),, where K₂ denotes the total number elements of B. Thcorem B.2 provides a sufficient condition for posterior consistency.
Theorem B.2. (Posterior consistency) Suppose Assumptions 3.1 3.7 hold. If the prior ma satisfies that
log(1/I,) = O(H, log + log L),
(8)
1 expl- log + log L + p.]).
Toll-Th. Thi mm] 2
(9)
Kn
1
Tb(l-n/a n/ll 2
(10)
Kn,"
- K, (1 (1811> M,)] > nc?
(11)
20

---

for some > 0, where
The < 1/|vTMM (n/HH, 2 (n) 2MM,H.HHHHH..
I/n < 1/lVnM/ (r./ H H,,) )2M (co Em. H.
with some cu > Th is the minimal density value of To within the interval E-E, - E, +1, and Mn is some sequence salisfying log( Mn) = O(lg(())) Then, there erists a sequence en satisfying ne = rn Hn log n + In log + 8n log Po + nwn and en < <1, such that
# = = =ddpp,p") 2 Dn) = Op- (e (e
for some large M > 0.
Proof. To prove this thcorem, suffices to check all three conditions listed in Lemma B.1 Checking condition (c):
Consider the set A = {B maxjer 18, 331 S wn, jax 13g 181 S w/.. where
CiEn
Wn MH Hn Hn Hny ) Ha (co En )2M.H. )n)Mi
w/ S M1 Kn (r/ Hn Hm En) En)di,H H2M
CiCn
we have |u( 9) p( 9r) 3c By the definition 4 we have
for some ei and co 1. If B E A, by Lemma A.2
[u(aM;,, n² (aMaa)) S 3 wn Finally for some s > 2 (for simplicity, we take s to be an even integer),
llg - Pa- ll. (/ (Iu(zMaa
1/8
(f (3c1en + wn) (o((MM))dd
1/s
S 30 + wn.
For any small b' > 0, condition (c) is satisfied as long as is sufficiently small, c. 2 Cow, for some large Co, and the prior satisfies - log (A) S yne Since
x (A) 2 (mo((En - Way Er, + wnl)) Tm x ((max 1ر S w/)))
387
z (2mjwn)"
w/r) ( Kn-
2 (2Ijwn) "m (1
(1-1/Km))
where the last inequality is due to the fact that w/, >> The Note that lim, 1/K,, e1, Since log(1/wn)) x H, log(En) + 2M, log(r,/ H.)) + log(1/en) + constant. O(H, log (note that log(1/en)
n)). then for sufficiently large n,
(1]
2)) + 2
log (A) S rn log
+ r. O(H, log(n)) + r. log
= Or,, Hr log(n) + log L).
Thus, the prior satisfies - (A) S yn for sufficiently large n, when ne > Co(r,H logn n + r, log L) for some sufficiently large constant Co. Thus condition (c) holds. Checking condition (a):
Let P,, denote the set of probability densities for the RNNs whose parameter vectors satisfy
B B, = {188] S Mn = (i: 1/:1 2 8',} satisfies [Tg] S kignn Lelin S Kn
21

---

where J denotes the number of input connections with the absolute weights greater than 8n ( n/rn) and K, n/s,, will be specified later, and
8, Mf Kn (kar (knn)/HH))))))) HH,,)H2,),, 2,, (cgMM)) )MM,
C1cn
for some 61 > 0, O. Let
c1en
on Mf Hn (karn/ Hn, ) H (co (co.Mn)2M.... )2M,
Consider two parameter vectors 9 and 8 in set Ba, such that there exists a structure y with |y| kn and hli S kn; and [ 88] S 6n for all j E 7, max(18j). (8)) s 6, for all $ y. Hence, by Lemma A.2 we have that |u(x1;M, 8*)) - p(z1. }²2 - 9cfe, e. For
two normal distributions N (1i, a and N(22, (2)2,0²), define the corresponding Kullback-Leibler divergence as
do(p Pwz) = 20 P2 m1l3.
Together with the fact that 2d Pn2) 1 du(pp) pp2)) have
dpp. pa & 4o(pa" .p=) - + (()
for some C > 0, given sufficiently small 61-
2M,
Given the above results, one can bound the packing number (I en) by Ejnr XH. 8n where X², denotes the number of all valid networks who has exact j connections and has no more than K, inputs. Since log XHr k' Kn log + log(k'sn Li + 2H, 22), log Mn Ollog n), ka I and Kn,nn I then
knrn 2Ma Ka
log N (n) S log hi knraxin
8n
log knrn + k 8n log Pr. kn log(2H,,) + 2k,r, log [L+ K 8m)
+ knr. log M, H (k vn/ / (2)) (n
CIEn
=k,,, (HH, log + log L) + sn log Pn. Pn.
We can choose kn and K, such that for sufficiently large n, knrn Hn log n + log L} x K, log x ne? and then log ((Pas S ne?
Checking condition (b):
Lemma B.3. (Theorem 1 of M5 Let X Binomial p) be Binomial random variable. For any < k < 1
1 - mmp)[2nH((ppp k/n))///))))
where & is the cumulative distribution function (CDF) of the standard Gaussian distribution and H(p,k/n) (kg(k)))())/p)) k/n) log[(1 k/n)///n p)].
Now,
T <(PM) < Pr (Binomial un) > knTw) + KnTo(13]] > Mn) + Pr( > K,,n)),
where Un = 1 - ( (1-6, 8,,)) By the condition on Ta and the fact that 6', > 1// it is easy to see that U, < (H, log n + log L. + log log K,) for some constant 7 > 0. Thus, by Lemma B.3 log Pr Pr(BBoomiaa(KK (K, v.,) > k,rw,) ) xs Tk,r.. HH,, log n + log + logPa] 2 nc?
due to the choice of kn, and - log Pr([[ g sn 2 (,) = W,, log + log L + pog log (pn Kn [T]] Z ne? due to the choice of k,, Thus, condition (b) holds as well. Then by
lemma B.1 the proof is completed.
22

---

B.2 Proof of Theorem 3.9
Proof. To prove Theorem 3.9 in the main text, it suffices to verify the four conditions on To listed in Theorem B. Let Mn = max 12 En). Condition can be verified by choosing 0n. such that E2/20f.n + log oogor.n = OHH, log + log L.. Con- ditions 9 and 10 can be verified by setting An = 1/(M," Kn 2M.H. (Ip.)]"> and Jo, < [M, vnn ( (cgM, ),M, Finally, condition 11 can be verified by Mn 2 2n and TH,, log n + log L + log Pm] + 22//20 2 n. Finally, based on the proof above we see that ne,, x Tr log n + rn log L sn log Pn now and lim En 0.
B.3 Posterior Consistency for Multilayer Perceptrons
As highlighted in Section 3 the MLP can be formulated as a regression problem. Here, the input is æ; 1/-1-i for some 1 Rí < n, with the corresponding output being y:- Thus, the dataset Dn can be represented as (x), yi Vi) We apply the same assumptions as for the MLP, specifically 3.1
3.4
and 3.7 We also use the same definitions and notations for
the MLP as those in 20
21
Leveraging the mathematical properties of sparse MLPs presented in 20 and the proof of sparse RNNs for a single time series discussed above, one can straightforwardly derive the following Corollary. Let d(p1.P2) denote the integrated Hellinger distance between two conditional densities p1 xy))) and p2(y(z). Let x(-(Dn) be the posterior probability of an event.
3.4 and 3. hold. If the misture Gaus- Corollary B.4. Suppose Assumptions 3.1
sian prior satisfies the conditions
K Knnn (Ipm))")) for some con- stant 0, Hn/(H, log + log L2 I 01, n² for some constant a, and Ou,- S min(1/lynK,, , /Hn)
An
O
HH, then there exists an an error se- quence en = ( + O(2)) such that lim En 0 and lim ne? 0c, and the posterior distribution satisfies
x(d(ps.Pa)) Cenl Da) = Pp (e-me))
(12)
for sufficiently large n and C > 0, where 22 = Hr log log &+ logpn//n, Pp denotes the underlying true data distribution, and Pa denotes the data distribution reconstructed by the Bayesian MLP based on its posterior samples.
C Asymptotic Normality of Connection Weights and Predictions
This section provides detailed assumptions and proofs for Theorem 3.12 and 3.13 For simplicity, we assume y0:1-n, is also given, and let Mi R₁ 1.
Let (.
log(ps
denote the likelihood function, and let (B) denote the den- (B)
ag,, .33 which denotes the d-th order partial derivatives. Let Hn(3) denote the Hessian matrix of (,),), and let ( d (B) denote the (i.j)-th component of H, (B) and H (B), respectively. Let BA A2/2 (B8)) (Bt) and bx.n = Vwn/nBB, For RNN with B, we define the weight truncation at the true model structure >" (B,,); = B, for i E 7* and (8y.): = for i i y For the mixture Gaussian prior (3), let Bs, (B*) = {B: : 18, - 3 < 6n. Vi E 1,, 8
sity of the mixture Gaussian prior (3). Let hi Aa
1 Vi }*}
n0 log Anooon
In addition, we make the following assumptions for Theorem 3.12 and Theorem 3.13 Assumption C.I. Assume the conditions of Lemma
1
3.10 hold with p(en) ol Kn the C 2 defined Assumption 3.4. For some on s.t. I' S 8n S Varn" In
and
2
let A(Gn,5m)
23

---

{B max 1 3 ön, d..d(p..Pu' S ) where is the posterior contraction rate as Assume there exists some constants C' and M > 0 such that
defined in Theorem 3.9
C.1 8" 13 (x,)) is generic, minkey |3 Côn and π (A((n, 8,)) n 20.
as
C.2 h(()))) < M. ( < M, 11 (B*) < M, (b) M, (1: < M hold for any i.j. k E 7" . 7", and B E B2s,. (8").
C.3 sup(LE((d'U3)] 11 8-1 - 1.2 = | BA = O1), where U Z Egg.- (Z), Z denotes a random variable drawn from neural network model parameterized by 8T,, and Eg,- denotes the mean of Z.
and
C.4 r 0 and the conditions for Theorem 2 of 56 hold.
Conditions C.1 to C.3 align with the assumptions made in 21. An additional assumption, C.4, is introduced to account for the dependency inherent in time series data. This assumption is crucial and employed in conjunction with C.3 to establish the consistency of the maximum likelihood estimator (MLE) By- for a given structure y². While the assumption rz/n + is sufficient for independent data, which is implied by Assumption 3.4 stronger restriction, specifically r//n n + 0, is necessary for dependent data such as time scrics. It worth noting that the conditions used in Theorem of 56 pertain specifically to the time series data itself.
8),
Let Pa (1 B) 03;,- ..,,,, denotes the d-th order partial derivative for some input.
Assumption C.2. |us(, < M, (+ < M, [mg(- 8) < M for any i,j E "rj=r''0 y",
and B E B28, (B*), where M as defined in Assumption C.1
Proof of Theorem 3. 1.12 and Theorem 3.13 The proof of these two theorems can be
conducted following the approach in [2I] The main difference is that in [21], they rely on Theorem 2.1 of [57], which assumes independent data. However, the same conclusions can be extended to time series data by using Theorem 2 of [56), taking into account assumption C.4. This allows for the consideration of dependent data in the analysis.
D Computation
Algorithm D gives the prior annealing procedure 21 In practice, the following implementation can be followed based on Algorithm I
For 0 Ti, perform initial training.
For Ti S $ T2, set o = o,n
init and gradually increase r7( = /
For Tı t T3, set (1)) 1 and gradually decrease of according to the formula:
(t) = T5-1: Ense + T- (o
d.n
For t > T3, set T7()) = 1, 0,n = end and gradually decrease the temperature 7 according to the formula: f . where constant.
(t)
Please refer to Appendix E for intuitive explanations of the prior annealing algorithm and further details on training a model to achieve the desired sparsity.
E Mixture Gaussian Prior
The mixture Gaussian prior imposes penalty the model parameters by acting as piecewise 1.2 penalty, applying varying degrees of penalty in different regions of the parameter space. Given values for no Lin, and Xn threshold value can be computed using Algorithm
24

---

Algorithm 1 Prior Annealing: Frequentist
[1] (Initial Training) Train a neural network satisfying condition (S*) such that a global optimal solution Bo = arg max 3 (3) reached, this stage can be accomplished by using SGD or Adam.
[2] (Prior Annealing) Initialize 8 at Bur and simulate from sequence of distributions
T ( T, (7)) (oo)) , x edt, n/ T (T3 for & 1.2..... . m, where 0 17()) < (2)
S mm 1, Th X,, (1-))) ()))))), and o 0.n
(1) (2)
0.n
(m) ond. This can be done by using stochastic gradient MCMC algorithms [S
After the stage m has been reached, continue to run the simulated annealing algorithm by gradually decreasing the temperature to a very small value. Denote the resulting model
by B.
[3] (Structure Sparsification) For each connection weight i E (1,2,.... K,), set 7 = 1,
V200,
An
T1 and 0 otherwise, where the threshold value is Ju,
18
log
Vot. 28.0
An
determined by solving ((>; = 119) 0.5, and 00, and Denote the the structure of the sparse RNN by 7.
[4] (Nonzero-weights Refining) Refine the nonzero weights of the sparse model (3, 5) by maximizing L, (3). Denote the resulting estimate by (3,
A, = 1e 7. 10 =.005
Init "free space"
Init "free space*
Negative probability
Mid "free space"
Mid "free space*
End "free space'
End "free space"
-3
-0.05
0.00
0.05
0.10
0.15
0.20
B
Figure Changes of the negative log prior density function in prior annealing of oi where the values of An and of,n are fixed.
I Parameters whose absolute values are below this threshold will receive large penalty, hence constituting the "penalty region', while parameters whose absolute values are above the threshold will receive relatively minimal penalty, forming the "free space'. The severity of the penalty the free space largely depends on the value of 01,a-
For instance, as depicted in Figure 3 based on the simple practice implementation detailed in Appendix D we fix A, 1e 7 7, of 0.05, and we set (oguit)² le 5, (a 1e and gradually reduce of from the initial value the end value. Initially, the free space is relatively small and the penalty region is relatively large. However, the penalty imposed on the parameters within the initial penalty region is also minor, making challenging to shrink these parameters to zero. As we progressively decrease o the free space enlarges,
25

---

the penalty region diminishes, and the penalty on the parameters within the penalty region intensifies simultancously. Once Th equals (o the parameters within the penalty region will be proximate to zero, and the parameters outside the penalty region can vary freely almost all areas the parameter space.
For model compression tasks, achieving the desired sparsity involves several steps post the initial training phase outlined in Algorithm First, determine the initial threshold value based on pre-set values for An, o... and (o Then, compute the proportion of the parameters in the initial model whose absolute values are lower than this threshold. This value serves as an estimate for anticipated sparsity. Adjust the (oope") value until the predicted sparsity aligns closely with the desired sparsity level.
F
Construction of Multi-Horizon Joint Prediction Intervals
To construct valid joint prediction intervals for multi-horizon forecasting, we estimate the individual variances for each time step the prediction horizon using procedures similar to one-step-ahead forecasting. We then adjust the critical value by dividing by the number of time steps m in the prediction horizon and use 2 (2/(2) as the adjusted critical value. This Bonferroni correction ensures the desired coverage probability across the entire prediction horizon.
As an example, we refer to the experiments conducted in 5.1.2 where we work with set of training sequences denoted by D, = (1) ..... VT -T 1 Here, (6) Vr+1 42....... VT represent the prediction horizon for each sequence, while T repre-
(4)
sents the length of the observed sequence.
Train a model by the proposed algorithm, and denote the trained model by (3, ()
Calculate o² estimator of o² 1=
1
(8) [u(3,v19) vitir+m) & (u(8.v13) v141IT+m)
where n1 V1.r)
11 1 and @ denotes elementwise product.
For test sequence V1.T calculate
(0)
Σ = V 1 (
(0)
Let & mx denote the vector formed by the diagonal elements of Σ.
The Bonferroni simultancous prediction intervals for all elements of 9/ 41:744 given by
(0)
1BB (0) 1 2a/(2m) (-s 0²2 of
where of represents the element-wise square root operation.
The Bayesian method is particularly advantageous when dealing with a large number of non-zero connection weights in (3,4). making the computation of the Hessian matrix of the log-likelihood function costly unfeasible. For detailed information utilizing the Bayesian method for constructing prediction intorvals, please refer to 21
G Numerical Experiments
G.1 French Electricity Spot Prices
Dataset. The given dataset contains the spot prices of electricity in France that were established over period of four years, from 2016 to 2020, using auction market. In
26

---

this market, producers and suppliers submit their orders for the next day's 24-hour period, specifying the electricity volume in MWh they intend to sell or purchase, along with the corresponding price in C/MWh. At midnight, the Euphemia algorithm, as described in (00 calculates the 24 hourly prices for the next day, based on the submitted orders and other constraints. This hourly dataset consists of 35064 observations, covering (3 x 365 + +366) 24 periods. Our main objective to predict the 24 prices for the next day by considering t--oftthh-week,
different explanatory variables such as the day-ahead forecast consumption, as well as the prices of the previous day and the same day in the previous week, as these variables are crucial in determining the spot prices of electricity. Refer to Figure for visual representation of the dataset.
The prediction models and training settings described below are the same for all 24 hours.
Bectricty Spot Prices
Unaadd
tainng
Celberor
3226 2017
2029
200
Figure 4: (a): French electricity spot prices (2016 to 2019). (b): Online learning method usod for ACI and AgACI in 31. where the prodiction range is assumed to between Tredd start and Tpred end.
Prediction Model. We use MLP with one hidden layer of size 100 and the sigmoid activation function as the underlying prediction model for all methods and all hours.
Training: Baselines. For all CP baselines, we train the MLP for 300 epochs using SGD with constant learning rate of 0.001 and momentum of 0.9, the batch size is set to be 100. For the ACI, we use the same list of y E (0. 0.000005, 0.00005, 0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0001, 1.0001,0.0002, 0.0002, 0.0003, 0.0004, 0.0005, O. 0.0006, 0.0007, 0.0008, 0.0009, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09} as 31
Training: Our Method. For our method PA, we train total of 300 epochs with the same learning rate, momentum, and batch size. We use Ti 150, T2 160, Ta 260. We run SGD with momentum for I < T1 and SGHMC with temperature= 1 for >= Tí. For the mixture Gaussian prior, we fix o. 0.01, (afnir)2 le 5, (af, le and
An = 77
Table 5: Uncertainty quantification results produced by different methods for the French electricity spot prices dataset. This table corresponds to the Figure
Methods
Coverage
Average PI Lengths (standard deviation)
Median PI Lengths (interquartile range)
PA offiine (ours)
90.06
20.19(3.23)
20.59(1 (.8)
AgACI online
91.05
23.27(11.19)
21.38(2.62)
ACI 7 = 0.01 online
90.21
22. 22.13(8.500
20.68(2.70)
ACI 7 = 0.0 online
92.89
24.39(9.17)
22.99(2.91
EnbPi V2 online
91.23
27. (7.99( 110.17)
26.55
NexCP online
91.18
24.41(10.40)
22.88(2.86)
G.2 EEG, MIMIC-III, and COVID-19
EEG The EEG dalaset, available at here, served as the primary source for the EEG signal time series. This dataset contains responses from both control and alcoholic subjects who were presented with visual stimuli of three different types. To maintain consistency with
27

---

previous work 88, we utilized the medium-sizee version, consisting of 10 control and 10 alcoholic subjects. We focused solely on the control subjects for our experiments, as the dataset summaries indicated their EFG responses were more difficult to predict. For detailed information on our processing steps, please refer to 33
COVID-19 We followed the processing steps of previous research [* and utilized COVID- 19 data from various regions within the same country to minimize potential distribution shifts while adhering to the exchangeability assumption. The lower tier local authority split provided us with 380 sequences, which we randomly allocated between the training. calibration, and test sets over multiple trials. The dataset can be found at here
MIMIC-III We collect patient data on the use of antibiotics (specifically Levofloxacin from the MIMIC-III dataset, [46]. However, the subset of the MIMIC-II dataset used in [3] not publicly available us, and the authors did not provide information on the processing
steps, such as SQL queries and Python code. Therefore, we follow the published processing steps provided in P+ We remove sequences with fewer than 5 visits more than 47 visits, resulting in a total of 2992 sequences. We randomly split the dataset into train, calibration, and test sets, with corresponding proportions of 43%, 47%, and 10%. We the white blood cell count (high) as the feature for the univariate time series from these sequences. Access to the MIMIC-III dataset requires PhysioNet credentialing
Table presents the training hyperparameters and prediction models used for the baselines in the three datasets, while Table presents the training details used for method.
We hypothesize that the mediocre performance on the FEG dataset due to the small size of the prediction model, which has only 2025 parameters and is substantially smaller than the number training sequences. Hence, we conducted experiments on the EEG dataset using slightly larger prediction model for different methods. The results are presented in Table [S] and the correspouding training details are provided in Table G Again, our method outperforms other baselines as well.
Table 6: Training details for the baselines, where LSTM-d refers to single hidden layer LSTM network with a size of d.
MIMIC-III
EEG
COVTD-19
model learning rate
LSTM-500
LSTM-20
LSTM-20
2x
0.01
0.01
Epochs
100
1000
Oplimizer
Adam
Adam
Adam
Batch size
150
150
150
Table 7: Training details for our method, where LSTM-d refers a single hidden layer LSTM network with a size of d.
MIMIC-I1I
EEG
COVID-19
Model Learning rate
LSTM-500
LSTM-20
LSTM-20
2c
0.01
0.01
Epochs
65
1000
Optimizer
Adam
Adam
Adam
Batch size
150
150
150
of
0.05
0.1
0.1
(oapil snsty y
le
le - 7
le 7
(oudy An
1e
1e 8
1e 8
-
1e 7
1e 7
Temperature
1
Ti
40
50
800
T2
40
55
850
Ta
55
950
28

---

Table Uncertainty quantification results produced by various methods for the EEG data using larger prediction model. The "coverage" represents the average joint coverage rate over different random seeds. The prediction interval (PI) lengths are averaged across prediction horizons and random seeds. The values in parentheses indicate the standard deviations of the respective means.
EEG
Model
Coverage
lengths
PA-RNN
91.0%(0.8%)
40.84(5.69)
CF-RNN
90.8%(1.8% 15.2%(2.5%)
43.4(6.79)
MQ-RNN
20.63(2.07)
DP-RNN
0.6%(0.1%)
5.02(0.53)
Table Training details for the EEG data with a slightly larger prediction model for all methods, where denotes that the information is not applicable.
PA-RNN
CF-RNN
MQ-RNN
DP-RNN
model
LSTM-100
LSTM-100
LSTM-100
LSTM-100
learning rate
le
1e-3 3
1e-3
le 1e-3
Epochs
150
150
150
Optimizer
Adam
Adam
Adam 150
Adam
Batch size of imit (ormm
150
150
150
0.01
le - 5
(0
An
temperature
1
Ti
100
12
105
T3
130
G.3 Autoregressive Order Selection
Model An Elman RNN with one hidden layer of size 1000. Different window sizes (i.e., or 15) will result in different total number of parameters..
Hyperparameters All training hyperparameters are given in Table 10
Table 10: Training details for the autoregressive order selection experiment.
PA-RNN
PA-RNN 15
RNN
RNN 15
Learning rate Trerations Optimizer Batch size
4e
1e 4
1e
25000
25000
25000
25000
SGHMC
SGHMC
SGHMC
SGHMC
36
36
36
36
Subsample size per iteration
50
50
50
50
Predicton horizon
1
of (ouut end (o X
0.05
0.05
4e 6
1e 7
le 7
le 7
le 7
Temperature
0.1
0.1
Ti
5000
5000
Th
10000
10000
T3
25000
25000
29

---

G.4 Large-Scale Model Compression
As pointed out by recent summary/suun /survey рapers on sparse deep learning 61 72 the lack of standardized benchmarks and metrics that provide guidance on model structure, task, dataset, and sparsity levels has caused difficulties in conducting fair and meaningful comparisons with previous works. For example, the task of compressing Penn Tree Bank
(PTB) word language model [8] a popular comparison task. However, many previous works 64 65
66 [7 have avoided comparison with the state-of-the-art method by either not using the standard baseline model, not reporting, or conducting comparisons at different sparsity levels. Therefore, we performed an extensive scarch of papers that reported performance on this task, and to the best of our knowledge, the state-of-the-ar method is the Automated Gradual Pruning (AGP) by 6S.
In our experiments, we train large stacked LSTM language models on the PTB dataset at different sparsity levels. The model architecture follows the same design as in 63 comprising an embedding layer, two stacked LSTM layers, and softmax layer. The vocabulary size for the model is 10000, the embedding layer size is 1500, and the hidden layer size is 1500 resulting a total of 66 million parameters. We compare our method with AGP at different sparsity levels, including 80%, 85%, 90%, 95%, and 97.5% as in 68 The results are summarized in Table [T]] and our method achieves better results consistently. For the AGP, numbers are taken directly from the original paper, and since they only provided one experimental result for each sparsity level, no standard deviation is reported. For our method, we run three independent trials and provide both the mean and standard deviation for each sparsity level. During the initial training stage, we follow the same training procedure as in [ The details of the prior annealing and fine-tuning stages of our method for different sparsity levels are provided below.
During the initial training stage, we follow the same training procedure as in 63 for all sparsity, and hence Ti 55.
During the prior anncaling stage, we train a total of 60 epochs using SGHMC. For all levels of sparsity considered,
0.9, minibatch size = 20, and we fix T2 = 5 + Ti = 60, T3 = T2 + 20 80. We set the initial temperature == 0.01 for t < T3 and and gradually decrease by t- T3
fix of = 0.5, An le 6, momentum
0.01 for t >> Tz.
During the fine tune stage, we apply similar training procedure as the init training stage, Le., we use SGD with gradient clipping and we decrease the learning rate by a constant factor after certain epochs. The minibatch size is set to be 20, and we apply early stopping based on the validation perplexity with a maximum of 30 epochs.
Table 12 gives all hyperparameters (not specified above) for different sparsity levels.
Note that, the mixture Gaussian prior by nature is regularization method, we lower the dropout ratio during the prior annealing and fine tune stage for models with relatively high sparsity.
30

---

Table 11: Test set perplexity for different methods and sparsity levels.
Method
Sparsity
Test Perplexity
baseline
0%
78.40
PA
80.87% + 0.09%
77.122 + 0.147
PA
85.83% = 0.05%
77.431 + 0.109
PA
90.53% + 0.05%
78.823 + 0.118
PA
95.40% + 0.00%
84.525 + 0.112
PA
97.78% = 0.05%
93.268 = 0.136
AGP
80%
77.52
AGP
85%
78.31
AGP
90%
80.24
AGP
95%
87.83
AGP
97.5%
103.20
Table 12: Word Language Model Compression: Hyperparameters for our method during the prior anncaling and fine tune stage. We denote the learning rate by LR, the prior annealing stage by PA, and the fine tune stage by FT.
Hyyperaarameteer/Spaaiity
80%
85%
90%
95%
97.5%
Dropout ratio
0.65
0.65
0.65
0.4
0.4
(otnt) )²
0.0005
0.0007
0.00097
0.00174
0.00248
(06,))
7
7
le
LR PA
0.004
0.004
0.008
0.008
0.01
LR FT
0.008
0.008
0.008
0.1
0.2
LRR. decay factor FT
1/0.95
1/0.95
1/0.95
1/0.95
1/0.95
LR decay epoch FT
5
5
5
5
5
Table 13: Training details for the additional autoregressive order selection experiment.
V (0.8 1.1 expe 50y -00g }) Wi-1 + 1
PA-RNN 1,3,5, 7, 10, 15
Learning rate Iterations
le
2000
Optimizer Batch size Subsample size per iteration
SCHMC
36 50
Predicton horizon
of (nitt end (o
0.05
le
1e 6
Xm
1e 6
Temperature
0.1
500
T2
1000
T
1500
31


Reviews:
[{"review_id": "hdoBvuyq1E", "reviewer": "Reviewer_AFD6", "summary": "The paper extends the sparse DNN theory from [1, 2] to time series data and proposes the sparsity method for RNNs. This method aims to improve uncertainty quantification on time-series tasks and provide state-of-the-art RNN compression.", "questions": "1. Could the authors address the novelty and model compression evaluation concerns? \n2. In my understanding, deep ensembles are s.o.t.a. for uncertainty estimation. Is it a reasonable baseline in the multi-horizon setting?", "limitations": "The authors addressed the limitations to some extent. One can also discuss more limitations of the proposed algorithm.", "rating": 5, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "* The combination of DNN sparsification and uncertainty quantification for time series data is an interesting direction.\n* The theoretical groundwork seems complete and convincing.", "weaknesses": "**Novelty**\n\n* **Theory** - The theoretical results for RNNs and time series seem almost identical to [1, 2] aside from minor time series specifics. I appreciate that the authors explicitly discuss the connection to these works. However, I have concerns if the theoretical contribution is significant enough. \n* **Computation** - the prior annealing algorithm (PA) for model sparsification is identical to [2]. Constructing the prediction intervals for time series might be novel, but it is a minor contribution from my perspective. \n\n**Large-scale model compression**\n\nThe authors state that one of their contributions is s.o.t.a. large-scale RNN compression method. \n\n1. In my opinion, the scope of the model compression area is much larger than the weight sparsification. One should either reformulate \"compression\" to \"weight sparsification\" or compare with low-rank approximation, quantization, and knowledge distillation methods. \n2. The proposed method is compared only with AGP[5] on PTB. The authors state that this area has problems with benchmarks and baselines. I agree that it makes the fair comparison challenging. On the other hand, I think the results on a single dataset against a single baseline from 2017 are not sufficient to support such a strong claim. I would consider at least one more dataset and compare with, e.g., GraNet[6], and some sparsification methods for RNNs: e.g., [7, 12]. Otherwise, please revise the claims and contributions.\n3. (Minor) All experiments are entirely deferred to the appendix. If the authors declared this as one of the contributions, I would expect it in the main paper. \n\n**Related work** \n\n* **Uncertainty quantification:** please add DropConnect[3] and deep ensembles[4] as uncertainty estimation methods and consider them for comparison in the multi-horizon setting.\n* **Sparse DNN** should be largely extended by discussing different before-training [8], post-training [9] and during-training [5,6,7,10,11] approaches. Regarding the RNN sparsification, I would consider citing and discussing the following works [12, 13, 14]. \n\n**Clarity**\n\nIt was hard to follow some theoretical parts of the paper due to large equations and overloaded notation, e.g., Theorem 3.8. I can see that the authors put effort into simplifying the notation and supporting the theory with discussions in some places. Nevertheless, I think it would be great to do an extra pass to make reading more clear. \n\n(Minor) Figure 1 takes a lot of space but just aims to explain M_l and R_l notation and the RNN setting. I would suggest making it more informative or maybe consider moving it to the appendix. The model compression results arguably seem more important than this figure. \n\n---\n\n[1] Consistent sparse deep learning: Theory and computation. JASA, 2021\n\n[2] Sparse Deep Learning: A New Framework Immune to Local Traps and Miscalibration, NeurIPS 2021\n\n[3] Dropconnect is effective in modeling uncertainty of bayesian deep networks 2019\n\n[4] Simple and scalable predictive uncertainty estimation using deep ensembles NIPS2017\n\n[5] To prune, or not to prune: exploring the efficacy of pruning for model compression, 2017\n\n[6] Sparse training via boosting pruning plasticity with neuroregeneration, NeurIPS2021 \n\n[7] Bayesian Compression for Natural Language Processing, EMNLP2018\n\n[8] SNIP: Single-shot network pruning based on connection sensitivity, ICLR2019\n\n[9] The Lottery Ticket Hypothesis: Finding Sparse, Trainable Neural Networks, ICLR2019\n\n[10] Learning Sparse Neural Networks through L0 Regularization, ICLR2018\n\n[11] Scalable training of artificial neural networks with adaptive sparse connectivity inspired by network science, Nature 2018\n\n[12] Spectral Pruning for Recurrent Neural Networks, AISTATS2020\n\n[13] Structured Sparsification of Gated Recurrent Neural Networks, AAAI2020\n\n[14] Stage-Wise Magnitude-Based Pruning for Recurrent Neural Networks. TNNLS2022.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review. We provide a point-by-point response to your comments below.\n\n$\\textbf{W1: Novelty}$\n\nThe seminal work [1] has established a general theoretical framework for studying the asymptotic behavior of posterior distributions and Bayesian estimators for high-dimensional statistical models. The posterior consistency theory established therein requires three conditions for the models and the prior given in section B of the appendix, lines 633-636. In words, these conditions require:\n(a) The set of models considered in the analysis cannot be too large (in terms of packing numbers).\n(b) The set of models almost encompasses the support of the prior distribution.\n(c) The prior distribution places a sufficient amount of mass near the true model.\nHowever, verifying conditions (a)-(c) for a specific problem is not trivial. For example, [2] addressed the issue for high-dimensional generalized linear models; [3] tackled the issue for deep neural networks with i.i.d. data; and [4] dealt with the issue for nonparametric Bayesian models. In this paper, we address the issue for RNNs with $\\alpha$-mixing sequences.\nAlthough we employ similar conditions as in [3] to constrain the set of neural networks used in data modeling, verification of conditions (a)-(c) remains non-trivial. This is detailed in the Appendix of the paper. Notably, some inequalities used for i.i.d. data no longer hold, such as the inequalities on $d_t$ difference (a generalization of the Kullback-Leibler divergence) used for proving Proposition 1 in [2]. To our knowledge, this paper provides the first theoretical study of sparse RNNs for time series data from the perspective of structure selection, parameter estimation, and prediction uncertainty quantification. Given the wide usefulness of RNNs for time series data [11] , we believe that their theory deserves special attention.\n\n$\\textbf{W2: Algorithm}$\n\nYes, the prior annealing algorithm is a direct application of the algorithm in [5]. Our main contribution is extending the theoretical results to RNN models for time series data. Additionally, our experiments demonstrate the advantages of the proposed method over prior works.\n\n$\\textbf{W3: Model sparsification and baselines}$\n\nThank you for the suggestion. In the revision, we will change the term `compression` to `sparsification` as suggested. We have conducted an additional experiment comparing our model sparsification approach to the baseline method proposed in [6]. Following their experimental setup exactly, we trained a 1 layer RNN with 128 hidden units on PTB data, using the same batch size, number of epochs, and 5 independent runs. As shown in table below, our approach achieves better test perplexity than [6] under similar (and higher) sparsity levels. All baseline results are directly adopted from their paper. \n\n|   Methods   | Test Perplexity  | Sparsity |\n| ----------- | ----------- | ----------- |\n| Baseline     |114.66 (0.35)  | 0%  |\n| $\\textbf{PA (ours)}$  | 117.80 (0.10)  | 70%  |\n| Spectral w/ rec [6] | 124.26 (0.39)  | 67%  |\n\nWe will add both the results and details of these additional experimental results in the revised version.\n\n$\\textbf{W4: Related work: uncertainty quantification}$\n\nFor the multi-horizon forecasting experiments, we followed the setup from [7] and used an LSTM as the base prediction model. Per your suggestion, we will add discussion of deep ensembles [8] and DropConnect [9] to the related work section. Those methods were originally developed for i.i.d. data like image classification. Therefore, adapting them to RNNs for multi-horzion time series forecasting would require careful tuning and an extension of their original methods with the Bonferroni correction to ensure a fair comparison.\n\n$\\textbf{W5: Related work: sparse DNN}$\n\nThe lottery ticket hypothesis [10] shows that for many vision tasks, there exist sparse networks that can be trained from scratch to achieve good performance, but the models must have special initialization. The work shed light on research on pruning before training. During training approaches typically add regularization during training to force parameters to go to 0. Post-training approaches operate on trained neural networks and attempt to remove network parameters based on some pruning criteria such as parameter magnitude, Hessian of the loss function, etc. We will add discussion of these works in the revision. From this perspective, our work is mostly aligned with the during-training approach.\n\n$\\textbf{W6: Clarity}$\n\nIn the revision, we will continue our efforts to simplify the notation and provide additional discussions for the theory to enhance the readability of the paper.\n\n$\\textbf{W7: Figure 1}$\n\nIn the revision, we will move Figure 1 to the appendix and add more model compression results to the main body of the paper.\n\n[1] Convergence rates of posterior distributions. Annals of Statistics, 2000.\n\n[2] Bayesian variable selection for high dimensional generalized linear models: convergence rates of the fitted densities. The Annals of Statistics, 2007.\n\n[3] Consistent sparse deep learning: Theory and computation. JASA2021.\n\n[4] Nonparametric bayesian model selection and averaging. 2008.\n\n[5] Sparse deep learning: A new framework immune to local traps and miscalibration. NeurIPS2021.\n\n[6] Spectral pruning for recurrent neural networks. AISTATS2020.\n\n[7] Conformal time-series forecasting. NeurIPS2021.\n\n[8] Simple and scalable predictive uncertainty estimation using deep ensembles. NeurIPS2017.\n\n[9] Dropconnect is effective in modeling uncertainty of bayesian deep networks. Scientific reports, 2021.\n\n[10] The Lottery Ticket Hypothesis: Finding Sparse, Trainable Neural Networks, ICLR2019.\n\n[11] Deep learning for twelve hour precipitation forecasts. Nature communications, 13(1):1–10, 2022."}, {"review_id": "QAPYTMCB2Z", "reviewer": "Reviewer_TwH7", "summary": "This paper proposes to extend sparse deep learning in the context of time series data. This context is different from the classical one as samples are not *i.i.d.* but dependent. The paper shows theoretical results -- posterior consistency and asymptotic normality of the weights --, a computation method based on posterior annealing and some numerical experiments that showcase the performance of the method on real datasets.", "questions": "- **Q1** - The theoretical results are very similar to the ones of [17], which considers *iid* data with DNN. As the proposed method consider the input of the network as fixed segments $x_{i:i-M_l+1}$, the considered networks can be assimilated as weight-tied DNN, so the class of architecture is not really different. And in `Theorem3.8`, it is not clear where the dependence between the samples intervenes. In particular, I am surprised that we get similar sample efficiency rates between the two settings, and I wonder what is lost by dropping the iid assumption. Could the authors comment on this point?\n- **Q2** -  Many details are deferred to the appendix or referred to in other articles. This makes the paper sometimes hard to read. In particular, `Section.4` is very hard to follow, in particular the part on the transformation of the Bayesian method into a regularization method. As this constitutes the core of the proposed method, I think this could be better explained in the main article.\n- **Q3** - In equation (4), it is unclear what are the conditions for  $\\varpi_n$ so that the argmin is not empty. In particular, as the size of the neural network in the set $\\mathcal G_n$ is constrained in `Assumption3.4`. This specific assumption that there exists a minimizer is used in Eq. bellow l.657, which relies on the fact that for all $n$, $\\|\\mu(X, \\beta, \\gamma) - \\mu^*\\|_{L^2(\\Omega)} \\le \\varpi_n$ is not empty. This is not guaranteed I think, in particular as $\\varpi_n$ goes to 0 at the same rate as the network size increases, while the data point to interpolate increases much faster. Note that a similar question arises in the proof for [17]. It seems that the results in [17] show that universal approximation holds for a specific class of functions even for sparse DNN, but it is unclear if the constraints in `Assumption3.4` allows for this. But maybe this is due to a misunderstanding on my side.\n- **Q4** The method relies on many hyperparameters, that are specified in the appendix. However, the process for their selection is not reported, and it would be nice to understand how sensitive the hyperparameter selection of the model is.\n\n**Minor remarks and typos:**\n- **m1** - `Section 5.2` - Isn't it possible to estimate the AR order of window size =1 by looking at the average dependency length? Or at least estimating if changing one sample changes the result?\n- **m2** - The proposed results seem to relate to the compress sensing and sparse coding theory, it would be nice to draw the connection between these two fields in the introduction. But I have no specific results in mind, so maybe it is too far-fetched.\n\n*unnumbered remarks don't call for answers*\n- Eq. (2) - Missing $+ v^{H_n}z^{H_n}_{i-1}$ at the end to show the structure that is not just feedforward?\n- l.123 - $\\gamma = \\{\\gamma_i \\in \\{0, 1\\} : j=1\\dots K_n\\}$ -> $\\gamma \\in \\{0, 1\\}^{K_n}$?\n- l.155 - an RNN of size ... has been\" -> \"... is\".\n- l.192 - mismatch between name and acronym `MIPP` and `marginal posterior inclusion probability`.\n- l.218 - $H_n$ is already the number of layers (l.113), finding a different notation would help the reader follow the arguments.\n- l305 - \"In particular, Our method\" -> \"our\".", "limitations": "* see weaknesses.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The paper extend sparse deep learning result to the challenging context of dependent data, applicable for time series.\n- The many reported results show improvement of the method compared to existing state-of-the-art methods and seems to yield interpretable models for time series.", "weaknesses": "- It is not totally clear what is really different between the results in [17] and the results presented here (see **Q1**).\n- The writing is not always clear (**Q2**), in particular for non-expert in sparse deep learning. It is not clear why `Eq.(4)` admits an non-empty set of minimizers (**Q3**).\n- There are many hyperparameters in the methods and while there value is reported, it is not clear how to select them (**Q4**).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review. We provide a point-by-point response to your comments below.\n\n$\\textbf{W1 and Q1: Novelty}$\n\nThe seminal work [1] established a general theoretical framework for studying the asymptotic behavior of posterior distributions and Bayesian estimators for high-dimensional statistical models. The posterior consistency theory established therein requires three conditions for the models and the prior, namely: (a) The set of models considered in the analysis cannot be too large (in terms of packing numbers).\n(b) The set of models almost encompasses the support of the prior distribution.\n(c) The prior distribution places a sufficient amount of mass near the true model. \nThe work [2] further extended this theory to non-i.i.d data under the three general conditions, leading to a similar posterior contraction rate. Dropping the i.i.d assumptions disables many inequalities used for the i.i.d case, necessitating specific considerations. However, the resulting contraction rates are not necessarily very different. In the same vein, we get a similar posterior contraction rate as reference [4] (cited in the paper).  Please also see our responses to reviewer AFD6 (W1)\n\n$\\textbf{W2 and Q2: Transformation of the Bayesian method}$\n\nIn the revision, we will add a lemma regarding applications of the Laplace approximation to Bayesian neural networks. Similar to Theorem 2.3 in [4], the main idea is to connect posterior consistency to the consistency of the maximum a posteriori (MAP) estimator. Then the prior annealing algorithm is employed for finding the optimum of the posterior.\n\n$\\textbf{W3 and Q3: Equation 4 and approximation error}$\n\nThis is a very thoughtful question. Lemma 4.1 of [3], which, through the trick of independent block sequence construction, shows that many properties of the i.i.d processes can be extended to mixing processes. While the lemma was proven for the case of $\\beta$-mixing, the author did mention her doubts about its applicability to $\\alpha$-mixing. Therefore, at least for the sequences of $\\beta$-mixing (which implies $\\alpha$-mixing or can be viewed as a subclass of $\\alpha$-mixing), the non-empty of the sparse RNN set can be guaranteed for many classes of functions. This issue will be mentioned in the revision.  Moreover, the contraction rate of our posterior consistency results is essentially a sum of the approximation error and estimation error, i.e. $\\epsilon_n^2 = O(\\varpi^2_n) + O(\\zeta_n^2)$ in Theorem 3.8, which provides the flexibility to be combined with other approximation theories to give get exact order of $\\epsilon_n^2$.\n\n$\\textbf{W4 and Q4: Hyperparameters}$\n\nThank you for your feedback. Regarding the four hyperparameters related to the mixture Gaussian priors $(\\lambda_n, \\sigma_{1,n}, \\sigma_{0,n}^{init}, \\sigma_{0,n}^{end})$, our algorithm is not sensitive to $\\lambda_n$ and $\\sigma_{0,n}^{end}$. $\\lambda_n$ is only used for adjusting the target sparsity if needed, and it is usually selected from {1e-6, 1e-7}. $\\sigma_{0,n}^{end}$ is generally set between {1e-5, 1e-6, 1e-7}, with the exact value depending on $\\sigma_{0,n}^{init}$. We gradually decrease $\\sigma_{0,n}^{init}$ to $\\sigma_{0,n}^{end}$ during prior annealing. Based on the number of iterations/epochs and the learning rate used in this stage, the difference between $\\sigma_{0,n}^{end}$ and $\\sigma_{0,n}^{init}$ should be adjusted so that sparsity remains stable.\nSelecting $\\sigma_{0,n}^{init}$ depends on the specific task. For example, as explained in Section E of the appendix, we select $\\sigma_{0,n}^{init}$ by achieving the target sparsity for model sparsification tasks. For other tasks like model selection and uncertainty quantification, $\\sigma_{0,n}^{init}$ is not as sensitive as long as sparsity is stable during annealing. Our algorithm is also robust to the temperature hyperparameter. We will add a section in the appendix instructing hyperparameter selection in the revised version.\n\n$\\textbf{Minors 1:}$\n\nThank you for your feedback. Through our selected network structure, the AR order will not be directly available in the case \nwindow size=1. However, in practice, it is still possible to estimate the AR order by looking at the average dependency length. This point will be mentioned in the revision. \n\n$\\textbf{Minors 2:}$\n\nThank you very much for your valuable suggestions! We will review compress sensing and sparse coding theory, and add them accordingly in the introduction.\n\n$\\textbf{Typos}$\n\nThank you for your feedback. We will fix these typos in the revised version.\n\n\n[1] Convergence rates of posterior distributions. Annals of Statistics, 2000.\n\n[2] Convergence rates of posterior distributions for non-i.i.d. observations. Ann. Statist. 35(1), 192–223\n\n[3] Rates of convergence for empirical processes of stationary mixing sequences. The Annals of Probability, 22, 94-116. \n\n[4] Consistent sparse deep learning: Theory and computation. JASA2021."}, {"review_id": "RQjJUvrzdR", "reviewer": "Reviewer_UrAN", "summary": "For iid data, sparse DL has been shown as a way towards consistency of input-output mappings and well understood distribution of model predictions. This theory is missing for time-series data however, which the authors thus introduce here. In particular, the following results for RNNs with Gaussian mixture parameters are presented:\n\n- Consistency of posteriors, structure selection and input-output mappings\n- (asymptotic) normality of predicted values", "questions": "\nMinor things:\n\nThere could be a bit more discussion about mixture priors: if and why they are a reasonable choice\n\nConsider adding a comma around in line 117 to make the meaning of \"only\" and \"or\" more clear in sentence \"the input $x_i$ can contain only $y_{i-1}$ or $y_{i-1:i-r}$ for some r>1\". Now it's not clear if only refers to \"only either of a or b\", or \"only a, or if not a then b is the only alternative\".\n\nI found it  confusing that $H$ was used to denote the layers and $L$ the number of hidden neurons. I would have found them more logical the other way around.\n\nAs a misc. note: I found the font of the paper odd, but that could be a local issue on my side.\n", "limitations": "In general limitations are discussed fairly well --  especially given the limited space. As always, this type of theoretical works is limited by assumptions that are made in particular. Here e.g. Gaussian mixture prior is assumed and I would be curious to hear more about the limitations of this choice of prior in more detail. What alternatives are there? The mixing assumptions can be limiting too, but the authors do acknowledge that much.\n\nAs mentioned above, it would also be nice to have bit more discussion on the assumptions and conditions of Theorems 3.8 and the other consequent theoretical results.\n\n", "rating": 7, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper was well organized and written. \n\nExtending the recent results from sparse DL work to time-series is valuable and a novel contribution.\n\nThe consistency results presented here are convincing and valuable (at least to this reviewer who is not an expert in this field) and should pave way for a lot further work in the area -- the assumptions made in the theoretical parts are fair for this nascent area of study and provides a good starting point.\n\n\nThe authors show how the method can be framed as regularization in the way of laplace approximation at the MAP estimator which is useful in practice. \n\nThe experiments provide evidence of superior performance in comparison to CP methods.\n", "weaknesses": "Comparison to other possible approaches to uncertainty quantification is lacking: there is no comparisons or discussion of this method to the alternatives mentioned in the introduction: \"include multi-horizon probabilistic forecasting [33], dropout-based methods [34], and recursive Bayesian methods [35].\". This is the case with the experiments, but also on theoretical side. Naturally theory for the other methods may not exist in the same vain, but it would have been useful if the authors tried to make some further comparisons. \n\nOne limiting assumption here is that $\\mu*{\\star}$ can be approximated by a sparse RNN arbitrarily well (i.e. ground-truth sparse model).   Especially since the connectivity of this *true* RNN model is assumed to be limited. Also, as the authors point out, universal approximation does not hold fully (but can still hold for some nonlinear function classes).  Relatedly, theorem 3.8 relates to the estimation error only, but what about the approximation error?\n\nThe assumptions and conditions required for Theorems 3.8 onward are not discussed in sufficient depth.\n\nThe AR order selection experiments, while promising, are quite narrow (only 5 datasets from a single model). More different datasets would be needed with different settings (of window sizes etc.)\n\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review. We provide a point-by-point response to your comments below.\n\n$\\textbf{W1: Comparison to other uncertainty quantification approaches}$\n\nWe apologize for the confusion. We actually included MQ-RNN [1] (multi-horizon probabilistic forecasting) and DP-RNN [2] (dropout-based method) as the two baselines for the set of time series experiments. We will make sure the baselines used in the experiments are clear by adding the corresponding citations in the table.\n\n$\\textbf{W2: Approximation error}$\n\nApologies for the confusion. The convergence rate of Theorem 3.8 actually depends on the approximation error. It is essentially a sum of the approximation error and the estimation error, i.e., $\\epsilon_n^2 = O(\\varpi^2_n) + O(\\zeta_n^2)$ in line 173 of the paper. Please also see our reply to W3Q3 to reviewer TwH7.\n\n$\\textbf{W3: Discussions about Theorems 3.8}$\n\nThank you for suggestions. We will add more discussions about the assumptions and conditions. Additional assumption for theorem 3.9 is essentially an identifiability assumption, which is commonly used for variable selection. For theorem 3.11 and 3.12, additional assumptions and discussions are given in section C in appendix.\n\n$\\textbf{W4: Additional AR order selection experiments}$\n\nThank you very much for your suggestion! As suggested, we have conducted additional AR order selection experiments with more window sizes (6 window sizes) using the following exponential autoregressive order model from [5].\n\n$$y_i = \\left(0.8-1.1\\exp\\\\{-50y_{i-1}^2\\\\}\\right)y_{i-1} + \\eta_i$$\n\nwhere where $\\eta_i \\sim N(0,1)$ are i.i.d. Gaussian random noises. We follow similar settings as section 5.2, i.e., generated 5 datasets, with the training sequence having a length of 10000, and both the validation and test sequences having a length of 1000. The results are given in the table below.\n\n|   Model   | Window size | FSR  | NSR | AR order | #hidden link | MSPE | MSFE |\n| ----------- | ----------- | ----------- | -----------| ----------- | ----------- | ----------- | ----------- |\n| PA-RNN | 1   | 0 | 0 | 1(0) | 0(0) | 1.003(0.005) | 1.004(0.004) |\n| PA-RNN | 3   | 0 | 0 | 1(0) | 0(0) | 1.006(0.005) | 0.999(0.004) |\n| PA-RNN | 5   | 0 | 0 | 1(0) | 0(0) | 1.007(0.005) | 1.000(0.004) |\n| PA-RNN | 7   | 0 | 0 | 1(0) | 0(0) | 1.006(0.005) | 1.000(0.003) |\n| PA-RNN | 10 | 0 | 0 | 1(0) | 0(0) | 1.002(0.006) | 1.002(0.004) |\n| PA-RNN | 15 | 0 | 0 | 1(0) | 0(0) | 1.002(0.007) | 1.001(0.004) |\n\nThis nonlinear autoregressive order model has an AR order of 1. We use the same hyperparameters for all window sizes, and our method managed to achieve perfect results for all window sizes with respect to all model selection metrics. This might be because this is a relatively simpler model compared with the one we used in Section 5.2. We will add both the results and details of these additional experimental results in the revised version.\n\n$\\textbf{Q1: Discussion about mixture priors}$\n\nThe mixture Gaussian prior can be viewed as a continuous relaxation of the spike-and-slab prior, which puts a point mass at 0. Compared to the spike-and-slab prior, a significant advantage of the mixture Gaussian prior is that it is differentiable, allowing stochastic gradient-based algorithm to be used for computation. And the mixture allows us to put enough weight on the neighborhood of 0, which essentially allows the prior to satisfy conditions (a)-(c) in section B of the appendix. \n\n$\\textbf{Q2: Line 117}$\n\nThank you, we will add a comma in the revised version to clarify the meaning as you suggested.\n\n$\\textbf{Q3: Notation H}$\n\nWe sincerely apologize for the confusion. We mainly followed the notations used in [3,4], but we will make sure it is more clear in the revised version.\n\n$\\textbf{Q4: Font of the paper}$\n\nWe will have this issue addressed in the revision.\n\n[1] A multi-horizon quantile recurrent forecaster. arXiv preprint, 2017.\n\n[2] Dropout as a bayesian approximation: Representing model uncertainty in deep learning. ICML2016.\n\n[3] Consistent sparse deep learning: Theory and computation. JASA2021.\n\n[4] Sparse deep learning: A new framework immune to local traps and miscalibration. NeurIPS2021.\n\n[5] Identification of nonlinear time series: First order characterization and order determination. Biometrika, 77(4):669–687, 1990."}, {"review_id": "HslTdCrX9X", "reviewer": "Reviewer_riun", "summary": "This paper focuses on the theoretical analysis of sparse deep learning to time series data.  Statistical propoerties of sparse RNNs are investigated including consistency and asymptotical behaviour.  The paper presents some numerical results showing that sparse deep learning outperforms existing methods in predicting uncertainty quantification for time series data, as well as in identifying the autoregressive order for time series data and large-scale model compression. The proposed method has practical implications in some fields, such as finance, healthcare, and energy.", "questions": "(1)The authors assume that the time series are $\\alpha$-mixing process, but in fact there are many other mixing processes, such as $\\beta$-mixing and $\\tau$-mixing. Please explain how the $\\alpha$-mixing affects the theoretical analysis of algorithms. Whether the different mixing will cause different results?\n\n(2)Whether the data used in the experiment satisfy the $\\alpha$-mixing? If yes, please provide some supporting materials. Some auto-regressive models are not $\\alpha$-mixed, e.g.,\\begin{equation}\\label{1}\n\\mathbf{x}_t=f_0(\\mathbf{x}_{t-1})+\\epsilon_t,~t\\in \\mathbb{Z},\n\\end{equation},\nwhere $f_0:[-K,K]^d\\rightarrow[-K,K]^d$ is a Lipschitz function with Lipschitz constant $\\leq$ 2 and $\\epsilon_t\\sim\\mathcal{B}(0.5)$ [4].\n\n(3)Please give the proof of Lemma B.1 for completeness.\n\n(4)What is the biggest difference between the algorithm the authors proposed in section 4 and the algorithm proposed in [5]？\n\n(5)This article is too similar to [5] and [2]. It seems that the result of this paper is an extension from i.i.d. observations to dependent data. Please clarify what are the main difficulties in the theoretical analysis for time series data compared with i.i.d. data.\n\n[2]Yan Sun, Qifan Song, and Faming Liang. Consistent sparse deep learning: Theory and computation. Journal of the American Statistical Association,  2021.\n[4]D. W. Andrews. Non-strong mixing autoregressive processes. Journal of Applied Probability, 1984.\n[5]Yan Sun, Wenjun Xiong, and Faming Liang. Sparse deep learning: A new framework immune to local traps and miscalibration. Advances in Neural Information Processing Systems, 2021.\n", "limitations": "This is a theoretical work. But part of the proof is a little vague such as proof of Theorem 3.11 and Theorem 3.12. Authors should highlight the key differences  of analysis techniques  between their work and [2] [5] (instead of just using the different tools).\n\n[2]Yan Sun, Qifan Song, and Faming Liang. Consistent sparse deep learning: Theory and Computation.\n[5]Yan Sun, Wenjnd computation. Journal of the American Statistical Association,  2021.un Xiong, and Faming Liang. Sparse deep learning: A new framework immune to local traps and miscalibration. Advances in Neural Information Processing Systems, 2021.", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "(1) The paper addresses an important and underexplored topic by studying sparse deep learning for dependent time series data. This fills a gap in existing research, which has mostly focused on i.i.d. data.\n\n(2) The paper presents the empirical results that show the superiority of sparse deep learning over existing methods in predicting uncertainty quantification and autoregressive order identification for time series data. This demonstrates the practical effectiveness of the proposed method.", "weaknesses": "(1)Some of the theorems in references are not given in the text and some symbols are confused, e.g.,\n 1)Please give the concrete content of Theorem 2 from [1] that the authors cite.\n 2)Please give the definition of $O_{P^*}$ on line 174.\n 3)The lemma used on line 575 should be lemma S1 of section 4 in [2]. The modifications of the above weaknesses can enhance the readability of the article.\n\n(2)There are related references that study the theoretical analysis of deep neural networks for temporally dependent observations, such as [3]. I encourage the authors to check them out and see if they should be included in the related work.\n\n[1]Pentti Saikkonen. Dependent versions of a central limit theorem for the squared length of a sample mean. Statistics & probability letters, 1995.\n[2]Yan Sun, Qifan Song, and Faming Liang. Consistent sparse deep learning: Theory and computation. Journal of the American Statistical Association, 2021.\n[3]M. Ma and A. Safikhani. Theoretical analysis of deep neural networks\nfor temporally dependent observations, in NeurIPS, 2022.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review. We provide a point-by-point response to your comments below.\n\n$\\textbf{W1: Readability: Theorems and Lemmas}$\n\nThank you very much for your thorough review of our paper. We sincerely appreciate your feedback, and in the camera-ready version, we will address the mentioned weaknesses as follows:\n\n(i) We will add the concrete content of Theorem 2 from [2] to the text for better clarity and reference.\n\n(ii) $O_{P^{*}}(a_n)$ is used to represents a sequence of random variables bounded in probability. Let $\\\\{X_n\\\\}$ be a sequence of random variables, and let $\\\\{a_n\\\\}$ be a sequence of strictly positive reals. We say $X_n/a_n$ is bounded in probability, if for every $\\epsilon > 0$, there exists $M_\\epsilon > 0$, such that \n\n$$P^{*}(|X_n/a_n| > M_\\\\epsilon) < \\\\epsilon$$\n\nfor all $n$, and we can write it as $X_n = O_{P^{*}}(a_n)$. We will add this definition in the revision.\n\n(iii) You are correct; the lemma used on line 575 should indeed refer to lemma S1 of section 4 in [4]. We will make this correction in the revised version of the paper.\n\n$\\textbf{W2: Related references}$\n\nThank you for providing such useful information! We greatly appreciate your suggestion, and we will add the work [3] to the related work section for discussion. Upon careful reading of [3], we observed that their focus is on establishing consistency rates for prediction error in the context of temporally dependent observations. In contrast, our work takes a Bayesian approach, which offers distinct statistical properties, to address a similar but inherently different problem. Please also see our responses to reviewer AFD6.\n\n$\\textbf{Q1: Different mixings}$\n\nOur results apply to other types of mixing sequences as implied by the following hierarchy of the five mixing conditions: \n(i) $\\psi$-mixing implies $\\phi$-mixing;\n(ii) $\\phi$-mixing implies both $\\rho$-mixing and $\\beta$-mixing; \n(iii) $\\rho$-mixing and $\\beta$-mixing each impliy $\\alpha$-mixing.\n\nSee reference [1] (Theorem 3.11) for theoretical justifications for the above hierarchy: \n\n$\\textbf{Q2: Data with mixing property}$\n\nThank you very much for your suggestion. We use those datasets in order to fairly compare with current state-of-the-art methods. We conducted additional AR order selection experiments using the following exponential autoregressive order model from [5].\n\n$$y_i = \\left(0.8-1.1\\exp\\\\{-50y_{i-1}^2\\\\}\\right)y_{i-1} + \\eta_i,$$\n\nwhere $\\eta_i \\sim N(0,1)$ are i.i.d. Gaussian random noises. This model is shown to be $\\alpha$-mixing according to [5], please see our new experimental details and results for this model in our responses to reviewer UrAN (W4).\n\n$\\textbf{Q3: Proof of Lemma B.1}$\n\nWe will give the complete proof of Lemma B.1 for completeness in the revision.\n\n$\\textbf{Q4: Algorithm}$\n\nPlease see our responses to reviewer AFD6 (W2).\n\n$\\textbf{Q5: Novelty}$\n\nPlease see our responses to reviewer AFD6 (W1).\n\n[1] R.C. Bradley (2007) Introduction to Strong Mixing Conditions. Vol. 1. Kendrick Press, Heber City (Utah).  \n\n[2] Dependent versions of a central limit theorem for the squared length of a sample mean. Statistics & probability letters, 1995.\n\n[3] Theoretical analysis of deep neural networks for temporally dependent observations. NeurIPS2022.\n\n[4] Consistent sparse deep learning: Theory and computation. JASA2021.\n\n[5] Identification of nonlinear time series: First order characterization and order determination. Biometrika, 77(4):669–687, 1990."}, {"review_id": "kHzZsCFPAz", "reviewer": "Reviewer_shpG", "summary": "This paper introduces an expansion of sparse deep learning theory, specifically tailored for the analysis of time series data. The authors provide a detailed explanation of the theoretical foundation behind this theory, employing recurrent neural networks (RNNs), particularly long short-term memory (LSTM) networks. The authors demonstrate that their novel approach yields competitive outcomes when compared to other cutting-edge methodologies from diverse research domains, including conformal prediction.", "questions": "Here are some questions/suggestions for the authors:\n\n- In Figure 1, the circles shouldn’t contain some variables or some naming? Are these the states of the RNN? This figure is a bit misleading to me.\n- In line 239, what is ******local-trap issue?:****** A bit more context could help the reader.\n- Is there any reason why you choose LSTM against other architectures? Some study on whether this applies to GRU or other RNN could be interesting.\n- Do you believe that the LSTM is providing good results because of the forget gate? Or do you think there is some alignment between your results and the use of LSTM, that might actually shadow the developed theroy?\n- In Table 3, when considering the autoregressive order selection, it appears that the errors of the proposed model are comparable to those of the baselines that only employ RNN. Notably, the absence of standard errors in the results raises some questions. It is possible that I may have overlooked something, or perhaps the authors could provide additional clarity regarding the significance of these findings. One noteworthy aspect is whether the reduced number of hidden links is the primary factor contributing to the performance of the proposed method.\n\nMinor questions:\n\n- What is PI length? Prediction Interval length?  It is just presented in the table with no extra explanation.\n- Why the method is called PA? When do you choose the name?", "limitations": "The authors did not explicitly address the limitations of their work, which raises concerns about the feasibility of deploying their proposed model in practical settings. It is crucial to understand the potential challenges and constraints associated with implementation. Considering the computational time and power required, it is essential to evaluate the trade-off involved and compare it to the relatively inexpensive nature of conformal prediction. Understanding these aspects would provide valuable insights into the practical implications and potential drawbacks of the proposed model. Plus, I lack some direct comments on the limitations with respect to conformal prediction, since my prior belief is that the their presented comparison w.r.t conformal prediction methods is scarce. \n\nWith all this in mind, I consider that the paper needs some further improvement, and I await for the author’s and other reviewers’ comments for the possibility of increasing the score.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well-written and grammatically correct with no evident typos. The narrative is clear: related work, theoretical grounding, computation and experiments. \n\nIt is clear that the authors have well investigated the theory of sparse deep learning and it is evident that they know it quite well the field. Hence, the derivations, in principle, seem correct — although I must admit that I found the ideas and derivations presented in the paper to be somewhat challenging to grasp, even though I devoted a considerable amount of time to understanding them. Plus, newcomers to the sparse deep learning theory can find the complete derivations in the Appendix, which is quite extensive.\n\nThe proposed method is tested in three relevant tasks related to time series data, and it is evaluated on at least four real-world datasets. This thorough evaluation approach provides strong empirical evidence to support the theoretical foundations of the method, demonstrating its validity and applicability.", "weaknesses": "Before outlining the weaknesses, I repeat again that I am a new reader of sparse deep learning theory. However, I still find that the paper might exhibit some weaknesses:\n\n1. First of all, while I believe the theoretical side of the paper is of good quality, I encountered difficulties in understanding the concepts due to the overwhelming presence of complex notation. Sometimes there are many variables, and many subscripts, which in temporal series can be, e.g. $y_{k-1:k-R_l}$. This makes the reading difficult and sometimes cumbersome.\n2. After reading the paper, I still cannot see why do you choose conformal prediction as the main comparison field. I lack some more context or motivation on the paper. The paper might be already interesting for the natural extension of sparse deep learning for time series data, but I also believe there should be some kind of justification of why these methods are a good alternative. \n3. The point above leads me to the following point: I worry that the literature review for Conformal Prediction might be scarce. Authors comment on the work by [1], but there is no reference to other works in the same line of this work [3], or other contemporary works such as [2]. \n4. And following this line, I wonder if this theory might be very challenging to apply in real-world scenarios compared to CP. The beauty of CP is that you can basically apply it on top of any ML model in a “plug-and-play” fashion, with the cost of very strong assumptions such as the ****************exchangeability**************** assumption. But it is more straightforward to apply and to easily understand than this sparse deep learning technique. I would like to know what the authors think about this.\n\n[1] Isaac Gibbs and Emmanuel Candes. Adaptive conformal inference under distribution shift. Advances in Neural Information Processing Systems, 34:1660–1672, 2021.\n\n[2] Barber, R. F., Candes, E. J., Ramdas, A., & Tibshirani, R. J. (2023). Conformal prediction beyond exchangeability. *The Annals of Statistics*, *51*(2), 816-845.\n\n[3] Romano, Y., Patterson, E., & Candes, E. (2019). Conformalized quantile regression. *Advances in neural information processing systems*, *32*.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you so much for your review. We provide a point-by-point response to your comments below.\n\n$\\textbf{W1: Complex notation}$\n\nIn the revision, we will simplify the notation and add more explanations to the assumptions and results, making the paper more accessible to the readers. \n\n$\\textbf{W2: Conformal prediction baselines}$\n\nWe chose conformal prediction as the main comparison field because conformal prediction methods have gained significant popularity and adoption as a widely-used approach for uncertainty quantification in various domains, including time series data. By selecting CP as the main comparison field, our aim is to assess the performance of our proposed method against established and widely-recognized approaches for uncertainty estimation.\n\nIn our experiments, we compared the proposed method with two state-of-the-art CP methods for time series data [1,4]. As suggested, we have included another recently published conformal baseline focused on time series (EnbPI V2) [5] for comparison in the single time series experiment (please see table below). We greatly appreciate the reviewer's input on the literature review. In the revised manuscript, we will provide a more comprehensive review, including references to [2,3].\n\nIt's worth noting that [1] presents an adaptive version of [3] specifically designed for data with potential distribution shifts, such as time series data. Additionally, [2] introduced a random swapping mechanism to address potentially non-exchangeable data. Their extension enables conformal prediction to be applied on top of a model trained with weighted samples. Their primary focus is providing theoretical justification of the gap in the coverage rate for the proposed method.\n\n|   Methods   | Coverage | Average PI Lengths (standard deviation)  | Median PI Lengths (interquartile range) |\n| -------- | -------- | -------- | -------- | \n| PA offline | 90.06 | 20.19(3.23)  | 20.59(1.82) |\n| AgACI online  | 91.05 | 23.27(14.49)  | 21.38(2.62) |\n| ACI $\\gamma=0.01$  online | 90.21  | 22.13(8.50)  | 20.68(2.70) |\n| ACI $\\gamma=0$ online | 92.89  | 24.39(9.17)  | 22.99(2.91) |\n| EnbPI V2 online [5] | 91.23  | 27.99(10.17)  | 26.55(5.09)|\n\n$\\textbf{W3: Our method applied in real-world scenarios compared to CP methods}$\n\nWe appreciate the reviewer's comments . While basic CP methods like split conformal [6] are easy to apply on top of any ML models to generate prediction sets with marginal coverage guarantee, the width or the usefulness of the prediction sets will depend on how well the ML model performs or how the non-conformity scores are defined. Our theoretical results provide convergence to the true model, so ML models will perform well under our assumptions. Moreover, advanced CP methods like [4] and [5] can require multiple runs and ensemble models, making them computationally intensive, especially for deep neural networks. In contrast, our sparse deep learning technique offers several advantages in real-world scenarios. Once the model is trained with the mixture Gaussian priors, which act as an alternative to L1 or L2 regularization, the prediction intervals become a byproduct of the training process. This means that no extra training or heavy computation is needed to produce prediction intervals.\n\n$\\textbf{Q1: Figure 1}$\n\nWe apologize for the confusion. The circles in Figure 1 represent the hidden states of a multi-layer RNN. To avoid ambiguity, we will include explicit labeling of the variables in the camera-ready version to enhance clarity.\n\n$\\textbf{Q2: What is local-trap issue?}$\n\nApologies for the confusion. In this context, the term 'local-trap' refers to the same phenomenon as 'local minima.' There has been substantial research examining the favorable properties of the optimization landscape for over-parameterized neural networks, e.g., [7], which aims to explain why gradient-based methods work well for highly non-convex neural network optimization problems. However, for our model, the mixture Gaussian prior introduces heavy penalties on parameters near $0$, so it is unclear how the optimization landscape would be affected and local traps could become an issue. That is the motivation behind applying the prior annealing approach.\n\n$\\textbf{Q3: Is there any reason why you choose LSTM?}$\n\nWe chose LSTM as the main network in our experimental setup to facilitate fair comparisons with current state-of-the-art methods for tasks like time series forecasting and NLP. Many of these existing methods employ LSTM, and by using the same network architecture, we can better assess the effectiveness of our proposed method relative to these benchmarks.\n\n$\\textbf{Q4: Table 3}$\n\nThank you for your feedback. In our table, we reported the corresponding standard deviations for all metrics except FSR and NSR. The reason for the absence of standard deviations for FSR and NSR is that their calculation requires variable selection results from all five datasets, making it infeasible to compute standard deviations for these metrics. If you require further clarification on any aspect, please feel free to let us know. \n\n$\\textbf{Q5 and Q6: PI and PA:}$\n\nYes, you are correct. \"PI length\" stands for Prediction Interval length. To clarify this in the revised version, we will add a brief explanation of \"PI length\" in the table caption. We chose the name 'PA' (Prior Annealing) for our method because a crucial aspect of our algorithm involves annealing the mixture Gaussian prior.\n\n[1] Adaptive conformal inference under distribution shift. NeurIPS2021.\n\n[2] Conformal prediction beyond exchangeability. The Annals of Statistics, 2023.\n\n[3] Conformalized quantile regression. NeurIPS2019.\n\n[4] Adaptive conformal predictions for time series. ICML2022.\n\n[5] Conformal prediction interval for dynamic time. ICML2021.\n\n[6] Distribution-free predictive inference for regression. JASA2018.\n\n[7] The loss surface of deep and wide. ICML2017."}]

Author rebuttal:

Meta-review: There was initially quite some uncertainty around the contribution of this submission and reviewers struggled with the limited clarity of the original exposition (partly also due to complex notation, etc). We strongly encourage the authors to work on the clarity of the exposition as suggested by the reviewers for the final version. The authors went through great length to generate additional results and provide clarifications in the rebuttal. We strongly encourage the authors to include these clarifications and additional results in the final version. After some back and forth, the reviewers were largely convinced (albeit still with relatively low certainty) of the relevance and novelty of this submission and overall I think it is an adequate contribution to NeurIPS.
Decision: Accept (poster)