Finite Sample Analysis Of Dynamic Regression Parameter Learning

TL;DR

This paper introduces finite sample guarantees for variance estimation in dynamic linear regression, analyzing the system operator spectrum to derive explicit convergence rates applicable to arbitrary sub-Gaussian noise.

Contribution

It provides the first finite sample complexity bounds for variance estimation in dynamic regression, linking convergence rates to spectral properties of the system operator.

Findings

Derived spectrum estimates of the system operator.

Established finite sample guarantees for variance estimators.

Validated approach on synthetic and real-world data.

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.

Figures

Equations

$$X_{t+1}$$

$$Y_t$$

$$C_{t+1}$$

$${\bar{x}}_{t+1}$$

$${\bar{x}}_{t+1}={\bar{x}}_t+\alpha \cdot u_{t+1}\cdot \left(Y_{t+1}-⟨{\bar{x}}_t,u_{t+1}⟩\right),$$

$$S^{\prime }(h_1,h_2,\dots ,h_T)=(h_1,h_1+h_2,\dots ,\sum _{i\le T-1}{h}_i,\sum _{i\le T}{h}_i).$$

$$Y=O_{u}Sh+z,$$

$$\frac{\mathbb{E}{\left|Y\right|}^2}{T^2}=\frac{{\Vert O_{u}S\Vert }_{HS}^2}{T^2}{\sigma }^2+\frac{T}{T^2}{\eta }^2.$$

$$O_{u}S=U\circ Diag({\gamma }_1,\dots ,{\gamma }_T)\circ W,$$

$$R=W^{\ast }\circ Diag({\chi }_T,{\chi }_{T-1},\dots ,{\chi }_2,{\chi }_1)\circ U^{\ast },$$

$$\frac{{\left|RY\right|}^2}{T}={\sigma }^2+\frac{{\Vert R\Vert }_{HS}^2}{T}{\eta }^2+\left(\frac{{\left|RY\right|}^2}{T}-\frac{\mathbb{E}{\left|RY\right|}^2}{T}\right).$$

$$R^{\prime }=W^{\ast }\circ Diag(0,0,\dots ,{\chi }_p,{\chi }_{p-1},\dots ,{\chi }_1)\circ U^{\ast }.$$

$$\frac{{\left|R^{\prime }Y\right|}^2}{p}={\sigma }^2+\frac{{\Vert R^{\prime }\Vert }_{HS}^2}{p}{\eta }^2+\left(\frac{{\left|R^{\prime }Y\right|}^2}{p}-\frac{\mathbb{E}{\left|R^{\prime }Y\right|}^2}{p}\right).$$

$$O_{u}S=U\circ Diag({\gamma }_1,\dots ,{\gamma }_T)\circ W,$$

$$R=W^{\ast }\circ Diag({\chi }_T,\dots ,{\chi }_1)\circ U^{\ast }$$

$$R^{\prime }=W^{\ast }\circ Diag(0,0,\dots ,0,{\chi }_p,\dots ,{\chi }_1)\circ U^{\ast }$$

$$\widehat{{\eta }^2}$$

$$\widehat{{\sigma }^2}$$

$$\mathbb{P}\left(\left|\frac{{\left|RY\right|}^2}{T}-\left({\sigma }^2+\frac{{\Vert R\Vert }_{HS}^2}{T}{\eta }^2\right)\right|\ge c\frac{B}{\sqrt{T}}\right)\le 4\delta ,$$

$$\mathbb{P}\left(\left|\frac{{\left|R^{\prime }Y\right|}^2}{p}-\left({\sigma }^2+\frac{{\Vert R^{\prime }\Vert }_{HS}^2}{p}{\eta }^2\right)\right|\ge c\frac{B}{\sqrt{p}}\right)\le 4\delta$$

$$B=\left(1+{\kappa }^2\right)\left(1+{\left|u_{min}\right|}^{-2}\right)\log \frac{1}{\delta }.$$

$$\left|\widehat{{\sigma }^2}-{\sigma }^2\right|$$

$$\left|\widehat{{\eta }^2}-{\eta }^2\right|$$

$$\frac{{\Vert R^{\prime }\Vert }_{HS}^2/p}{{\Vert R\Vert }_{HS}^2/T}\ge 1+const$$

$${\Vert O_{u}S\Vert }_{nuc}=\sum _{t\le T}{\lambda }_t(O_{u}S)\le 4n\left|u_{max}\right|T\log T.$$

$$c\frac{1}{n\left|u_{max}\right|\log T}$$

$$c\frac{1}{n^2{\left|u_{max}\right|}^2{\log }^{2}T}T$$

$$c\frac{1}{n^4{\left|u_{max}\right|}^4{\log }^{4}T}T$$

$${\lambda }_1(O_{u}S)\le \left|u_{max}\right|\cdot T\text{mbox}and{\lambda }_T(O_{u}S)\ge \frac{1}{2}\left|u_{min}\right|$$

$${\Vert O_{u}S\Vert }_{nuc}=\sum _{t\le T}{\lambda }_t(O_{u}S)\le 4n\left|u_{max}\right|T\log T$$

$$\frac{1}{4}{\left|u_{min}\right|}^2{T}^2\le {\Vert O_{u}S\Vert }_{HS}^2=\sum _{t\le T}t{\left|u_t\right|}^2\le \frac{1}{2}{\left|u_{max}\right|}^2{T}^2.$$

$$2\sin \left(\frac{\pi (T-t)}{2(T+1)}\right)\le {\lambda }_t(S_T^{\prime -1})\le 2\sin \left(\frac{\pi (T+1-t)}{2(T+1)}\right)$$

$$\frac{1}{T}\le {\lambda }_T(S_T^{\prime -1})\le 2\sin \left(\frac{\pi }{2(T+1)}\right).$$

$${\Vert S\Vert }_{op}\le T\text{ and }\left|S_T^{\prime }x\right|\ge \frac{1}{2}\left|x\right|\text{ for all x∈RT}.$$

$$\frac{1}{2}{B}_{V_1}\subset (P_{V^{\prime }}S)(B_{Tn}),$$

Full text

1 Introduction

A dynamic linear regression (West and Harrison,, 1997, Chapter 3), or non-stationary regression, is a situation where we are given a sequence of scalar observations $\left\{Y_{t}\right\}_{t\leq T}\subset\mathbb{R}$, and observation vectors $\left\{u_{t}\right\}_{t\leq T}\subset\mathbb{R}^{n}$ such that $Y_{t}=\left<X_{t},u_{t}\right>+z_{t}$ where $X_{t}\in\mathbb{R}^{n}$ is a regressor vector, and $z_{t}$ a random noise term. In contrast to a standard linear regression, the vector $X_{t}$ may change with time. One common objective for this problem is at time $T$, to estimate the trajectory of $X_{t}$ for $t\leq T$, given the observation vectors and observations, $\left\{u_{t}\right\}_{t\leq T}$, $\left\{Y_{t}\right\}_{t\leq T}$, and possibly to forecast $Y_{T+1}$ if $u_{T+1}$ is also known.

In this paper we model the problem as follows:

[TABLE]

where $\left<\cdot,\cdot\right>$ is the standard inner product on $\mathbb{R}^{n}$, $z_{t}$, the observation noise, are zero-mean sub Gaussian random variables, with variance $\eta^{2}$, and the process noise variables $h_{t}$ take values in $\mathbb{R}^{n}$, such that coordinates of $h_{t}$ are zero-mean sub Gaussian, independent, and have variance $\sigma^{2}$. All $h_{t}$ and $z_{t}$ variables are assumed to be mutually independent. The vectors $u_{t}$ are an arbitrary sequence in $\mathbb{R}^{n}$, and the observed, known, quantities at time $T$ are $\left\{Y_{t}\right\}_{t\leq T}$ and $\left\{u_{t}\right\}_{t\leq T}$.

The system (1)-(2) is a special case of a Linear Dynamical System (LDS). As is well known, when the parameters $\sigma,\eta$ are given, the mean-squared loss optimal forecast for $Y_{T+1}$ and estimate for $X_{T}$ are obtained by the Kalman Filter (Anderson and Moore,, 1979; Hamilton,, 1994; Chui and Chen,, 2017). In this paper we are concerned with estimators for $\sigma,\eta$, and finite sample complexity guarantees for these estimators.

Let us first make a few remarks about the particular system (1)-(2). First, as a natural model of time varying regression, this system is useful in a considerable variety of applications. We refer to West and Harrison, (1997), 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 $\sigma,\eta$ in (1)-(2) as a pure case of finding the optimal learning rate for $X_{t}$. Indeed, the Kalman filter equations for (1)-(2), are given by (3)-(4) below, where (3) describes the filtered covariance update and (4) the filtered state update. Here $\bar{x}_{t}$ is the estimated state, given the observations $Y_{1},\ldots,Y_{t}$, see West and Harrison, (1997).

[TABLE]

In particular, following (4), the role of $\sigma$ and $\eta$ may be interpreted as regulating how much the estimate of $\bar{x}_{t+1}$ is influenced, via the operator $\frac{C_{t+1}}{\eta^{2}}$, by the most recent observation and input $Y_{t+1},u_{t+1}$. Roughly speaking, higher values of $\sigma$ or lower values of $\eta$ would imply that the past observations are given less weight, and result in an overfit of the forecast to the most recent observation. On the other hand, very low $\sigma$ or high $\eta$ would make the problem similar to the standard linear regression, where all observations are given equal weight, and result in a lag of the forecast. See Figure 3 in Supplementary Material Section A for an illustration.

Finally, it is worth mentioning that the system (1)-(2) is closely related to the study of online gradient (OG) methods (Zinkevich,, 2003; Hazan,, 2016). In this field, assuming quadratic cost, one considers the update

[TABLE]

where $\alpha$ is the learning rate, and studies the performance guarantees of the forecaster $\left<\bar{x}_{t},u_{t+1}\right>$. Compared to (4), the update (5) is simpler, and uses a scalar rate $\alpha$ instead of the input-dependent operator rate C_{t+1}\big{/}\eta^{2} of the Kalman filter. However, due to the similarity, every domain of applicability of the OG methods is also a natural candidate for the model (1)-(2) and vice-versa. As an illustration, we compare the OG to Kalman filter based methods with learned $\sigma$,$\eta$ in the experiments section.

In this paper we introduce a new estimation algorithm for $\sigma,\eta$, termed STVE (Spectrum Thresholding Variance Estimator), and prove finite sample complexity bounds for it. In particular, our bounds are an explicit function of the parameters $T$ and $\left\{u_{t}\right\}_{t=1}^{T}$ for any finite $T$, and indicate that the estimation error decays roughly as $T^{-\frac{1}{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 2, most existing estimation methods for LDSs, such as subspace identification (van Overschee and de Moor,, 1996; Qin,, 2006), or improper learning (Anava et al.,, 2013; Hazan et al.,, 2017; Kozdoba et al.,, 2019), do not apply to the system (1)-(2), due to non-stationarity. On the other hand, the methods that do apply to (1)-(2) either lack guarantees, or have only asymptotic analysis which in addition relies strongly on Guassianity of the noises.

Moreover, our approach differs significantly from the existing methods. We show that the structure of equations (1)-(2) 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 $\sigma$ and $\eta$ 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: The related work is discussed in Section 2 and Section 3 contains the necessary definitions. In Section 4 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 5. In Section 6 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.

2 Literature

We refer to Chui and Chen, (2017); Hamilton, (1994); Anderson and Moore, (1979); Shumway and Stoffer, (2011) for a general background on LDSs, the Kalman Filter and maximum likelihood estimation.

Existing approaches to the variance estimation problem may be divided into three categories: (i) General methods for parameter identification in LDS, either via maximum likelihood estimation (MLE) (Hamilton,, 1994), or via subspace identification (van Overschee and de Moor,, 1996; Qin,, 2006). In particular, finite sample bounds for system identification were given in (Campi and Weyer,, 2005; Vidyasagar and Karandikar,, 2006) and in the recent work Tsiamis 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 auto-correlation function, such as the classical Mehra, (1970); Belanger, (1974), or for instance more recent Wang et al., (2017); Dunik et al., (2018). (iii) Improper Learning methods, such as Anava et al., (2013); Hazan et al., (2017); Kozdoba et al., (2019). In these approaches, one does not learn the LDS directly, but instead learns a model from a 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 (1)-(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, asymptotically stationary systems, where the observation operator ($u_{t}$ in our notation) is constant and the Kalman gain (or, equivalently $C_{t}u_{t}$ in eq. (3)) 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 $u_{t}$ 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 to system (1)-(2) .

Among the methods that do apply to (1)-(2) are the general MLE estimation, and some of the auto-correlation methods (Belanger,, 1974; Dunik et al.,, 2018). On one hand, both types of approaches may be applicable to systems apriori more general than (1)-(2). 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 as a result the whole method has no guarantees. The results in Belanger, (1974); Dunik et al., (2018) do have asymptotic consistency guarantees. However, these rely on some explicit and implicit assumptions about the system, the sequence $u_{t}$ in our case, which can not be easily verified. In particular, Belanger, (1974); Dunik et al., (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 $u_{t}$. 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 $\left\{u_{t}\right\}_{t=1}^{T}$.

In addition, as mentioned earlier, MLE methods by definition must assume that the noises are Gaussian (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 of the noises. Therefore, to achieve such an extension, one would only need to correspondingly extend the main probabilistic tool, the Hanson-Wright inequality (Hanson and Wright,, 1971; Rudelson et al.,, 2013, see also Section 4 and Supplementary Material Section E). One such extension, for vectors with the convex concentration property, was recently obtained in Adamczak, (2015).

3 Notation

We refer to Bhatia, (1997) and Vershynin, (2018) as general references on the notation introduced below, for operators and sub Gaussian variables, respectively.

Let $A:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ be an operator with a singular value decomposition $A=U\cdot Diag(\lambda_{1},\ldots,\lambda_{s})\cdot W$, where $s\leq\min\left\{m,n\right\}$ and $\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{s}>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 $A=U\cdot Diag(\lambda_{1},\ldots,\lambda_{s})\cdot W$). The Hilbert-Schmidt (Frobenius) norm is defined as $\left\|A\right\|_{HS}=\sqrt{\sum_{i=1}^{s}\lambda_{i}^{2}}$. The nuclear and the operator norms are given by $\left\|A\right\|_{nuc}=\sum_{i=1}^{s}\lambda_{i}$ and $\left\|A\right\|_{op}=\lambda_{1}$ respectively.

A centered ($\mathbb{E}X=0$) scalar random variable $X$ is sub-Gaussian with constant $\kappa$, denoted $X\sim SG(\kappa)$, if for all $t>0$ it satisfies $\mathbb{P}\left(\left|X\right|>t\right)\leq 2\exp\left(-\frac{t^{2}}{\kappa^{2}}\right)$. A random vector $X=(X_{1},\ldots,X_{m})$ is $\kappa$ sub-Gaussian, denoted $X\sim SG_{m}(\kappa)$, if for every $v\in\mathbb{R}^{m}$ with $\left|v\right|=1$ the random variable $\left<v,X\right>$ is $\kappa$ sub-Gaussian. A random vector $X$ is $\sigma$-isotropic if for every $v\in\mathbb{R}^{m}$ with $\left|v\right|=1$, $\mathbb{E}\left<v,X\right>=\sigma^{2}$.

Finally, a random vector $X=(X_{1},\ldots,X_{m})$ is $\sigma$-isotropically $\kappa$ sub-Gaussian with independent components, denoted $X\sim ISG_{m}(\sigma,\kappa)$ if $X_{i}$ are independent, and for all $i\leq m$, $\mathbb{E}X_{i}=0$, $\mathbb{E}X_{i}^{2}=\sigma^{2}$ and $X_{i}\sim SG(\kappa)$. Clearly, if $X\sim ISG_{m}(\sigma,\kappa)$ then $X$ is $\sigma$-isotropic. Recall also that $X\sim ISG_{m}(\sigma,\kappa)$ implies $X\sim SG_{m}(\kappa)$ (Vershynin,, 2018). The noise variables we discuss in this paper are $ISG(\kappa,\sigma)$.

Throughout the paper, absolute constants are denoted by $c,c^{\prime},c^{\prime\prime},\ldots$. etc. Their values may change from line to line.

4 Overview of the approach

We begin by rewriting (1)-(2) in a vector form. To this end, we first encode sequences of $T$ vectors in $\mathbb{R}^{n}$, $\left\{a_{t}\right\}_{t\leq T}\subset\mathbb{R}^{n}$, as a vector $a\in\mathbb{R}^{Tn}$, constructed by concatenation of $a_{t}$’s. Next, we define the summation operator $S^{\prime}:\mathbb{R}^{T}\rightarrow\mathbb{R}^{T}$ which acts on any vector $(h_{1},h_{2},\ldots,h_{T})\in\mathbb{R}^{T}$ by

[TABLE]

Note that $S^{\prime}$ is an invertible operator. Next, we similarly define the summation operator $S:\mathbb{R}^{Tn}\rightarrow\mathbb{R}^{Tn}$, an $n$-dimensional extension of $S^{\prime}$, which sums $n$-dimensional vectors. Formally, for $(h_{l})_{l=1}^{Tn}\in\mathbb{R}^{Tn}$, and for $1\leq j\leq n,1\leq t\leq T$, $(Sh)_{(t-1)\cdot n+j}=\sum_{i\leq t}h_{(i-1)\cdot n+j}$. Observe that if the sequence of process noise terms $h_{1},\ldots,h_{T}\in\mathbb{R}^{n}$ is viewed as a vector $h\in\mathbb{R}^{Tn}$, then by definition $Sh$ is the $\mathbb{R}^{Tn}$ encoding of the sequence $X_{t}$.

Next, given a sequence of observation vectors $u_{1},\ldots,u_{T}\in\mathbb{R}^{n}$, we define the observation operator $O_{u}:\mathbb{R}^{Tn}\rightarrow\mathbb{R}^{T}$ by $(O_{u}x)_{t}=\left<u_{t},\left(x_{(t-1)\cdot n+1},\ldots,x_{(t-1)\cdot n+n}\right)\right>$. In words, coordinate $t$ of $O_{u}x$ is the inner product between $u_{t}$ and $t$-th part of the vector $x\in\mathbb{R}^{Tn}$. Define also $Y=(Y_{1},...,Y_{T})\in\mathbb{R}^{T}$ to be the concatenation of $Y_{1},...,Y_{T}$. With this notation, one may equivalently rewrite the system (1)-(2) as follows:

[TABLE]

where $h$ and $z$ are independent zero-mean random vectors in $\mathbb{R}^{Tn}$ and $\mathbb{R}^{T}$ respectively, with independent sub Gaussian coordinates. The variance of each coordinate of $h$ is $\sigma^{2}$ and each coordinate of $z$ has variance $\eta^{2}$.

Up to now, we have reformulated our data model as a single vector equation. Note that in that equation, the observations $Y$ and both operators $O_{u}$ and $S$ are known to us. Our problem may now be reformulated as follows: Given $Y\in\mathbb{R}^{T}$, assuming $Y$ was generated by (7), provide estimates of $\sigma,\eta$.

As a motivation, we first consider taking the expectation of the norm squared of eq. (7). For any operator $A:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ and zero-mean vector $h$ with independent coordinates and coordinate variance $\sigma^{2}$, we have $\mathbb{E}\left|Ah\right|^{2}=\left\|A\right\|_{HS}^{2}\sigma^{2}$, where $\left\|A\right\|_{HS}$ is the Hilbert-Schmidt (or Frobenius) norm of $A$. Taking the norm and expectation of (7), and dividing by $T^{2}$, we thus obtain

[TABLE]

Next, note that $\left\|O_{u}S\right\|_{HS}^{2}$ is known, and an elementary computation shows that $\frac{\left\|O_{u}S\right\|_{HS}^{2}}{T^{2}}$ is of constant order (as a function of $T$; see (25)), while the coefficient of $\eta^{2}$ is $\frac{1}{T}$. Thus, if the quantity $\frac{\left|Y\right|^{2}}{T^{2}}$ were close enough to its expectation with high probability, we could take this quantity as a (slightly biased) estimator of $\sigma^{2}$. However, as it will become apparent later, the deviations of $\frac{\left|Y\right|^{2}}{T^{2}}$ around the expectation are also of constant order, and thus $\frac{\left|Y\right|^{2}}{T^{2}}$ can not be used as an estimator. The reason for these high deviations of $\frac{\left|Y\right|^{2}}{T^{2}}$ is that the spectrum of $O_{u}S$ is extremely peaked. The highest squared singular value of $O_{u}S$ is of order $T^{2}$, the same order as sum of all of them, $\left\|O_{u}S\right\|_{HS}^{2}$. Contrast this with the case of identity operator, $Id:\mathbb{R}^{Tn}\rightarrow\mathbb{R}^{Tn}$: We have $\mathbb{E}\left|Id(h)\right|^{2}=\mathbb{E}\left|h\right|^{2}=Tn\sigma^{2}$, and one can also easily show that, for instance, $Var\left|Id(h)\right|^{2}=Tn\sigma^{2}$, and thus the deviations are of order $\sqrt{Tn}\sigma$ – a smaller order than $\mathbb{E}\left|Id(h)\right|^{2}$. While for the identity operator the computation is elementary, for a general operator $A$ the situation is significantly more involved, and the bounds on the deviations of $\left|Y\right|^{2}$ will be obtained from the Hanson-Wright inequality (Hanson and Wright,, 1971, see also Rudelson et al., (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_{u}S$ by taking the pseudo-inverse. Let $R:\mathbb{R}^{T}\rightarrow\mathbb{R}^{Tn}$ be the pseudo-inverse, or Moore-Penrose inverse of $O_{u}S$. Specifically, let

[TABLE]

be the singular value decomposition of $O_{u}S$, where $\gamma_{1}\geq\gamma_{2}\geq\ldots\geq\gamma_{T}$ are the singular values.

For the rest of the paper, we will assume that all of the observation vectors $u_{t}$ 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 $S$ is invertible and $O_{u}$ has rank $T$, we have $\lambda_{t}>0$ for all $t\leq T$. For $i\leq T$, denote $\chi_{i}=\gamma_{T+1-i}^{-1}$. Then $\chi_{i}$ are the singular values of $R$, arranged in a non-increasing order, and we have by definition

[TABLE]

where $W^{*},U^{*}$ denote the transposed matrices of $U,V$, defined in (9).

Similarly to Eq. (8), we apply $R$ to (7), and since $\left\|RO_{u}S\right\|_{HS}=T$, by taking the expectation of the squared norm we obtain

[TABLE]

In this equation, the deviation term $\left(\frac{\left|RY\right|^{2}}{T}-\frac{\mathbb{E}\left|RY\right|^{2}}{T}\right)$ is of order $O(\frac{1}{\sqrt{T}})$ with high probability (Theorem 1). Moreover, the coefficient of $\sigma^{2}$ is $1$, and the coefficient of $\eta^{2}$, which is $\frac{\left\|R\right\|_{HS}^{2}}{T}$, is of order at least $\Omega(\frac{1}{\log^{2}T})$ (Theorem 3, see Section 5 for additional details) – much larger order than $\frac{1}{\sqrt{T}}$. Since $\left|RY\right|^{2}$ and $\left\|R\right\|_{HS}^{2}$ are known, it follows that we have obtained one equation satisfied by $\sigma^{2}$ and $\eta^{2}$ up to an error of $\frac{1}{\sqrt{T}}$, where both coefficients are of order larger than the error.

Next, we would like to obtain another linear relation between $\sigma^{2}$, $\eta^{2}$. To this end, choose some $p=\alpha T$, where $0<\alpha<1$ is of constant order. The possible choices of $p$ are discussed later in this section. We define an operator $R^{\prime}:\mathbb{R}^{T}\rightarrow\mathbb{R}^{Tn}$ to be a version of $R$ truncated to the first $p$ singular values. If (10) is the SVD decomposition of $R$, then

[TABLE]

Similarly to the case for $R$, we have

[TABLE]

The deviations in (12) are also described by Theorem 1. Note also that since $\left\|R^{\prime}\right\|_{HS}^{2}$ is the sum of $p$ largest squared singular values of $R$, by definition it follows that $\frac{\left\|R^{\prime}\right\|_{HS}^{2}}{p}\geq\frac{\left\|R\right\|_{HS}^{2}}{T}$.

Now, given two equations in two unknowns, we can solve the system to obtain the estimates $\widehat{\sigma^{2}}$ and $\widehat{\eta^{2}}$. The full procedure is summarized in Algorithm 1, and the bounds implied by Theorem 1 on the estimators $\widehat{\sigma^{2}}$ and $\widehat{\eta^{2}}$ are given in Corollary 2. We first state these results, and then discuss in detail the various parameters appearing in the bounds.

Theorem 1**.**

Consider a random vector $Y\in\mathbb{R}^{T}$ of the form $Y=O_{u}Sh+z$ where $h\sim ISG_{Tn}(\sigma,\kappa)$ and $z\sim ISG_{T}(\eta,\kappa)$. Set $\left|u_{min}\right|=\min_{t}\left|u_{t}\right|$. Then for any $0<\delta<1$,

[TABLE]

where $B$ is given by

[TABLE]

The bounds on the estimators of Algorithm 1 are given in the following Corollary. As discussed below, in addition to Theorem 1, the key to the derivation of this Corollary are the estimates of the spectrum of $R$, given in Theorem 3, Section 5.

Corollary 2**.**

Let $\widehat{\sigma^{2}},\widehat{\eta^{2}}$ be the estimators of $\sigma^{2},\eta^{2}$ obtained form Algorithm 1 with $p\geq\frac{1}{4}T$. Set $\left|u_{max}\right|=\max_{t}\left|u_{t}\right|$. Then for any $0<\delta<1$, with probability at least $1-8\delta$,

[TABLE]

with $B$ given by (15).

We first discuss the assumption $\left|u_{min}\right|>0$. This assumption is made solely for notational convenience, as detailed below. To begin, note that some form of lower bound on the norms of the observation vectors $u_{t}$ must appear in the bounds. This is simply because if one had $u_{t}=0$ for all $T$, then clearly no estimate of $\sigma$ would have been possible. On the other hand, our use of the smallest value $\left|u_{\min}\right|$ may seem restrictive at first. We note however, that instead of considering the observation operator $O_{u}:\mathbb{R}^{Tn}\rightarrow\mathbb{R}^{T}$, one may consider the operator $O_{\bar{u}}:\mathbb{R}^{Tn}\rightarrow\mathbb{R}^{\bar{T}}$ for any subsequence $\left\{\bar{u}_{\bar{t}}\right\}_{\bar{t}=1}^{\bar{T}}$. The observation vector $Y$ would be correspondingly restricted to the subsequence of indices. This allows us to treat missing values and to exclude any outlier $u_{t}$ with small norms. All the arguments in Theorems 1 and 3 hold for this modified $O_{\bar{u}}$ without change. The only price that will be paid is that $T$ will be replaced by $\bar{T}$ in the bounds. Moreover, we note that typically we have $\left|u_{t}\right|\geq 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, (11)-(12), in two unknowns, $\sigma^{2},\eta^{2}$. Note that in order to be able to obtain $\eta^{2}$ from these equations, at least one of the coefficients of $\eta^{2}$, either $\frac{\left\|R\right\|_{HS}^{2}}{T}$ or $\frac{\left\|R^{\prime}\right\|_{HS}^{2}}{p}$ must be of larger order than $\frac{1}{\sqrt{T}}$, the order of deviations. Providing lower 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 to study the spectrum of $R$, and our approach proceeds indirectly via the analysis of the nuclear norm of $O_{u}S$. These results are stated in Theorem 3. In particular, we show that $\frac{\left\|R\right\|_{HS}^{2}}{T}$ is $\Omega(\frac{1}{\log^{2}T})$, which is the source of the log factor in (17).

Finally, in order to solve the equations (11)-(12), not only the equations must have large enough coefficients, but the equations must be different. This is reflected by the term $\left(1-\frac{p\left\|R\right\|_{HS}^{2}}{T\left\|R^{\prime}\right\|_{HS}^{2}}\right)^{-1}$ in (16), (17). Equivalently, while $\frac{\left\|R^{\prime}\right\|_{HS}^{2}/p}{\left\|R\right\|_{HS}^{2}/T}\geq 1$ by definition, we would like to have

[TABLE]

for the bounds (16), (17) to be stable. Note that since both $\left\|R^{\prime}\right\|_{HS}^{2}$ and $\left\|R\right\|_{HS}^{2}$ are computed in Algorithm 1, the condition (18) can simply be verified before the estimators $\widehat{\sigma^{2}},\widehat{\eta^{2}}$ are returned.

It is worth emphasizing that for simple choices of $p$, say $p=\frac{1}{4}T$, the condition (18) does hold in practice. Note that, for any $p$, we can have $\frac{\left\|R^{\prime}\right\|_{HS}^{2}/p}{\left\|R\right\|_{HS}^{2}/T}=1$ only if the spectrum of $R$ is constant. Thus (18) amounts to stating that the spectrum of $R$ exhibits some decay. As we show in experiments below, the spectrum (squared) of $R$, for $u_{t}$ derived from daily temperature features, or for random Gaussian $u_{t}$, indeed decays. See Section 6, Figures 1(a) and 1(b). In particular, in both cases (18) holds with $const>1$. Additional bounds on the quantity $\left(1-\frac{p\left\|R\right\|_{HS}^{2}}{T\left\|R^{\prime}\right\|_{HS}^{2}}\right)^{-1}$ under various assumptions on the sequence $u_{t}$ are given in Section G of the Supplementary Material.

5 Properties of $O_{u}S$ and $R$

As discussed in Section 4 (see the discussion following eq. (11)), one of the crucial points enabling Algorithm 1 and its analysis is the fact that the quantity $\frac{\left\|R\right\|_{HS}^{2}}{T}$ is bounded below by an expression that is of much higher order than the noise magnitude $\frac{1}{\sqrt{T}}$.

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_{u}S$ (Lemma 4, Supplementary Material Section B). Recall that the nuclear norm was defined in Section 3, and that for a sequence $u_{t}$ we set $\left|u_{max}\right|=\max_{t}\left|u_{t}\right|$ and $\left|u_{min}\right|=\min_{t}\left|u_{t}\right|$. Then:

[TABLE]

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 of the Laplacian are known precisely, satisfying $\lambda_{l}=2\sin\left(\frac{\pi(T-l)}{2T}\right)$. Next, we state the lower (and upper) bounds for $R$.

Theorem 3**.**

Let $R:\mathbb{R}^{T}\rightarrow\mathbb{R}^{Tn}$ be the pseudoinverse of $O_{u}S$. Then

[TABLE]

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 3 is nevertheless obtaining lower bounds on $\left\|R\right\|_{HS}^{2}$. Our approach is rather indirect, and we obtain these bounds from the nuclear norm bound (19) via a Markov type inequality on the eigenvalues.

6 Experiments

6.1 Synthetic Data

In this section the performance of STVE is evaluated on synthetic data. The data was generated by the LDS (1)-(2), using Gaussian noises with $\sigma^{2}=0.5,\eta^{2}=2$. The input dimension was $n=5$, and the input sequence $u_{t}$ sampled from the Gaussian $N(0,I_{n})$.

We run the STVE algorithm for different values of $T$, where for each $T$ we sampled the data $150$ times. Figure 2(a) shows the average (over 150 runs) estimation error for both process and observation noise variances for various values of $T$. As expected from the bounds in Corollary 2, it may be observed in Figure 2(a) that the estimation errors decay roughly at the rate of $const/\sqrt{T}$. A typical spectrum of $R$ is shown in Figure 1(b). For larger $T$, the spectra also exhibits similar decay.

6.2 Temperatures and Electricity Consumption

In this section we examine the relation between daily temperatures and electricity consumption in the data from Hong et al., (2014) (see also Hong, (2016)). The following forecasting methods are compared: a stationary regression, an online gradient, and a Kalman filter for a dynamic regression, with parameters learned via MLE or STVE. We find that the Kalman filter methods provide the best performance, with no significant difference between STVE and MLE derived systems.

The data consists of total daily electricity consumption (load) $y_{t}$, and the average daily temperature, $v_{t}$, 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 I. 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 H.

An elementary inspection of the data reveals that the load may be reasonably approximated as a quadratic function of the temperature, $y_{t}=x_{t,1}\cdot 1+x_{t,2}\cdot v_{t}+x_{t,3}\cdot v_{t}^{2}$, where $u_{t}=(1,v_{t},v_{t}^{2})$ is the observation vector (features), and $x_{t}=(x_{t,1},x_{t,2},x_{t,3})$ is the possibly time varying regression vector. This is shown in Figure 2(b), where we fit a stationary (time invariant) regression of the above form, using either only the first or 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.

We use the first half of the data (train set) to learn the parameters $\sigma,\eta$ of the online regression (1)-(2) via MLE optimization and using STVE. We also use the train set to find the optimal learning rate $\alpha$ 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 on the train set. In addition, we learn a 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 (1)-(2) are made with a Kalman filter (at time $t$, we use the filtered state estimate $\tilde{x}_{t}$, which depends only on $y_{1},\ldots,y_{t}$ and $u_{1},\ldots,u_{t}$, and make the prediction $\tilde{y}_{t+1}=\left<\tilde{x}_{t},u_{t+1}\right>$).

Daily squared prediction errors (that is, $\left(y_{t}-\tilde{y}_{t}\right)^{2}$) are shown in Figure 2(c) (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 $\left|u_{min}\right|\geq 1$ in this experiment, due to the constant $1$ coordinate, and also $\left|u_{max}\right|\leq 5$, due to the normalization. Since these operations are typical for any regression problem, we conclude that the direct influence of $\left|u_{min}\right|$ and $\left|u_{max}\right|$ on the bounds in Theorem 1 and Corollary 2 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 (1)-(2), and obtained the first sample complexity guarantees for such estimators. We have also shown how the shape of the spectrum of $R$ 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 1 and demonstrated in Section 6, the system (1)-(2) is of independent interest in applications. 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).

Checklist

For all authors…

1. (a)

Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope? [Yes] 2. (b)

Did you describe the limitations of your work? [Yes] 3. (c)

Did you discuss any potential negative societal impacts of your work? [N/A] The main focus of this work is a theoretical analysis. 4. (d)

Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] 2. 2.

If you are including theoretical results…

1. (a)

Did you state the full set of assumptions of all theoretical results? [Yes] 2. (b)

Did you include complete proofs of all theoretical results? [Yes] Thus Supplementary Materail contains all the proofs. Proofs outlines are given in the main body of the paper. 3. 3.

If you ran experiments…

1. (a)

Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] The data is publically available, see the references in Section 6. The short code is not provided at the moment, but can be fully derived from Algorithm 1. 2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Sections 6 and I. 3. (c)

Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [Yes] 4. (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. 4.

If you are using existing assets (e.g., code, data, models) or curating/releasing new assets…

1. (a)

If your work uses existing assets, did you cite the creators? [Yes] 2. (b)

Did you mention the license of the assets? [N/A] The link to the data author’s documentation is provided. 3. (c)

Did you include any new assets either in the supplemental material or as a URL? [No] 4. (d)

Did you discuss whether and how consent was obtained from people whose data you’re using/curating? [N/A] 5. (e)

Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] 5. 5.

If you used crowdsourcing or conducted research with human subjects…

1. (a)

Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A] 2. (b)

Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A] 3. (c)

Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]

Finite Sample Analysis Of Dynamic Regression Parameter Learning - Supplementary Material

Appendix A Outline

This Supplementary Material is organized as follows: The proofs of the results of Section 5, including Theorem 3, as well as proofs of Theorems 1, and Corollary 2, are given in Sections C to F. In Section G we state and prove additional bounds on the quantity $\left(1-\frac{p\left\|R\right\|_{HS}^{2}}{T\left\|R^{\prime}\right\|_{HS}^{2}}\right)^{-1}$ that appears in Corollary 2. Section H contains the discussion of missing values in STVE. Additional details on the electricity consumption experiment are given in Section I.

Finally, we describe Figure 3. The data (observations, red) was generated from the system (1)-(2) with one dimensional state ($n=1$), and we set $u_{t}=1$. The states produced by the Kalman filter with ground truth values of $\sigma,\eta$ are shown in green, while states obtained from other choices of parameters are shown in black and blue.

Appendix B Bounds on $O_{u}S$, Lemma 4

In this Section we prove the following bounds on the spectrum of $O_{u}S$. See Section 5 for a discussion of these bounds.

Lemma 4**.**

The singular values of $O_{u}S$ satisfy the following:

[TABLE]

The proof of this lemma uses the following auxiliary result: Let $D_{T}:\mathbb{R}^{T}\rightarrow\mathbb{R}^{T-1}$ be the difference operator, (D_{T}x)_{t}=x_{t+1}-x_{t}\mbox{ for t\leq T-1}. In the field of Finite Difference methods, the operator $D_{T}D_{T}^{*}$ is known as the Laplacian on the line, or as the discrete derivative with Dirichlet boundary conditions, and is well studied. The eigenvalues of $D_{T}D_{T}^{*}$ may be derived by a direct computation, and correspond to the roots of the Chebyshev polynomial of second kind of order $T$. In particular, the following holds:

Lemma 5**.**

The operator $D_{T}$ has kernel of dimension $1$ and singular values $\lambda_{l}(D_{T})=2\sin\left(\frac{\pi(T-l)}{2T}\right)$ for $l=1,\ldots,T-1$.

We refer to Mitchell and Griffiths, (1980) for the proof of Lemma 5. Next, in Lemma 6, we show that the inverse of the operator $S^{\prime}$, defined in (6), is a one dimensional perturbation of $D_{T}$, which implies bounds on singular values of $S^{\prime-1}$.

Lemma 6**.**

The singular values of $S^{\prime-1}_{T}$ satisfy

[TABLE]

for $1\leq t\leq T-1$, and

[TABLE]

The proof of this is given in the next section. With the estimates of Lemma 6, the bounds in Lemma 4 follow. Proof of Lemma 4:

Proof.

By Lemma 6,

[TABLE]

Note also that $S$ is by definition a collection of $n$ independent copies of $S^{\prime}_{T}$, and therefore the spectrum of $S$ is that of $S^{\prime}_{T}$, but each singular value is taken with multiplicity $n$. In particular it follows that (28) holds also for $S$ itself. Since clearly $\left\|O_{u}\right\|_{op}\leq\left|u_{max}\right|$, the upper bound on $\left\|O_{u}S\right\|_{op}$ in (23) follows from (28).

For the lower bound, denote by $V^{\prime}$ the orthogonal complement to the kernel of $O_{u}$, $V^{\prime}=\left(Ker(O_{u})\right)^{\perp}$. Denote by $P_{V^{\prime}}:\mathbb{R}^{Tn}\rightarrow V^{\prime}$ the orthogonal projection onto $V^{\prime}$. We have in particular that $O_{u}S=O_{u}P_{V^{\prime}}S$. Next, the operator $S$ maps the unit ball $B_{Tn}$ of $\mathbb{R}^{Tn}$ into an ellipsoid $\mathcal{E}_{T}$, and by (28), we have $\frac{1}{2}B_{Tn}\subset\mathcal{E}_{T}$. It therefore follows that

[TABLE]

where $B_{V_{1}}$ is the unit ball of $V_{1}$. It remains to observe that for every $x\in V_{1}$, we have

[TABLE]

Combining (29) and (30), we obtain the lower bound in (23).

To derive (24), recall that the nuclear norm is sub-multiplicative with respect to the operator norm:

[TABLE]

This follows for instance from the characterization of the nuclear norm as trace dual of the operator norm (Bhatia,, 1997, Propositions IV.2.11, IV.2.12). Next, since the spectrum of $S$ is the spectrum of $S^{\prime}$ taken with multiplicity $n$, we have

[TABLE]

and it remains to bound the nuclear norm of $S^{\prime}$.

Using the inequality

[TABLE]

and Lemma 6, we have

[TABLE]

Combining (31), (32) and (34), the inequality (24) follows.

It remains to estimate the Hilbert-Schmidt norm of $O_{u}S$, which can be done by a direct computation. Recall that for any operator $A:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ the Hilbert-Schmidt norm satisfies

[TABLE]

for any orthonormal basis $\phi_{i}$ in $\mathbb{R}^{m}$. Let $e_{ti}$ be the standard basis vector in $\mathbb{R}^{Tn}$ corresponding to coordinate $i\leq n$ at time $t\leq T$. Let $e_{t^{\prime}}$, $t^{\prime}\leq T$ denote the standard basis in $\mathbb{R}^{T}$. Then

[TABLE]

where $u_{t^{\prime}i}$ is the $i$-th coordinate of $u_{t^{\prime}}$. It follows that

[TABLE]

and hence

[TABLE]

The bounds (25) follow directly from (38). ∎

Appendix C Proof of Lemma 6

Proof.

Recall that the operator $S^{\prime-1}_{T}$ is given by $S^{\prime-1}_{T}x=(x_{1},x_{2}-x_{1},\ldots,x_{T}-x_{T-1})$ and the operator $D=D_{T}:\mathbb{R}^{T}\rightarrow\mathbb{R}^{T-1}$ is given by $Dx=(x_{2}-x_{1},\ldots,x_{T}-x_{T-1})$. Let $V=\text{span}\left\{e_{2},\ldots,e_{T}\right\}$ be the subspace spanned by all but the first coordinate. Let $P_{V}:\mathbb{R}^{T}\rightarrow\mathbb{R}^{T}$ be the projection onto $V$, i.e. a restriction to second to $T$’th coordinate. Observe that the action of $DP_{V}$ is equivalent to that of $S^{\prime-1}_{T-1}$. Therefore the singular values of $S^{\prime-1}_{T-1}$ are identical to those of $DP_{V}$. To obtain bounds on the singular values of $DP_{V}$, note that $(DP_{V})^{*}DP_{V}=P_{V}D^{*}DP_{V}$ – that is, $P_{V}D^{*}DP_{V}$ is a compression of $D^{*}D$. Thus, by the Cauchy’s Interlacing Theorem (Bhatia,, 1997, Corollary III.1.5),

[TABLE]

for all $t\leq T-1$. In conjunction with Lemma 5 this provides us with the estimates for all but the smallest singular value of $S^{\prime-1}_{T-1}$ (since $\lambda_{T}(D^{*}D)=0$). We therefore estimate $\lambda_{T-1}(S^{\prime-1}_{T-1})=\lambda_{T}(DP_{V})$ directly, by bounding the norm of $S^{\prime}_{T-1}$. Indeed, for any $T$, by the Cauchy-Schwartz inequality,

[TABLE]

Thus we have $\left\|S^{\prime}_{T}\right\|_{op}\leq T$, which concludes the proof of the Lemma. ∎

Appendix D Proof of Theorem 3

Proof.

The bound on $\left\|R\right\|_{op}$ follows directly from the lower bound on the singular values of $O_{u}S$ in (23). Since $R$ is of rank $T$, the upper bounds on $\left\|R\right\|_{HS}$ follow directly from the $\left\|R\right\|_{op}$ bound.

The lower bounds on $\left\|R\right\|_{HS}^{2}$ and $\left\|R^{*}R\right\|_{HS}^{2}$ follow from the upper bounds on the nuclear norm in Lemma 4. Note that this argument would not have worked if we only had upper bounds on the Hilbert-Schmidt norm, rather than the nuclear norm in Lemma 4. Denote $\gamma_{i}=\lambda_{i}(O_{u}S)$. From (24) in Lemma 4, the number of $\gamma_{i}$ that are larger than $4n\left|u_{max}\right|\log T$ satisfies

[TABLE]

Since there are total of $T$ singular values $\gamma_{i}$ overall, we can equivalently rewrite (42) as

[TABLE]

This immediately implies the lower bounds on $\left\|R\right\|_{HS}^{2}$ and $\left\|R^{*}R\right\|_{HS}^{2}$. ∎

Appendix E Proof of Theorem 1

The two main probabilistic tools that we use are the Hanson-Wright inequality (Hanson and Wright,, 1971), and a classical norm deviation inequality for sub Gaussian vectors, as follows:

Theorem 7** (Hanson-Wright Inequality).**

Let $X=(X_{1},\ldots,X_{m})\in\mathbb{R}^{m}$ be a random vector such that the components $X_{i}$ are independent and $X_{i}\sim SG(\kappa)$ for all $i\leq m$. Let A be an $m\times m$ matrix. Then, for every $t\geq 0$,

[TABLE]

In particular, recall that we are interested in concentration of $\left|RY\right|^{2}$. We may write:

[TABLE]

The deviations of the first two terms may be bounded via Theorem 7. For the third term, we use the following:

Lemma 8**.**

For any $X\sim SG_{m}(\kappa)$ we have:

[TABLE]

Note that in Lemma 8 we do not require the coordinates of $X$ to be independent. This will be important in what follows. Lemma 8 is standard and can be proved via covering number estimates of the Euclidean ball, see for instance Vershynin, (2018), Section 4.4.

In addition, the following observation is used throughout the text:

Lemma 9**.**

Let $A:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ be an operator, and let $h=(h_{1},\ldots,h_{m})$ have independent coordinates with $\mathbb{E}h_{i}^{2}=\sigma^{2}$. Denote by $\left\{\lambda_{i}\right\}_{i=1}^{k}$, $k\leq\min(m,n)$, the singular values of $A$. Then

[TABLE]

The elementary proof is omitted.

We now prove Theorem 1.

Proof.

Let $0<\delta<1$ be given. We first bound the deviations from the expectation for the first term in (45), $\left|RO_{u}Sh\right|^{2}$. We apply the Hanson-Wright inequality with $X=h$ and $A=(RO_{u}S)^{*}RO_{u}S$. By definition, $A$ has a single eigenvalue $1$ with multiplicity $T$. Thus clearly $\left\|A\right\|_{HS}^{2}=T$ and $\left\|A\right\|_{op}=1$.

For an appropriate constant $c^{\prime}>0$, set $t=c^{\prime}\kappa^{2}\sqrt{T}\log\frac{1}{\delta}$. Then,

[TABLE]

The deviation of the second term in (45) is similarly bounded using $A=R^{*}R$. Recall that by Theorem 3 we have $\left\|R\right\|_{op}\leq 2\left|u_{min}\right|^{-1}$, and note that $\left\|R^{*}R\right\|_{HS}^{2}\leq T\left\|R\right\|_{op}^{4}\leq cT\left|u_{min}\right|^{-4}$. Set $t=c^{\prime}\kappa^{2}\left|u_{min}\right|^{-2}\sqrt{T}\log\frac{1}{\delta}$. With this choice it follows that both terms in the minimum in (44) are larger than $c\log\frac{1}{\delta}$ and we have

[TABLE]

Finally, we bound the third term in (45). Denote by $D$ the event

[TABLE]

and by $E$ the event

[TABLE]

By Lemma 8 applied to $z$, and using the fact that $\left\|(RO_{u}S)^{*}R\right\|_{op}\leq 2\left|u_{min}\right|^{-1}$,

[TABLE]

Next, using independence of $h,z$ and $h\sim SG_{Tn}(\kappa)$,

[TABLE]

where $D^{c}$ is the complement of $D$. Therefore, combining (52) and (53),

[TABLE]

Combining (48), (49) and (54), we obtain via the union bound:

[TABLE]

Similarly, we obtain a bound for the equations involving $R^{\prime}$:

[TABLE]

The only difference in the derivation of (56) compared to (55) is in the application of Lemma 8. In the later case, to replace $\sqrt{T}$ with $\sqrt{p}$ in (52), we apply Lemma 8 with $X=P_{V}z$ rather than with $X=z$, where $P_{V}$ is the projection onto the range of $R^{\prime}$ – a $p$-dimensional space. Note that $P_{V}z$ does not necessarily have a structure of $p$ independent coordinates, but is sub Gaussian and isotropic. Therefore Lemma 8 still applies.

∎

Appendix F Proof of Corollary 2

We now turn to prove Corollary 2.

Proof.

Denote

[TABLE]

Using (56) and (55) we may write

[TABLE]

where $e_{1},e_{2}$ are error terms such that $\left|e_{1}\right|\leq E(\sqrt{T})$ and $\left|e_{2}\right|\leq E(\sqrt{p})$ holds with probability at least $1-8\delta$.

It follows that

[TABLE]

and

[TABLE]

Thus

[TABLE]

It remains to observe that

[TABLE]

where the inequality follows from eq. (43) in the proof of Theorem 3. ∎

Appendix G Thresholding Gap Analysis

The main result of this section is the following Proposition.

Proposition 10**.**

Given the sequence $\left\{u_{t}\right\}_{t=1}^{T}\subset\mathbb{R}^{n}$, define the scalar sequence $\tilde{u}_{t}=\sum_{i\leq n}u_{ti}$ and set

[TABLE]

Then for $p=\frac{1}{4}T$,

[TABLE]

The general idea behind the proof of Proposition 10 is to show that for $n=1$, the ratio $\left(1-\frac{p\left\|R\right\|_{HS}^{2}}{T\left\|R^{\prime}\right\|_{HS}^{2}}\right)^{-1}$ can be controlled. This is done in Lemmas 11 and 12 below. In particular, Lemma 11 is a general statement about integrals of monotone real functions under certain order constraints, and Lemma 12 provides a relation of the spectrum of $R$ to that of $S^{\prime-1}$. It is then shown that for arbitrary $n$, $O_{u}S$ contains a certain copy of an $n=1$-dimensional operator with parameters $\tilde{u}_{t}$, which implies the bounds.

For the case $n=1$, stated in Lemma 12, the argument consists of showing that the spectrum of $R$ is upper and lower bounded by appropriately decaying functions, and therefore can not be “too constant”. We first obtain general estimates for the integrals of such upper and lower bounded functions in the following Lemma:

Lemma 11**.**

Let $f:[0,1]\rightarrow\mathbb{R}$ be a monotone non-increasing function such that for all $x\in[0,1]$,

[TABLE]

for some $M\geq 1$. Set $t_{0}=\frac{1}{4}$. Then

[TABLE]

Proof.

Denote

[TABLE]

Write

[TABLE]

and set $v:=f(t_{0})$. Then, among all $f$ that satisfy (67) and $f(t_{0})=v$, $r(f)$ is minimized on $f$ with maximal $I(t_{0},1,f)$ and minimal $I(0,t_{0},f)$. Due to the form of the constraint (67) and monotonicity, this minimizer is given by

[TABLE]

where $t_{v}^{-}:=\max\left\{0,1-\sqrt{v}\right\}$ and $t_{v}^{+}:=1-\sqrt{\frac{v}{M}}$. Our problem is therefore now reduced from a minimization of $r(f)$ over the function space to a problem of minimizing $r(v):=r(\tilde{f_{v}})$ over a single scalar variable $v$.

To this end, first note that we can assume w.l.o.g that $M\geq(1-t_{0})^{-2}$. Indeed, for larger $M$, the lower bound on $r(f)$ can only become smaller, since $r(f)$ would be minimized over a larger set. By construction, the value $v$ must satisfy $(1-t_{0})^{2}\leq v\leq M(1-t_{0})^{2}$, and for $M\geq(1-t_{0})^{-2}$, the value $v=1$ satisfies these inequalities.

Next, by direct computation (taking the derivative in $v$) one verifies that for any $M\geq(1-t_{0})^{-2}$, $r(v)$ as a function of $v$ is minimized at $v=1$. It remains to observe that $r(1)=1-\frac{1}{3\sqrt{M}}$, which yields the statement of the Lemma. ∎

We can now treat the $n=1$ case.

Lemma 12**.**

Consider the case $n=1$, and $p=\frac{1}{4}T$. Then

[TABLE]

Proof.

For any two operators $A,B$ we have

[TABLE]

see Bhatia, (1997); Gohberg and Krein, (1969). Note that in the case $n=1$ and $\left|u_{min}\right|>0$, $O_{u}S$ is invertible. Thus the singular values of $O_{u}S^{\prime}$ satisfy

[TABLE]

where the first inequality follows by applying (73) to $\left(O_{u}S^{\prime}\right)^{-1}$, and the second by considering $O_{u}S^{\prime}$ itself. Equivalently, for $i=1,\ldots,T$

[TABLE]

and using Lemma 6,

[TABLE]

By changing the index and taking squares, we have

[TABLE]

Now, using an elementary discretization argument, for $p=\frac{1}{4}T$ by Lemma 11 it follows that

[TABLE]

which implies (72). ∎

Finally, we prove Proposition 10.

Proof.

As discussed earlier, the key to a statement such as (66) is to show that the spectrum of $R$ is non constant, and in particular has enough small singular values. This is equivalent to providing appropriate lower bounds on the spectrum of $O_{u}S$. Here we derive such bounds by comparison with an $n=1$ case as follows: Let $V\subset\mathbb{R}^{Tn}$ be a $T$ dimensional subspace spanned by vectors for which for every time $t$, all coordinates at time $t$ are identical. Formally, for $x\in V$, for every $t$ we require that $x_{(t-1)n+i}=x_{(t-1)n+j}$ for all $i,j\leq n$.

Observe that the restriction of $O_{u}S$ to $V$, the operator $O_{u}SP_{V}$, acts equivalently to the $n=1$ case operator defined by $O_{\tilde{u}}S^{\prime}$. In particular, similarly to the argument in Lemma 12, it follows that

[TABLE]

Moreover, note that $\left(O_{u}SP_{V}\right)^{*}O_{u}SP_{V}=P_{V}S^{*}O_{u}^{*}O_{u}SP_{V}$ and $\left(O_{u}S\right)^{*}O_{u}S$ are non-negative operators and $\left(O_{u}S\right)^{*}O_{u}S\geq\left(O_{u}SP_{V}\right)^{*}O_{u}SP_{V}$ (that is, $\left(O_{u}S\right)^{*}O_{u}S-\left(O_{u}SP_{V}\right)^{*}O_{u}SP_{V}$ is non-negative). It follows that for all $i\leq T$,

[TABLE]

The upper bounds on the spectrum $\lambda_{i}(O_{u}S)$ may again be obtained via (73). Indeed,

[TABLE]

where the first inequality follows from (73) while the second is due to the multiplicity $n$ of each singular value of $S^{\prime}$ in $S$. The rest of the argument proceeds as in Lemma 12. ∎

Appendix H Missing Values in STVE

We now discuss the treatment of missing values in STVE. Recall that our starting point is the vector form of the system, (7), which we rewrite here:

[TABLE]

where $h\in\mathbb{R}^{Tn}$ and $z\in\mathbb{R}^{T}$ are sub Gaussian vectors and $Y=(y_{1},\ldots,y_{T})\in\mathbb{R}^{T}$ is the observation vector. Suppose that $M$ out of $T$ observation values are missing, at times $t_{1},\ldots,t_{M}$. Set $T^{\prime}=T-M$ and define a projection operator $P_{A}:\mathbb{R}^{T}\rightarrow\mathbb{R}^{T}$ as the operator that omits the coordinates $t_{1},\ldots,t_{M}$. Formally, let $e(t)$, $t\leq T$, be the standard basis vector in $\mathbb{R}^{T}$, with $1$ at coordinate $t$ and zeros elsewhere. Then

[TABLE]

We can then rewrite (82) as

[TABLE]

Note that the vector $P_{A}Y$ contains only available, non-missing values of $y_{t}$. Similarly to the case with no missing values, we define $R$ as the Moore-Penrose inverse of $P_{A}O_{u}S$. Then we have

[TABLE]

Note that since a Moore-Penrose inverse of $R$ only acts on the image of $R$, we have $RP_{A}=R$ and therefore

[TABLE]

similarly to the case with no missing values. The only difference here will be that $R$ will now have $T^{\prime}$ rather than $T$ non-zero singular values. $R^{\prime}$ will be defined similarly. The analysis of the spectrum of $R$ concerns only non-zero singular values of $R$ and holds with no change other than that $T$ should be replaced by $T^{\prime}$. Consequently, the whole approach works identically with $T$ replaced by $T^{\prime}$ everywhere, and in the bounds of Theorems 1 and 3 in particular.

Appendix I Electricity Consumption and Temperatures Data

The raw data in Hong et al., (2014) contains electricity consumption levels for 20 nearby areas, also referred to as zones. Each zone had slightly different consumption characteristics. All figures in Section 6.2 refer to Zone 1. However, the results are similar for all zones. In particular, while Figure 2(c) shows the (smoothed) prediction errors of each method over time, the total error on the test set, normalized by the number of days, for each method and each zone, is shown in Figure 4(b). As with zone 1, adaptive methods perform better for all zones. STVE based estimator is better than MLE in half the cases, but the differences between STVE and MLE performance are negligible compared to errors in the stationary regression or OG.

In the rest of this section we discuss the preprocessing of the data. The raw data contains hourly consumption levels for the 20 zones, and hourly temperatures from 11 weather stations with unspecified locations in the region. Here we are only interested in the total daily consumption – the total consumption over all hours of the day – for each zone. For each station we also consider the average daily temperature for the day. Moreover, since the average daily temperatures at different stations are strongly correlated ($\rho\geq 0.97$ for all pairs of stations), we only consider a single daily number, $v_{t}$ – an average temperature across all hours of the day and all stations.

From the total of about 234 weeks in the data, 9 non-consecutive weeks have missing loads (consumption) data, 4 in the first half of the period (train set), and 5 in the second half. For the purposes of training the stationary regression, we excluded the data points with missing values from the train set. The STVE, MLE and OG methods were trained with missing values. The prediction results of all methods were evaluated only at the points where actual values were known.

The temperatures and all the loads are normalized to have zero-mean and standard deviation 1 on the first half of the data (train set).

The MLE estimates of the parameters $\sigma^{2},\eta^{2}$ were obtained using the DLM package of the R environment (Petris,, 2010). The package uses an L-BGFS-B algorithm with numerically approximated derivatives as the underlying optimization procedure.

One immediate reason for the dependence of the load on temperature to change with time (as shown in Figure 2(b) for Zone 1) is that the load simply grows with time (assuming the temperatures do not exhibit a trend of the same magnitude). Indeed, it is easy to verify that there is an upward trend in the overall load, as shown in Figure 4(a). However, the upward trend is non-uniform in temperature, since the parabolas in Figure 2(b) are not a shift of each other by a constant.

Finally, in Figure 4(b) average prediction errors are shown for each method and for each zone in the dataset. Zone 9 is known to be an industrial zone, (Hong,, 2016), where the load does not significantly depend on the temperature, hence the higher error rates for all the methods. For the other zones the situation is similar to Zone 1, with STVE and MLE performing comparably, and better than either OG or the Stationary regression.