Finite Sample Analysis Of Dynamic Regression Parameter Learning
Authors: Mark Kozdoba
Edward Moroshko
Shie Mannor
Koby Crammer
Venue: NeurIPS 2022 Accept (2022)
Keywords: dynamic regression
online regression
variance estimation
kalman filter

Abstract: We consider the dynamic linear regression problem, where the predictor vector may vary with time. This problem can be modeled as a linear dynamical system, with non-constant observation operator, where the parameters that need to be learned are the variance of both the process noise and the observation noise. While variance estimation for dynamic regression is a natural problem, with a variety of applications, existing approaches to this problem either lack guarantees altogether, or only have asymptotic guarantees without explicit rates. In particular, existing literature does not provide any clues to the following  fundamental question: In terms of data characteristics, what does the convergence rate depend on?  In this paper we study the global system operator -- the operator that maps the  noise vectors to the output. We obtain estimates on its spectrum, and as a result derive the first known variance estimators with finite sample complexity guarantees. The proposed bounds depend on the shape of a certain spectrum related to the system operator, and thus provide the first known explicit geometric parameter of the data that can be used to bound estimation errors. In addition, the results hold for arbitrary sub Gaussian distributions of noise terms.  We evaluate the approach on synthetic and real-world benchmarks.

Finite Sample Analysis Of Dynamic Regression Parameter Learning
Mark Kozdoba
Edward Moroshko
Technion, Israel Institute of Technology markkōef technion.ac.
Technion, Israel Institute of Technology edvard.noroshho0@mmal. com
Shie Mannor
Technion, Israel Institute of Technology and NVIDIA Research
shieβee. twchniin......ee ac. il
Koby Crammer
Technion, Israel Institute of Technology kobydee. technion. ac.
Abstract
We consider the dynamic linear regression problem, where the predictor vector may vary with time. This problem be modeled linear dynamical sys- tem, with non-constant observation operator, where the parameters that need to be learned are the variance of both the process noise and the observation noise. While variance estimation for dynamic regression is a natural problem, with a variety of applications, existing approaches this problem either lack guarantees altogether, or only have asymptotic guarantees without explicit rates. In particular, existing literature does not provide any clues the following fundamental question: In terms of data characteristics, what does the convergence rate depend on? In this paper we study the global system operator the operator that maps the noise vee- tors to the output. We obtain estimates on its spectrum, and result derive the first known variance estimators with finite sample complexity guarantees. The proposed bounds depend on the shape of certain spectrum related to the system operator, and thus provide the first known explicit geometric parameter of the data that can be used to bound estimation errors. In addition, the results hold for bitrary sub Gaussian distributions noise terms. We evaluate the approach on synthetic and real-world benchmarks.
Introduction
A dynamic linear regression (West and Harrison. [197), Chapter 3), non-stationary regression, is a situation where we are given sequence scalar observations {Yihes C R, and observation vectors [wtreT C R such that Y, = (X,, wy) + 2, where X, E R is regressor vector, and 2, random noise term. In contrast to a standard linear regression, the vector X, may change with time. One common objective for this problem is at time T, to estimate the trajectory of X, for t < T. given the observation vectors and observations, {wceesT (YthesT and possibly to forecast Yr+1 WT+1 is also known.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

In this paper we model the problem as follows:
Xe+1 Xx + l1
(I)
Yi = (Xe, ue) + St,
(2)
where (.,-) is the standard inner product on R=, =t, the observation noise, are zero-mean sub Gaus- sian random variables, with variance r/², and the process noise variables h2 take values in R", such that coordinates of fyy are zero-mean sub Gaussian, independent, and have variance o². All f, and 2, variables are assumed to be mutually independent. The vectors Ut are an arbitrary sequence in IR, and the observed, known, quantities at time ' are {YeheT and {we)erT
The system (1) a special case of a Linear Dynamical System (LDS). As is well known, when the parameters 01 are given, the mean-squared loss optimal forecast for Y++1 and estimate for Xr are obtained by the Kalman Filter (Anderson and Moore, 1979 Hamioon 1994 Chu and Chen 20T7). In this paper are concerned with estimators for o, 7), and finite sample complexity guarantees for these estimators.
Let us first make few remarks about the particular system (I))).. First, as natural model of time varying regression, this system is useful in a considerable variety of applications. We refer to West and Harrison (I197)), Chapter 3, for numerous examples. In addition, an application to electricity consumption time series as a function of the temperature is provided in the experiments section of this paper. Second, one may regard the problem of estimating o, [) in (I)((Z) as pure case of finding the optimal learning rule for X,. Indeed, the Kalman filter equations for (I)((2), are given by (G)-(E) below, where (E) describes the filtered covariance update and (e) the filtered state update. Here T, is the estimated state, given the observations Yi...., Yt, see West and Harrson (1997)).
n
C++1 ((C, + o2 121 uo+1) + (Co + 021)
(3)
C++1 u+1 (V4+1 - (Tt, U(+1)).
T+41 = Et r/2
(4)
In particular, following (9), the role of a and y may be interpreted as regulating how much the estimate of ēs+1 is influenced, via the operator C441 by the most recent observation and input x²
Yi+1, Mr+11 Roughly speaking, higher values of or lower values of 1) would imply that the past observations are given less weight, and result in an overfit of the forecast to the most recent observa- tion. On the other hand, very low a or high 7) would make the problem similar to the standard linear regression, where all observations are given equal weight, and result lag of the forecast. See Figure B in Supplementary Material Section A for an illustration.
Finally, it is worth mentioning that the system (I) (2) is closely related the study of online gradient (OG) methods (Zinkevich, 2003; Hazan, 2016). In this field, assuming quadratic cost, one considers the update
Fi41 = Fy + a U4+1 (Y2+1 (Tee M(+1));
where is the learning rate, and studies the performance guarantees of the forecaster (āz, u4+1). Compared to (9)) the update (5) is simpler, and uses a scalar rate a instead of the -eput-dependent operator rate C++1/1" of the Kalman filter. However, due to the similarity, every domain of appli- cability of the OG methods is also a natural candidate for the model (I))(2 and vice-versa. As an illustration, we compare the OG to Kalman filter based methods with learned o.1y in the experiments section.
In this paper we introduce a new estimation algorithm for a, 17, termed STVE (Spectrum Threshold- ing Variance Estimator), and prove finite sample complexity bounds for it. In particular, our bounds are an explicit function of the parameters T and {ur) for any finite T, and indicate that the es- timation error decays roughly as - 2,, with high probability. To the best of our knowledge, these are the first bounds of this kind. As we discuss in detail in Section D. most existing estimation meth- ods for LDSs such as subspace identification (van Overschee and de Moor. 1996, Qin, 2006), or improper learning (Anava et all, 2013; Hazan all 2017: Kozdoba et all, 2019. do not apply to the system (m),,,, due to non-stationarity. On the other hand, the methods that do apply to (II) either lack guarantees, or have only asymptotic analysis which in addition relies strongly on Guassianity of the noises
2

---

Moreover, our approach differs significantly from the existing methods. We show that the structure of equations (M) is closely related to, and inherits several important properties from, the classical discrete Laplacian operator on the line leading to new arguments that have not been explored in the literature. In particular, we use this connection to show that an appropriate inversion of the system produce estimators that are concentrated enough so that and may be recovered. The heart of the paper is the new definition of the estimators that exploits explicitly the shape of a certain data dependent operator, and the subsequent concentration analysis. In particular, this approach yields the first known geometric parameters of the data that can be used to bound convergence rates.
The rest of the paper is organized as follows: G contains the necessary definitions. In Section we describe in general lines the methods and the main results of this paper. The technical estimates on certain operator spectra, that are critical to the analysis and may be of independent interest, are stated in Section H Section we present experimental results on synthetic and real data. Due to space constraints, while we outline the main arguments in the text, the full proofs are deferred to the Supplementary Material.
The related work is discussed in Section D and Section
2 Literature
We refer to Chu and Chen (201);; Hamilton (1994); Anderson and Moore (11979); Shumway and Stotter (20TI) for a general background on LDSs, the Kalman Filter and maximum likelihood esti- mation.
Existing approaches to the variance estimation problem may be divided into three categories: (i)
General methods for parameter identification in L.S,, either via maximum likelihood estimation (MLE) (Hamilion, (199)), via subspace identification (ban Overschee and de Moor, [196; Qun.
2006) In particular, finite sample bounds for system identification were given in (Campi and Weyer
2003 Vidyasagar and Karandikar, 2006) and in the recent work siamis and Pappas (2019). (ii)
Methods designed specifically to learn the noise parameters of the system, developed primarily in the control theory community, in particular via the innovation correlattnl function, such the
classical Mchra (1970); Belanger (1974), or for instance more recent Wang et all (2017); Dunik al (2018). (iii) Improper Learning methods, such as Anava all (2013); [Hazan et all (2017); Kozdobal al (2019). In these approaches, one does not learn the LDS directly, but instead learns a model from certain auxiliary class and shows that this auxillary model produces forecasts that are as good as the forecasts of an LDS with "optimal" parameters.
Despite the apparent simplicity of the system (M)-2),, most of the above methods do not apply to this system. This is due to the fact that most of the methods are designed for time invariant, asymp- totically stationary systems, where the observation operator (u, in our notation) is constant and the Kalman gain (or, equivalently Cru in eq. (G)) converges with t. In particular this limitation exists in all the system identification results cited above, and is essential to the approaches taken there. However, if the observation vector sequence ut changes with time a necessary property for the dynamic regression problem the system will no longer be asymptotically stationary. In particular, due to this reason, neither the subspace identification methods, nor any of the improper learning approaches above apply system (I)((()
Among the methods that do apply to (- are the general MLE estimation, and some of the auto- correlation methods (Belanger, 1974 Dunik et all 2018). On one hand, both types of approaches may be applicable to systems apriori more general than (()) On the other hand, the situation with consistency guarantees- the guarantee that one recovers true parameters given enough observations for these methods is somewhat complicated. Due to the non-convexity of the likelihood function, the MLE method is not guaranteed to find the true maximum, and result the whole method has no guarantees. The results in Belanger (19974;; Dunikeral (2018) do have asymptotic consistency guarantees. However, these rely some explicit and implicit assumptions about the system, the sequence u, in our case, which can not be easily verified. In particular, [eelaner (19724); Dunn call (2018) assume uniform observability of the system, which we do not assume, and in addition rely on certain implicit assumption about invertibility and condition number of the matrices related to the sequence Ue. Moreover, even if one assumes that the assumptions hold, the results are purely asymptotic, and for any finite T, do not provide a bound of the expected estimation error as a function of T and (ur)/=1-
3

---

In addition, mentioned earlier, MLE methods by definition must assume that the noises are Gaus- sian (or belong to some other predetermined parametric family) and autocorrelation based methods also strongly use the Gaussianity assumption. Our approach, on the other hand, requires only sub Gaussian noises with independent coordinates. We note that there are straightforward extensions of our methods to certain cases with dependencies. Indeed, the operator analysis part of this paper does not depend on the distribution the noises. Therefore, to achieve such an extension, one would only need to correspondingly extend the main probabilistic tool, the Manson-Wright inequality (Hanson and Wright [197,, Rudelson et all 2013, see also Section M and Supplementary Material Section (E). One such extension, for vectors with the convex concentration property, was recently obtained in AdamczaR (20115).
3
Notation
We refer to Bhatra (1997)) and Vershynin (2018) general references on the notation introduced below, for operators and sub Gaussian variables, respectively.
Let A R" - R be an operator with singular value decomposition A U .... Ag) J,),W, where 8 S min {m, n and x 2 X2 2 2 A > 0. Note that singular values are strictly positive by definition (that is, vectors corresponding to the kernel of A do not participate in the decomposition == U Di...)A....,,... as ||A||
W). The Hilbert-Schmidt (Frobenius) norm is defined
HS = Ei= XZ. The nuclear and the operator norms are given by ||Al|nue Σi= A and
|| = \1 respectively.
A centered (EX 0) scalar random variable X is sub-Gaussian with constant K, denoted X ~ SG(K), if for all 1 > O satisfies P (|X) > t) I 2exp xz
A random vector X = (X ...... Xm) is A. sub-Gaussian, denoted X ~ SG (1s). for every v E R. with |o| = 1 the random variable (v. X) is x sub-Gaussian. A random vector X is a-isotropic if for every E R with |0| 1, IE X) = o².
Finally, a random vector X = (X Xm) ) is o-isotropically x sub-Gaussian with independent components, denoted X ISGm if X, are independent, and for all S m, EX,, EX² 0² and X, ~ SG(R). Clearly, if X ~ (SGm(a, k)) then X is o-isotropic. Recall also that X ISG (o, k) implies X SG, (x) (Vershynin 2018). The noise variables we discuss in this paper are ISG(K,
Throughout the paper, absolute constants are denoted by e, G, c", ,..... etc. Their values may change from line to line.
4 Overview of the approach
We begin by rewriting (M) vector form. To this end, we first encode sequences of T vectors in " {athir C RR, as a vector E RTT, constructed by concatenation of 0z 's. Next, we define the summation operator S' RT FT which acts on any vector (h1, hh,....,hu hr) IRT by
S' ((11,2 he...., hr) (h1, h1 h2
Σn, h4 N,.DM hd).
(6)
SST-1
EST
Note that S' is invertible operator. Next, we similarly define the summation operator S RR IT', an n-dimensional extension of S'. which sums n--imenniiona vectors. Formally, for (hx) e RTT, and for 1 S j S n, < < T, (Sh) sequence of process noise terms fa,...., .a....,hhy hr IR" is viewed as vector e RT",, then by definition encoding of the sequence X₁- Sh is the I
Eige m) Observe that the
In words, coordinate of O,, Next, given a sequence of observation vectors uI,..., Me E R", we define the observation operator
O, : TT T by (O,æ), (u,,
the inner product between u, and t-th part of the vector T" Define also (Yi,..., Yr) e RT to be the concatenation of Y,..., Yr. With this notation, one may equivalently rewrite the system (M)( as follows:
Y O,,Sh +

---

where h and are independent zero-mean random vectors in IR"" and R" respectively, with inde- pendent sub Gaussian coordinates. The variance of each coordinate of is r² and each coordinate of has variance rf².
Up now, have reformulated our data model as single vector equation. Note that in that equation, the observations Y and both operators On and S are known to us. Our problem may now be reformulated as follows: Given Y R², assuming Y was generated by ((A), provide estimates of o, 1)-
As motivation, we first consider taking the expectation of the norm squared of eq- (τ). For any operator A : Rm + Rm and zero-mean vector h with independent coordinates and coordinate variance o², we have | = || s o², where [|A]| is the Hilbert-Schmidt (or Frobenius) norm of A. Taking the norm and expectation of (.), and dividing by T we thus obtain
E E
10.5 HS T
0 T Tz"
T
Next, note that |O,S||As is known, and an elementary computation shows that ||0., stant order (as function of T; see (25)), while the coefficient of ry² is 1 Thus, if the quantity lY] T. were close enough its expectation with high probability, we could take this quantity as a (slightly
is of con-
T
biased) estimator of o². However, as it will become apparent later, the deviations of [Y]2
" around the
expectation are also of constant order, and thus |X| T can not be used as estimator. The reason for
these high deviations of |Y] T is that the spectrum of O,,S' is extremely peaked. The highest squared singular value of S, is of order 72, the same order as sum of all of them, 1 Contrast this with the case of identity operator, Id : IR - RTT, We have E [ = E [h|² = Tno? and one can also easily show that, for instance, V Var Id(h)) Tno², and thus the deviations of order /Tnio- smaller order than E [dd(h)) While for the identity operator the computation is elementary, for general operator A the situation is significantly more involved, and the bounds on the deviations of y will be obtained from the Hanson- Wright inequality (Hanson and Wright I971, see also Rudelson all (2013)), combined with standard norm deviation bounds for isotropic sub Gaussian vectors.
With these observations in mind, we proceed to flatten the spectrum of O., by taking the pseudo- inverse. Let R : RT + TT be the pseudo-inverse, or Moore-Penrose inverse of O,S. Specifically, let
O,S = U a Diag()...... be the singular value decomposition of O,S, where n 2 2 PT are the singular values.
W,
(9)
For the rest of the paper, we will assume that all of the observation vectors uN are non-zero. This assumption is made solely for notational simplicity and may easily be avoided, as discussed later in this section. Under this assumption, since invertible and O₂ has rank T, we have X > 0 for all S T. For T, denote Xi - Then Xi are the singular values of R, arranged in non-increasing order, and we have by definition
R V' o Diag(Xr, -or-OXXX--------... Xr-1,..., X2, X2,X1) OU",
(10)
where W U denote the transposed matrices of U, V, defined in (5))
Similarly to Eq. (,, we apply Rt to (7), and since || Hs = T, by taking the expectation of the squared norm we obtain
|RY T
|
E EIRY12)
=
(11)
T
In this equation, the deviation term
is of order O ( with high probability (Theorem M. Moreover, the coefficient of o² is 1, and the coefficient of n72 which is IR is of order at least 2 7) (Thcorem EL see Section bl for additional details) much larger order than Since 2 and |[]]]
[RY|2 BIRVI)
i are known, it follows that we have obtained one equation satisfied by 2 and 1 up to an error of where both coefficients are of order larger than the error.
5

---

Algorithm Spectrum Thresholding Variance Estimator (STVE)
1: Input: Observations Y, observation vectors ut, with < T, and aT.
2: Compute the SVD of O,S,
OnS = Uo Diag(?,...... , W,
where Vi 2 72 2 ...
YT > Denote X TT+- for SisT.
3: Construct the operators
R= W Dia.
and
R' = W
x1) U
4: Produce the estimates:
|RY|
IRY 22)
|
[ iss
T
T
|
|R] sp.
0
T
Next, we would like to obtain another linear relation between 0², 1/2. To this end, choose some p p= oT' where 0 < a < is of constant order. The possible choices of p are discussed later this section. We define operator R' RR + R to be version of R truncated to the first singular values. (IT)) is the SVD decomposition R, then
R' = W* Diag(O, O.... Xpr xx))
Similarly to the case for R, we have
R'Y o²
IR'SS |R
HS 7²
|RY
E R
(12)
p
p
The deviations in (122) are also described by Theorem M. Note also that since
| | HS
is the sum of | T
p largest squared singular values R, by definition it follows that || 1
Now, given two equations in two unknowns, we can solve the system to obtain the estimates 22 and - The full procedure is summarized in Algorithm M. and the bounds implied by Theorem M on the estimators x² and n7 are given in Corollary Z. We first state these results, and then discuss in detail the various parameters appearing in the bounds.
Theorem 1. Consider random vector Y TR. the form Y O2,Sh+ where h ISGTa(O, and 2 ~ (S, 1s). Set [umin = min, [ue]. Then for any 0 < 5 < 1,
RY 2
|
B
2 +
2
C
I 48
(13)
T
T
VT
Y
B Vp
P
0² +
HS 1/2
2 c
S 48
(14)
P
where B is given by
B (1 + i22) (1 [yman -2)
1 log
(15)
8
The bounds on the estimators of Algorithm are given in the following Corollary. As discussed below, in addition to Theorem MI the key to the derivation of this Corollary are the estimates of the spectrum of R, given Theorem B, Section B

---

I' spectrum xx"
2.0
p
T/4
R' spectrum x ( his/ |
( hslp W
T T/4
24
1.0
0.5
160
300
(a) Spectrum (x)) of R for electricity data (blue). (p) / as function of p (orange). The mean RR] Hs /T (green).
(b) Similar figure for synthetic data, T 500.
Figure 1: Spectra of R
Corollary 2. Lel 02, 2p be the estimators of o² ,²,12² obtained form Algorithm 0 with p 2 I.. Set [umaz| max, |ue|. Then for any O < 8 < 1, with probability at least 86,
B
p |RI
2
S
(16)
VT
T |
B
p |R] HS
2
2
<
n² Umaa Plog2 T,
(17)
VT
T || HS
with B given by (IS))
We first discuss the assumption [Umin| 0. This assumption is made solely for notational con- venience, as detailed below. To begin, note that some form of lower bound the norms of the observation vectors Ue must appear in the bounds. This is simply because if one had u, = 0 for all T, then clearly no estimate of would have been possible. On the other hand, smallest value |umin| may seem restrictive at first. We note however, that instead of considering the observation operator O,, : TT + R one may consider the operator On T + IR for any subsequence ( The observation vector would be correspondingly restricted the subse- quence of indices. This allows us to treat missing values and to exclude any outlier u, with small norms. All the arguments
of the
Theorems M and & hold for this modified O₂ without change. The only price that will be paid is that 7 will be replaced by 7 in the bounds. Moreover, we note that typically we have |ul] 2 1 by construction, see for instance Section 6.2 Additional discussion of missing values may be found in Supplementary Material Section H
Next, up to this point, we have obtained two equations, that in order to be able to obtain 1/2 from these equations, at least one of the coefficients of r7²,
(12), in two unknowns. 0² o7. Note
either | or || must be larger order than L. the order deviations. Providing lower vT' bounds on these quantities is one of the main technical contributions of this work. This analysis uses the connection between the operator S and the Laplacian on the line, and resolves the issue of translating spectrum estimates for the Laplacian into the spectral estimates for R. We note that there are no standard tools study the spectrum of R, and our approach proceeds indirectly via the analysis of the nuclear norm of O, These results are stated in Theorem B. In particular, we show that || is ( 1og2 which the source of the log factor in (17).
T
Finally, order to solve the equations (ITI) (12), not only the equations must have large enough
coefficients, but the equations must be different. This is reflected by the term (1 z R|
in
TI 1
| /v 2 1 by definition, we would like to have
Equivalently, while
|R Hs /p
2 + const
(18)
HS
7

---

for the bounds (L))) (TZ)) to be stable. Note that since both | and |Rliss are computed in Algorithm M. the condition (IK) can simply be verified before the estimators o², 172 are returned.
It is worth emphasizing that for simple choices of say p T, the condition (I)) does hold in
| Rs
practice. Note that, for any p. we can have
= only if the spectrum of R is constant.
[| it
Thus (T)) amounts to stating that the spectrum of R exhibits some decay. As we show in experiments below, the spectrum (squared) of R. for derived from daily temperature features, for random Gaussian ut, indeed decays. See Section & Figures Ia and IT. In particular, in both cases (TS) holds
with const > 1. Additional bounds on the quantity (1 plRRi under various assumptions on the sequence u, are given in Section G the Supplementary Material.
TIIR" lis
Properties of S,S and
As discussed in Section A (see the discussion following eq. (II),, one of the crucial points enabling Algorithm and its analysis is the fact that the quantity |R]e is bounded below by an expression that is of much higher order than the noise magnitude VT'
In this section we provide the formal statement of this and other associated results, and discuss the related arguments. First, we obtain the following bound on the spectrum of O (Lemma Supplementary Material Section E). Recall that the nuclear norm was defined in Section B, and that
for sequence u, we set [Mmax max; |ug| and |umin = min, |ue]. Then:
|10,51 Σ A) I 4n Umaa T log
(19)
The proof of this bound exploits the connection between S and the Laplacian on the line. In particular, we use the fact that the eigenvalues the Laplacian are known precisely, satisfying Next, state the lower (and upper) bounds for R.
λ sin (TT-D)
Theorem 3. Let R R + RT be the pseudoinverse of S. Then
S || S 2 |mmn
(20)
n Uman log T
T
|| HS S Ummii|
T.
(21)
n Umax
2 log" 2
T
1
|| R S [Mmin]
(22)
n |mar
log 4 T
Due to the complicated structure of R as a pseudo-inverse of a composition of operators, there are no direct ways to control individual eigenvalues of R. Thus the main technical issue resolved in Theorem & is nevertheless obtaining lower bounds || iss Our approach is rather indirect, and we obtain these bounds from the nuclear norm bound (LL) via a Markov type inequality on the eigenvalues.
6 Experiments
6.1 Synthetic Data
In this section the performance of STVE is evaluated synthetic data. The data generated by the LDS (M),,, using Gaussian noises with o² = 0.5, n7² = 2. The input dimension was n = and the input sequence 14 sampled from the Gaussian T.)
We run the STVE algorithm for different values of T, where for each we sampled the data 150 times. Figure Ed shows the average (over 1 150 runs) estimation error for both process and observation noise variances for various values observed in Figure de that the estimation emrors decay roughly the rate of const/ /T. typical spectrum of R. is shown in Figure [R. For larger T, the spectra also exhibits similar decay.
As expected from the bounds in Corollary E may be

---

Fdll Dea Pog Secoel
Tarkeay Prgresiin
(a) STVE variance estimation errors T, and the T decay.
(b) Load (y-axis) against Tempera- (c) Smoothed prediction errors. Sta- ture (x-axis), both axis normalized. tionary Regression trained on first The full data (blue points), regres- half (blue), MLE (orange), STVE sion learned on the first half (or- (green), OG (red).
ange), regression learned the sec- ond half (green).
Figure 2: Evaluation
6.2 Temperatures and Electricity Consumption
In this section we examine the relation between daily temperatures and electricity consumption in the data from Hong all (201E) (see also Hong (2016)). The following forecasting methods are compared: a stationary regression, an online gradient, and Kalman filter for dynamic regression, with parameters learned via MLE STVE. We find that the Kalman filter methods provide the best performance, with significant difference between STVE and MLE derived systems.
The data consists of total daily electricity consumption (load) yt, and the average daily temperature, ut, for a certain region, for the period Jan-2004 to Jun-2008, Full details on the preprocessing of the data, as well as additional details on the experiments, are given in Supplementary Material Section II Here we note that the data contains missing load values, for 9 non-consecutive weeks (out of about 234 weeks total). All methods discussed here, including STVE, can naturally incorporate missing values, as discussed in Supplementary Material Section F.
An elementary inspection of the data reveals that the load may be reasonably approximated quadratic function of the temperature, = It,1 + F2,2 Ut + Ft,3 - oz, where 4, = (1, or, vz) is the observation vector (features), and 2t (are, 24,2, vector. This is shown in Figure RE, where we a stationary (time invariant) regression of the above form, using either only the first only the second half of the data. We note that these regressions differ the regression vector changes with time. It is therefore of interest to track it via online regression.
the possibly time varying regression
We use the first half of the data (train set) to learn the parameters a, 17 of the online regression (II) via MLE optimization and using ST We also the train set to find the optimal learning rate a for the OG forecaster described by the update equation (5)). This learning rate is chosen as the rate that yields smallest least squares forecast error the train In addition, learn time independent, stationary regression on the first half of the data.
We then employ the learned parameters to make predictions of the load given the temperature, by all four methods. The predictions for the system (W)-L) are made with a Kalman filter (at time t, we use the filtered state estimate ž., which depends only on y1...., Vr and uy,..., Ue, and make the prediction I (ā,, 1t)))).
Daily squared prediction errors (that is, (yr - T.) are shown in Figure (smoothed with a moving average of 50 days). We see that the adaptive models (MLE, STVE. OG) outperform the stationary regression already on the train set (first half of the data), and that the difference in performance becomes dramatic on the second half (test). It is also interesting to note that the performance of the Kalman filter based methods (MLE, STVE) is practically identical, but both are somewhat better than the simpler OG approach.
We also note that by construction, we have [Wmin| 2 this experiment, due to the constant coordinate, and also |Umaa] < 5, due to the normalization. Since these operations are typical for any regression problem, we conclude that the direct influence of [mmin| and |umarl on the bounds in Theorem M and Corollary Z will not usually be significant.

---

7 Conclusion And Future Work
In this work we introduced the STVE algorithm for estimating the variance parameters of LDSs of type (I),)), and obtained the first sample complexity guarantees for such estimators. We have also shown how the shape of the spectrum of can be exploited to obtain the estimators and the related bounds, thus providing the first explicit geometric parameter of the data that affects the bounds.
As discussed in Section and demonstrated in Section & the system (I) of independent interest in upplications. However, we also believe that the analysis presented here is an important first step towards a finite time data-dependent quantitative understanding of general LDSs, and perhaps even non-linear dynamical systems.
Acknowledgments and Disclosure of Funding
This research was supported by the ISRAEL SCIENCE FOUNDATION (grant No. 2199/20).
References
Adamczak, R. (2015). note on the wright inequality for random vectors with dependencies. Electronic Communications Probability, 20.
Anava, O., Hazan, E., Mannor, S., and Shamir, O. (2013). Online learning for time series prediction. In COLT 2013 The 26th Annual Conference Learning Theory, June 12-14, 2013, Princeton University, NJ, USA.
Anderson, B. and Moore, J. (1979). Optimal Filtering. Prentice Hall.
Belanger, R. (1974). Estimation of noise covariance matrices for linear time-varying stochastic process. Automatica, 10(3).
Bhatia, R. (1997). Matrix Analysis. Graduate Texts in Mathematics. Springer New York.
Campi. M. and Weyer, E. (2005). Guaranteed non-asymptotic confidence regions in system identification. Automatica, 41(10):151-1174.
Chui, C. and Chen, G. (2017). Kalman Filtering: with Real-Time Applications. Springer Interna- tional Publishing.
Dunik, J., Kost, O., and Straka, O. (2018). Design of measurement difference autocovariance method for estimation of process and measurement noise covariances. Automatica, 90.
Gohberg, I. and Krein, M. (1969). Introduction to the Theory of Linear Nonselfadjoint Operators. Translations mathematical monographs, American Mathematical Society.
Hamilton, J. (1994). Time Series Analysis. Princeton University Press.
Hanson, D. L. and Wright, E. T. (1971). A bound on tail probabilities for quadratic forms in inde- pendent random variables. Ann. Math. Statist., 42.
Hazan, E. (2016). Introduction to online convex optimization. Found. Trends Optim.
Hazan, E., Singh, K., and Zhang, C. (2017). Online learning of linear dynamical systems. In Advances in Neural Information Processing Systems, pages 6686-6696.
Hong. T. (2016).
Tao hongs's blog.
http://blog. drhongtao. comp:ttoo.ccm00000000 comm2016/07 data. data. Accessed: 1/8/2019.
get www.o2002--ddd-oeeaaaaaaaaa
Hong, T., Pinson, P., and Fan, S. (2014). Global energy forecasting competition 2012. International Journal of Forecasting, 30:357-363.
Kozdoba, M., Marecek, J., Tehrakian, T., and Mannor, S. (2019). On-line learning of linear
dynamical systems: Exponential forgetting in Kalman filters. AAAI.
10

---

Mehra, R. (1970). On the identification of variances and adaptive kalman filtering. IEEE Transac- tions on Automatic Control, 15(2).
Mitchell, A. and Griffiths, D. (1980), The finite difference method in partial differential equations. Inter--nnerrccience publication. pubbicattionn.. Wiley.
Petris, G. (2010). An package for dynamic linear models. Journal of Statistical Software.
Qin, J. (2006). An overview of subspace identification. Computers chemical engineering, 30(10-12):100--5513
Rudelson, M., Vershynin, R., et al. (2013). Hanson-wright inequality and sub sub-gaussian concentra- tion. Electronic Communications Probability, 18.
Shumway, R. and Stoffer, D. (2011). Time Series Analysis and Its Applications (3rd ed.).
Tsiamis, A. and Pappas, G. J. (2019). Finite sample analysis of stochastic system identification. In 2019 IEEE 58th Conference on Decision and Control (CDC), pages 3648-3654. IEEE.
van Overschee, P. and de Moor, L.. (1996). Subspace identification for linear systems: theory, imple- mentation, applications. Kluwer Academic Publishers.
Vershynin, R. (2018). High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
Vidyasagar, M. Karandikar, R. L. (2006). learning theory approach to system identification and stochastic adaptive control. Probabilistic and randomized methods for design under uncer- tainty, pages 265-302.
Wang, H., Deng. Z., Feng, B., Ma, H., and Xia, Y. (2017). An adaptive kalman filter estimating process noise covariance. Neurocomputing, 223.
West, M. and Harrison, (1997). Bayesian Forecasting and Dynamic Models (2nd ed.). Springer- Verlag.
Zinkevich, M. (2003). Online convex programming and generalized infinitesimal gradient ascent. ICML.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [N/A] The
main focus of this work is theoretical analysis.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] Thus Supplementary Materail contains all the proofs. Proofs outlines given in the main body of the paper.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main exper- imental results (either in the supplemental material or as URL)? [Yes] The data is publically available, see the references in Section The short code is not provided the moment, but can be fully derived from Algorithm MI

---

(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Sections f and I.
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes]
4. If you are using existing assets (c.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [N/A] The link to the data author's docu- mentation is provided.
(c) Did you include any assets either in the supplemental material or URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifi- able information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals. if applicable? [N/A]
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A]
12


Reviews:
[{"review_id": "nQ1L-_Kzx6C", "reviewer": "Reviewer_p7x7", "summary": "This paper studies estimation of process and observation noise variance in a subclass of linear dynamical systems whose states evolve according to $X_{t+1} = X_t + h_t$ and observations follow $Y_{t} = \\langle X_t, u_t \\rangle + z_t$, in which $u_t$ are known observation vectors. The paper presents the first algorithm for the estimation of variance in this class of dynamical systems that enjoys finite-sample guarantees. The approach revolves around forming variance estimators that exploit data-dependent operators. ", "questions": "Additional related work suggestions that seem relevant:\n- Rashidiejad, Paria, Jiantao Jiao, and Stuart Russell. \"SLIP: Learning to Predict in Unknown Dynamical Systems with Long-Term Memory.\" 34th Conference on Neural Information Processing Systems (NeurIPS 2020), Vancouver, Canada. 2020.\n- Tsiamis, Anastasios, and George J. Pappas. \"Linear systems can be hard to learn.\" In 2021 60th IEEE Conference on Decision and Control (CDC), pp. 2903-2910. IEEE, 2021.", "limitations": "Limitations are discussed. The work is mainly theoretical.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths** \n- The paper presents the first variance estimation for the LDS subclass with finite-sample guarantees. The resulting algorithm is simple.\n- As far as I know, techniques used in the derivation of the approach are original.\n- The paper is reasonably well-written. The authors provide intuition regarding their choices and approach and compare them with alternatives.\n\n**Weaknesses**\n- The system considered in this work is a subclass of LDS with scalar observations, which can limit the applicability of the approach.\n- The experiment section is a bit weak (although the paper is mainly theoretical). At least based on the experiment in Section 6.2, STVE does not perform better than MLE (followed by the Kalman filter). Additional empirical evaluation that highlights the benefits of STVE (e.g. in handling non-Gaussian noise) would strengthen the paper). ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "MQI1AQ1PsxB", "reviewer": "Reviewer_1GcG", "summary": "This paper provides novel finite-sample guarantees for parameter estimation in the dynamic linear regression setting, where the true linear model’s coefficients vary across time. By recharacterizing the dynamic regression problem into a linear dynamical system, the authors note that existing methods, such as the Kalman Filter, rely on knowing certain variance parameters. While one could estimate these variance parameters using maximum likelihood estimation, the accompanying Gaussian noise assumption is perhaps too strong and there are only asymptotic guarantees of the validity of the estimates. The authors of this paper have developed a procedure with finite sample guarantees for estimating the variance parameters with a less onerous sub-Gaussian noise assumption for the dynamical system. With a simulation, the authors demonstrate that the novel estimator’s error decays at a rate in line with the theoretical guarantees. In an empirical example of predicting electricity consumption over time, the authors show that their estimator, once inputted into a Kalman Filter algorithm, performs comparably to using maximum likelihood estimations as the inputs. ", "questions": "The field of data assimilation often deals with methods like the Kalman filter with unknown parameters; I was surprised not to see this mentioned and it would be helpful for the authors to explicitly address their contributions relative to work in that field. For instance, how does the submission compare with https://arxiv.org/abs/1903.09122, https://arxiv.org/abs/1912.12309, or https://arxiv.org/abs/2001.06270?", "limitations": "The paper discusses limitations naturally for the most part, given the extensive discussion of comparable methods and highlighting the new assumptions the novel approach relies on. The authors do note, briefly, at the end of the Experiments section that their novel approach performs “practically identical” to the maximum-likelihood estimation approach. While this is not at odds with the main claim of the authors of finally providing an estimator with finite sample guarantees, they are essentially acknowledging a limitation of the new estimator’s performance versus existing methods.\n\nAs far as the discussion on societal impact and the like, the paper eschews such discussion due to the theoretical nature of research.\n", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n\nThe paper provides novel finite-sample guarantees in a useful generalization of linear regression. The proofs seem innovative and non-trivial, providing bounds that may have previously seemed out-of-reach. In particular, the use of existing theory on the discrete Laplacian operator on the line seems like an especially novel approach that other theoreticians will find intriguing and useful.\n\nThe authors do a great job of showing how the problem set-up, which starts very general, can be recharacterized into a more specific problem that has been well-studied. No less, the authors home in on an aspect that has previously been untouched, providing bounds on estimating variance parameters in the dynamical system.\n\nThe clarity of the paper’s writing and the cogent overview of the literature were also appreciated. Via the problem set-up section, the authors provided very clear exposition on the links to previous literature and what the novel contributions of this paper is, aided by uniting the notation of their approach with others. The writing was clear throughout, and the authors provided ample details and justification for the proof sketches in the body of the paper.\n\nWeaknesses\n\nThe two main ways to improve the paper are in enriching the experimental section and making more clear the importance of the novel finite sample guarantees.\n\nIn the empirical sections, it seems like MLE often does not differ much from the novel approach. While the authors briefly note this, it would be useful to have them discuss the implications of this more. For example, under what circumstances do they differ, and can we learn anything from those settings? In addition, for the performance comparisons on the real world data example, it would be helpful to list MSE or some other test metrics. While plotting the predictions versus the truth was certainly a step in the right direction, it could seem like obfuscation to not provide the metrics, especially since it seems like the novel estimator may have performed worse than the MLE. For the simulation study, the sensitivity analysis provided on the number of samples was certainly of utmost importance. However, the authors should also note why they chose some of the settings (e.g., why the number of covariates was five or the choices for the ground truth variance parameters). It was not clear whether these settings were chosen due to computational constraints, to enhance the signal-to-noise ratio, to adhere to some realistic empirical examples, or if they were just drawn out of a hat.\n\nIt may be useful to provide more color on the usefulness of better variance parameter estimates and these guarantees. Are the parameters only useful insofar as inputs to the Kalman Filter and related approaches to get better dynamic coefficient estimates? Or are these variance parameters intrinsically interesting?\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "H7ZaSfvvwa5", "reviewer": "Reviewer_Wf8E", "summary": "The paper considers considers non-stationary linear regression with sub-Gaussian observation noise z (of variance \\eta^2), where the regression vector X undergoes at each time t an additive modification by a process-noise vector h, also sub-Gaussian and zero mean with independent components of variance \\sigma^2. \nEstimation of the true trajectory and forecasting can be done by Kalman-filtering, once the true noise-parameters are known, and the principal contribution of the paper is to provide and analyze an algorithm to estimate the variances \\sigma^2 and  \\eta^2. The approach is an alternative to maximum likelihood estimation (MLE).\n\nA compressed form of the process equations leads to a linear equation in  \\sigma^2 and  \\eta^2 modulo a fluctuation term. Applications of the Hanson-Wright inequalitiy bound these fluctuations and show that they are of smaller order than the coefficients in the equation. A second equally structured equation is obtained from a truncated time-series, and the two equations are solved for \\sigma^2 and  \\eta^2. The paper is careful to demonstrate that this system of equations is in general well conditioned relative to the magnitude of fluctuations.\n\nAn experimental section compares the algorithm to competing methods and works about well as MLE. ", "questions": "Some typos:\n\nl 146, square missing in expectation\n\nl 185 and l 186, \"the\" missing before \"sum\" and \"identity\"\n\nl 199, I believe it should be \\gamma instead of \\lambda\n\np 19, equation (78). u_max and u_min seem exchanged.", "limitations": "There is no section on limitations, but I am also not clear on what it should contain for this paper.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n1. In contrast to MLE the approach has performance guarantees for the more general sub-Gaussian noise model.\n2. The paper contains a number of interesting mathematical tricks.\n\nWeaknesses:\n1. The presentation is a bit opaque. Section 4 could perhaps be a lot more readable if portioned into sub-sections: Obtaining the representation (7) - why using the pseudo-inverse and obtaining (11) - obtaining the second equation (12) - etc.\n2. The role of Theorem 3 and Section G could be better explained. The most problematic part of the error bounds in Corollary 2 (inverse parenthesis) can already be observed when executing the algorithm, so from a practical point of view Theorem 3 and Section G are superfluous. Their interest lies in showing that generically the division in step 4 of the algorithm is unproblematic.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept