# Cluster-and-Conquer: A Framework for Time-Series Forecasting

Reese Pathak1 Rajat Sen2 Nikhil Rao3 N. Benjamin Erichson4,5 Michael I.
Jordan1,5 Inderjit S. Dhillon3,6

(1Department of Electrical Engineering and Computer Sciences, UC Berkeley
2Google Research
3Amazon
4School of Engineering, University of Pittsburgh
5Department of Statistics, UC Berkeley
6Department of Computer Science, UT Austin
)

###### Abstract

We propose a three-stage framework for forecasting high-dimensional time-
series data. Our method first estimates parameters for each univariate time
series. Next, we use these parameters to cluster the time series. These
clusters can be viewed as multivariate time series, for which we then compute
parameters. The forecasted values of a single time series can depend on the
history of other time series in the same cluster, accounting for intra-cluster
similarity while minimizing potential noise in predictions by ignoring inter-
cluster effects. Our framework—which we refer to as “cluster-and-conquer”—is
highly general, allowing for any time-series forecasting and clustering method
to be used in each step. It is computationally efficient and embarrassingly
parallel.We motivate our framework with a theoretical analysis in an idealized
mixed linear regression setting, where we provide guarantees on the quality of
the estimates. We accompany these guarantees with experimental results that
demonstrate the advantages of our framework: when instantiated with simple
linear autoregressive models, we are able to achieve state-of-the-art results
on several benchmark datasets, sometimes outperforming deep-learning-based
approaches.

##  1 Introduction

High-dimensional time-series forecasting is crucial in applications such as
finance [45], e-commerce [17], medical data analysis [31], fluid flows and
climatology [13, 2]. Modern time-series datasets can have millions of
correlated time series; e.g., demand forecasting for an online store like
Amazon might require predicting the future sales of millions of items, and
these demands are generally inter-related. Thus time-series forecasting
methods should be scalable and equipped with the ability to handle inter-time-
series correlation.

Most time-series forecasting methods operate at the two extremes: (1)
Univariate modeling where the future prediction of a time series depends on
its own past and associated covariates; (2) Multivariate modeling where the
future prediction of a time series can potentially depend on the past of all
the time series in the dataset. Univariate models are unable to capture inter
time-series correlations during predictions, yet are computationally efficient
and/or easy to parallelize. On the other hand, multivariate time-series models
are better able to capture these time-series correlations, yet can suffer from
computational inefficiencies. We give a more extensive comparison of prior
work within these two extremes, in Section 2.

The disadvantages of univariate and multivariate time-series modeling point to
a need to come up with forecasting methods that can achieve a middle ground,
i.e, the future values of a time series are modeled as a function of its own
past and only a few other related time series. Therefore, in this paper we ask
the following question: _For high-dimensional time-series problems, is it
possible to develop statistically and computationally efficient methods that
exploit only relevant inter-time-series relationships?_. Prior works like [5]
have indeed shown some promising empirical results in this direction, but
detailed models of this problem with theoretical grounding is still lacking.

In this paper, we tackle this question in both theory and practice, and answer
it in the affirmative. We propose a ‘cluster-and-conquer’ framework. Given $n$
time series, we first use highly scalable univariate models (such as scalar
autoregression) to learn parameters of each time series. Then, we apply
efficient clustering algorithms on these parameters, to partition the
collection of time series into $k$ clusters. Finally, for each cluster, we
apply a joint forecasting approach (using all the time series only in the
cluster).

The main contributions of this paper are a framework for time-series modelling
with the following properties:

  * (i)

By varying the number of clusters, we are able to interpolate between
univariate and multivariate modelling, and reap the statistical and
computational benefits of each paradigm.

  * (ii)

When employed with standard linear autoregressive models, our method shows
strong empirical performance. Specifically, we are able to outperform many
previously proposed time-series modelling approaches over a wide variety of
datasets.

  * (iii)

Our framework—when analyzed in a mixture-of-regressions setting—enjoys
theoretical guarantees. Namely, we can show that this framework is able to
ensure recovery of model parameters at standard parametric rates.

Notation. A multivariate time-series dataset with $n$ time series and $T$ time
points can be represented as $\mathbf{X}\in\mathbf{R}^{n\times T}$. For
integers, $i\leqslant j$, $i:j:=\\{i,i+1,\cdots,j\\}$. For a positive integer
$n$, $[n]$ denotes the set $\\{1,2,\cdots,n\\}$. The notation
$x_{i}^{(\mathcal{J})}$ denotes the values at the time points in $\mathcal{J}$
for time series $i$ represented as a column vector, in ascending order of the
indices in $\mathcal{J}$. $\|\cdot\|_{p}$ denotes the $\ell_{p}$ norm for
vectors. The operator and Frobenius norms of matrices are denoted by
$\left\|\cdot\right\|_{\rm op}$ and $\left\|\cdot\right\|_{\rm F}$
respectively. For two vectors $u\in\mathbf{R}^{p},v\in\mathbf{R}^{q}$, we
write $u\otimes v$ for their tensor (outer) product, which is a $p\times q$
matrix with entries $u_{i}v_{j}$, $i\in[p],j\in[q]$.

\begin{overpic}[width=411.93767pt]{figs/overview.pdf}
\put(10.0,-1.0){\small{$1$}} \put(27.0,-1.0){\small{$2$}}
\put(44.0,-1.0){\small{$3$}} \put(84.0,-1.0){\small{$N=\left\lfloor
T/b\right\rfloor$}} \par\put(96.0,23.0){\small{$x_{i}^{s_{N}}$}}
\put(80.8,23.0){\small{$x_{i}^{s_{N}-d:s_{N}-1}$}}
\par\put(17.0,23.0){\small{$x_{i}^{s_{1}}$}}
\put(2.0,23.0){\small{$x_{i}^{s_{1}-d:s_{1}-1}$}}
\par\put(3.0,6.5){\small{$\widehat{c}(i)$}}
\par\put(10.0,5.5){\Large{$\vdots$}} \put(27.0,5.5){\Large{$\vdots$}}
\put(44.0,5.5){\Large{$\vdots$}} \put(89.0,5.5){\Large{$\vdots$}}
\put(65.0,6.5){\Large{$\ldots$}} \end{overpic} Figure 1: Illustration of
forming regression datasets from time-series data blocks. Since $b$ is a
multiple of the mixing time of the time series, the samples for the regression
problem will be approximately be independent and identically distributed.
Above, $\\{s_{i}\\}_{i=1}^{N}$ are the sample time points selected within each
chunk.

##  2 Related Work

There is a vast body of work in classical time-series forecasting; we refer
the reader to Hyndman and Athanasopoulos [21] for an overview. Here, we mainly
focus on the time-series forecasting literature in the machine learning
community; see also Benidis et al. [6].

#### Univariate modeling.

Basic forecasting techniques train one model per dataset but the prediction
using the trained model is done for each time series and its covariates
individually. Specifically, in univariate modeling, we typically use
predictions of the form

$\widehat{x}_{i}^{(T+1)}=\widehat{f}_{i}(x_{i}^{(1:T)}),\quad i=1,\dots,n.$

Autoregressive (AR) models are a classic example that are computationally
scalable and enjoy low variance, but these models typically suffer from high
statistical bias in that they cannot account for complicated inter-time-series
relationships. In contrast, recent deep learning methods [35, 8, 15, 29, 34]
aim to address this problem by training shared weights of a large model on the
whole dataset—effectively selecting $\widehat{f}_{i}$ based on the other time
series $\\{x^{(t)}_{i}\\}_{t\leqslant T,j\neq i}$. At prediction time the
estimate of future values of each univariate time series is modelled by a
function of its own past and covariates.

#### Multivariate modeling.

More advanced forecasting techniques aim to capture dependencies across
different time series in the dataset. Prediciton for each time series is
computed as a function of the history of _all_ the other time series.
Specifically, in multivariate modeling, we typically use predictions of the
form

$\widehat{x}_{i}^{(T+1)}=\widehat{f}_{i}(x_{1}^{(1:T)},,x_{2}^{(1:T)},\dots,x_{n}^{(1:T)}).$

Models that use this paradigm include dense vector autoregression (VAR),
dynamic linear models (DLMs), TRMF [43], DeepGLO [36] and graph neural network
based models [28, 9, 42, 40, 20, 3, 14]. However, some of these models are
computationally inefficient. For instance, VAR does not scale beyond even
moderately sized datasets. The complexity to fit a VAR model on a dataset of
$n$ time series scales as $O(n^{6})$ [43]. Similarly, DLMs typically scale as
$O(n^{3})$, which is again prohibitive for large-$n$ problems. Additionally,
these models typically have a large number of parameters, and can therefore
suffer from high model variance. At prediction time, forecasts can be made
using contributions from potentially irrelevant time series which can lead to
higher statistical error.

#### Time-series clustering.

The vast literature on time-series clustering is also relevant to our work. We
refer the reader to the surveys by Warren Liao [39] and Aghabozorgi et al. [1]
for an in-depth treatment. Model-based clustering techniques that assume that
similar time series are governed by similar model parameters are especially
relevant to our work [22, 12, 38]. In recent work [5], Bandara and co-authors
consider a similar cluster-based time-series prediction framework. However,
there are two main differences: first, the authors do not present theory to
justify or accompany their modelling framework. Secondly, their within-cluster
model does not forecast individual time series based on the entire history of
the other time series, and could miss important inter-cluster time-series
correlations.

##  3 Cluster-and-Conquer Framework

Suppose we are given a history of $T$ timepoints of $n$ univariate time
series, $\\{x_{i}^{(t)}\\}_{i\leqslant n,t\leqslant T}$. We would like to
produce a prediction, for the next time point, which is at time $T+1$.222For
simplicity we write all our equations for predicting one time point in the
future, however all the ideas trivially extend to multi-step forecasting. We
want our predictions to satisfy $\widehat{x}_{i}^{(T+1)}\approx
x_{i}^{(T+1)}$, for all time series $i=1,\dots,n$. Now, to overcome some of
the issues discussed in the related work section, we propose a three-stage
procedure for time-series forecasting.

Our framework consists of the following three steps:

  * •

Local step:  We first divide the time-series dataset into $N=\left\lfloor
T/b\right\rfloor$ chunks of contiguous times as shown in Figure 1; within each
of the $N$ chunks, we select a time point $s_{i},i\in[N]$, for which we sample
the data. At a high level, we can think of selecting $b$ as a multiple of the
(largest) mixing time of the time series. We can now extract one sample for a
local autoregressive regression problem per time series from each block. For
example for the $i$-th time series the sample from the $j$-th block can be

$(\text{covariates},\text{target})=(x^{(s_{j}-d:s_{j}-1)}_{i},x_{i}^{(s_{j})}).$

For each time series we collect all the samples and solve a least-squares
problem, and hence this step is efficient. The AR parameters obtained are
$\\{\widehat{\theta}^{\rm AR}_{i}\\}_{i=1}^{n}$.

  * •

Clustering step:  Next, we cluster the parameters of the AR model. We can use
a clustering algorithm, for instance spectral clustering [23], or METIS [22]
with inputs $\\{\widehat{\theta}^{\rm AR}_{i}\\}_{i=1}^{n}$ to obtain a
clustering of time series into $k$ clusters. Let $\widehat{\mathsf{c}}(i)$ be
the estimated cluster identity of time series $i$.

  * •

Global step:  Given a partition, we can estimate VAR parameters for every
cluster. We can reuse the block structure from above and extract one sample
per block for the VAR problem associated with each time series. For instance,
for the $i$-th time series the sample for the $j$-th block can be

covariates
$\displaystyle=x^{(s_{j}-d:s_{j}-1)}_{i},\\{x^{(s_{j}-d:s_{j}-1)}_{l}\\}_{\widehat{\mathsf{c}}(l)=\widehat{\mathsf{c}}(i)},$
target $\displaystyle=x^{(s_{j})}_{i}.$

The covariates are the history of the $i$th time series, together with all the
histories of the time series assigned to the same cluster—that is, time series
$l$ satisfying $\widehat{\mathsf{c}}(l)=\widehat{\mathsf{c}}(i)$. Therefore,
in this model the future of a time series can depend on the past of all time
series in the same cluster, through the learned VAR coefficients. Note that we
use VAR here for exposition and to make it amenable to theoretical analysis.
In practice, this step can be solved using more complex low-rank [43] or deep-
learning [36] approaches.

Each of the samples for the two regression problems above should be close to
i.i.d as $b$ is a multiple of the mixing time of the processes. Therefore we
will analyze an i.i.d analogue of the above system which is similar to a mixed
linear regression problem. The mixed linear regression problem is a useful
surrogate of our problem in an i.i.d setting that lets us analyze the salient
features of our framework—parameter recovery of local AR models, recovery of
clusters using clustering on the recovered AR parameters and recovery of VAR
parameters within cluster—while avoiding the complexities arising from
non-i.i.d data. Note that risk bounds and parameter recovery bounds for simple
AR models are still an area of active research [19, 32, 25] and we would like
to decouple the hardness of analyzing our framework from the issues arising
from a non-i.i.d dataset. These issues can be carefully dealt with using tools
from mixing time literature and are left for future work. Analogous i.i.d
settings have been used before in the literature to analyze interesting
properties of time-series clustering algorithms [38, 12].

##  4 Theoretical Results

In this section we provide our formal algorithm and associated parameter
recovery guarantees in the Mixed Linear Regression (MLR) setting. Although
idealized, this formal model allows us to state formal guarantees about
parameter recovery. Additionally, we refer the reader to Appendix D.1, where
we further explain how this setting is closely related to the original time-
series problem.

MLR model:  We consider the following fixed design mixed linear regression
problem. There are $n$ regression designs each with $T$ samples. $x_{i}^{(t)}$
denotes the $t$-th covariate from the $i$-th design. We will use the notation
$X^{(t)}=\\{x_{i}^{(t)}\\}_{i=1}^{n}$, where $X^{(t)}\in\mathbf{R}^{n\times
d}$.

We assume that for some positive integer $k$, we have (unknown) cluster
labels, $\mathsf{c}^{\star}(1),\dots,\mathsf{c}^{\star}(n)\in[k]^{n}$, for
each regression problem, where for all $i\in[n],t\in[T]$:

$\theta_{i}\stackrel{{\scriptstyle\textnormal{ind.}}}{{\sim}}\mathsf{N}_{d}(\overline{\theta}_{\mathsf{c}^{\star}(i)},\nu^{2}I_{d}).$

We then observe $\\{y_{i}^{(t)}\\}_{i=1}^{n}$ for $t=1,\dots,T$ according to
the following linear observational model:

$y_{i}^{(t)}=\theta_{i}^{\mathsf{T}}x_{i}^{(t)}+\sum_{\begin{subarray}{c}j\neq
i\\\
\mathsf{c}^{\star}(j)=\mathsf{c}^{\star}(i)\end{subarray}}\gamma_{ij}^{\mathsf{T}}x_{j}^{(t)}+\varepsilon_{i}^{(t)},\quad\mbox{for}~{}i\in[n].$
(1)

For simplicity, we assume that

$\varepsilon_{i}^{(t)}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}\mathsf{N}(0,\sigma^{2})\quad\mbox{and}\quad\gamma_{ij}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}\mathsf{N}_{d}(0,\tau^{2}I_{d}).$

The target $y_{i}^{(t)}$ can be thought of as the future value of $i$th time
series that needs to be predicted. It depends linearly on its own past
$x_{i}^{(t)}$ and also the past of the other time series in the same cluster
$\\{x_{j}^{(t)}\\}$. The local AR parameters $\theta_{i}$ for each cluster
come from the same Gaussian distribution.

###  4.1 Algorithm for Mixed Linear Regression

Our framework for the mixed linear regression problem is instantiated as
Algorithm 4.1.

  * Algorithm 1  _Cluster & Conquer for Mixed Linear Regression._

input: $y_{i}^{(t)},x_{i}^{(t)}$ for $i\in[n],t\in[T].$
---
|  1\. _Train one AR model per time series
$\\{\widehat{\theta}_{i}^{(AR)}\\}_{i=1}^{n}$,_
| _according to Equation ( 2)._
|  2\. _Estimate cluster indices $\\{\widehat{\mathsf{c}}(i)\\}_{i=1}^{n}$
using Algorithm 4.4._
|  3\. _Solve VAR models for each cluster as in ( 5)_.

###  4.2 Overview of Guarantees

We now analyze Algorithm 4.1 in the mixed linear regression setting defined in
Section 4, under the following assumptions.

######  Assumption 1.

The number of clusters $k\geqslant\rho n$ for some constant $\rho\in(0,1)$.

######  Assumption 2.

The cluster sizes are uniform: each cluster has $m=n/k$ time series associated
to it.

######  Assumption 3.

Within each cluster, the design is isotropic. More formally, for cluster
$\iota\in[k]$, suppose that time series $i_{1},\dots,i_{m}\in[n]$ belong to
this cluster (i.e., $\mathsf{c}^{\star}(i_{j})=\iota$ for $j\in[m]$). Define
for each $t\in[T]$

$x^{(t)}_{\iota}=(x_{i_{1}}^{(t)},\dots,x_{i_{m}}^{(t)})\in\mathbf{R}^{md}$

Then we assume $\frac{1}{T}\sum_{t=1}^{\mathsf{T}}x_{\iota}^{(t)}\otimes
x_{\iota}^{(t)}=I_{md},~{}~{}\mbox{for all}~{}\iota\in[k].$

Our analysis will aim to demonstrate that Algorithm 4.1 provides good
estimates of the autoregression parameters. To do so, we begin in Section 4.3
by providing a stochastic decomposition for the scalar AR parameters estimated
for each time series. Having done so, we will establish that the scalar AR
parameters are nearly a Gaussian mixture model. We are then able to extend
existing guarantees for Gaussian mixtures to our setting, giving explicit
conditions under which we can exactly recover the cluster labels for the mixed
linear regression model via spectral clustering on the scalar AR estimates.
These guarantees are provided in Section 4.4. We then study the properties of
a single cluster-wide vector autoregression, and provide a high probability
bound on the errors incurred from model noise. By union bounding over the
clusters and the failure probability in our exact clustering guarantee, we are
able to give a formal statement on the recovery guarantees of Algorithm 4.1.

###  4.3 Scalar Autoregression

As the first step of our procedure we propose fitting scalar autoregression
parameter estimates according to a $\mathrm{AR}(d)$ model. Formally, we
estimate $\widehat{\theta}^{\rm AR}_{i}\in\mathbf{R}^{d}$, for $i\in[n]$ via
the following empirical risk minimization problem,

$\widehat{\theta}^{\rm AR}_{i}\coloneq\mathop{\rm
arg\,min}_{\theta\in\mathbf{R}^{d}}\bigg{\\{}\frac{1}{T}\sum_{t=1}^{T}\big{(}\theta^{\mathsf{T}}x_{i}^{(t)}-y_{i}^{(t)}\big{)}^{2}\bigg{\\}}.$
(2)

Our analysis begins with a stochastic decomposition of $\widehat{\theta}^{\rm
AR}_{i}$. To state it, we introduce the following notation for the (weighted)
empirical covariance matrix, for a given time series. For any weights
$\rho=(\rho^{(1)},\dots,\rho^{(T)})\in\mathbf{R}^{T}$, we define

$\widehat{\Sigma}_{i}(\rho)\coloneq\frac{1}{T}\sum_{t=1}\mathsf{T}(\rho^{(t)})^{2}x_{i}^{(t)}\otimes
x_{i}^{(t)}.$

In the special case $\rho=\mathbf{1}$, we write
$\widehat{\Sigma}_{i}\coloneq\widehat{\Sigma}_{i}(\mathbf{1}).$ With this
notation in hand, we can now state the decomposition mentioned above.

######  Lemma 1.

For each $i\in[n]$, we have $\widehat{\theta}^{\rm
AR}_{i}\stackrel{{\scriptstyle\textnormal{i.i.d.}}}{{\sim}}\mathsf{N}_{d}(\overline{\theta}_{\mathsf{c}^{\star}(i)},\Lambda_{i})$,
where

$\Lambda_{i}\coloneq\nu^{2}I_{d}+\frac{\sigma^{2}}{T}\widehat{\Sigma}_{i}^{-1}+\frac{\tau^{2}}{T}\widehat{\Sigma}_{i}^{-2}\widehat{\Sigma}_{i}(\rho_{i}),$

where $\rho_{i}=(\rho_{i}^{(1)},\dots,\rho_{i}^{(T)})$, with
$\rho_{i}^{(t)}\coloneq\sqrt{\sum_{j:\mathsf{c}^{\star}(j)=\mathsf{c}^{\star}(i)}\|x_{j}^{(t)}\|_{2}^{2}}$,
for $t\in[T]$.

Whereas the true AR parameters $\theta_{i}$ follow a Gaussian mixture model,
Lemma 1 shows that $\widehat{\theta}^{\rm AR}_{i}$, while not a Gaussian
mixture due to the different covariances $\Lambda_{i}$, is at least centered
at the cluster mean $\overline{\theta}_{\mathsf{c}^{\star}(i)}$ for each
$i\in[n]$. The proof has been deferred to the appendix in Section B.

We note that under Assumption 3, the covariances above can be further
simplified to

$\Lambda_{i}\coloneq\Big{(}\nu^{2}+\frac{\sigma^{2}}{T}\Big{)}I_{d}+\frac{\tau^{2}}{T}\widehat{\Sigma}_{i}(\rho_{i}).$
(3)

This will be a useful identity in the sequel.

###  4.4 A Result for Near Mixtures of Gaussians

As we demonstrated in Lemma 1, after computing autoregression estimates
$\widehat{\theta}^{\rm AR}_{i}$, we obtain a near mixture of Gaussians. In
order to obtain guarantees for the vector autoregression, we first show that
the clustering estimates are consistent. We begin with the following general
result for near mixtures of Gaussians.

#### Setting.

To state the result, we begin by formally defining our setting. Suppose there
are $k$ clusters, and $n$ data points $X_{1},\dots,X_{n}\in\mathbf{R}^{d}$.
These data are associated to unknown cluster labels $\mathsf{c}^{\star}(i)$,
where we view $\mathsf{c}^{\star}\colon[n]\to[k]$ as the labelling function.
We assume that there are unknown cluster centers,
$\theta^{\star}_{1},\dots,\theta^{\star}_{k}\in\mathbf{R}^{d}$. Further assume
we observe

$X_{i}\stackrel{{\scriptstyle\textnormal{ind.}}}{{\sim}}\mathsf{N}_{d}(\theta^{\star}_{\mathsf{c}^{\star}(i)},\Sigma_{i}),\quad
i=1,\ldots,n.$

We additionally assume that the sample covariances
$\Sigma_{1},\dots,\Sigma_{n}$ are unknown. Note that this is not quite a
mixture of Gaussians, since we allow that $\Sigma_{i}\neq\Sigma_{j}$, even if
$\mathsf{c}^{\star}(i)=\mathsf{c}^{\star}(j)$.

#### Clustering near Gaussian mixtures.

We study the spectral clustering algorithm, Algorithm 4.4 and provide the
following recovery guarantee for near mixture of Gaussians.

######  Theorem 1.

Suppose that $\Sigma_{i}\preccurlyeq\gamma I_{d}$, for $i\in[n]$. Suppose
Assumption 1 holds. Assume that the centers are well separated:

$\min_{i\neq j\in[k]}\|\theta^{\star}_{i}-\theta^{\star}_{j}\|_{2}\geqslant
32\sqrt{\gamma}k\sqrt{\frac{1+\tfrac{d}{n}}{\beta}}\max\left\\{1,\frac{\beta}{\rho}\right\\}.$

Then, with probability at least $1-\exp(-0.08n)$, there exists a permutation
$\phi\colon[k]\to[k]$ such that $\widehat{\mathsf{c}}_{\rm
spec}(i)=\phi(\mathsf{c}^{\star}(i))$ for all $i\in[n]$.

The proof is presented in Section A of the appendix. The above theorem implies
the following corollary for our mixture of regression problem.

######  Corollary 1.

Suppose that Assumptions 1, 2, and 3 hold. Let

$\lambda\coloneq\nu^{2}+\frac{\sigma^{2}+\tfrac{\tau^{2}}{\rho}}{T}.$

If the following separation condition holds,

$\min_{i\neq j\in[k]}\|\overline{\theta}_{i}-\overline{\theta}_{j}\|_{2}\
\geqslant
32\sqrt{\lambda}k\sqrt{\frac{1+\tfrac{d}{n}}{\beta}}\max\left\\{1,\frac{\beta}{\rho}\right\\},$

then spectral clustering on $\\{\widehat{\theta}^{\rm AR}_{i}\\}_{i=1}^{n}$
produces an estimate $\widehat{\mathsf{c}}_{\rm spec}$ which is exactly equal
to the true clustering $\mathsf{c}^{\star}$ (up to relabellings of
$\\{1,\dots,k\\}$) with probability $1-\exp(-0.08n)$.

  * Algorithm 2 [Spectral clustering] _Algorithm for estimating clusters from samples[23]._

input: sample matrix $X\in\mathbf{R}^{d\times n}$, number of clusters, $k$
---
|  1\. _Compute truncated SVD of sample matrix._
|
$X_{k}\coloneq\sum_{i=1}^{k}\widehat{\sigma}_{i}\widehat{u}_{i}\widehat{v}_{i}^{\mathsf{T}}\eqcolon\widehat{U}\widehat{\Sigma}\widehat{V}^{\mathsf{T}}\in\mathbf{R}^{d\times
n}$
|  2\. _Compute projection onto top- $k$ left singular vectors._
|
$Y\coloneq\widehat{U}^{\mathsf{T}}X_{k}=\widehat{\Sigma}\widehat{V}^{\mathsf{T}}\in\mathbf{R}^{k\times
n}$
|  3\. _Solve $k$-means on the columns of $Y$._ Solve:
|
$\mbox{minimize}\quad\frac{1}{2n}\sum_{i=1}^{n}\|Y_{i}-\eta_{\kappa_{i}}\|_{2}^{2}$,
|  with variables $\eta_{1},\dots,\eta_{k}\in\mathbf{R}^{k}$ and
$\kappa_{1},\dots,\kappa_{n}\in[k]$.
output: labelling $\widehat{\mathsf{c}}_{\rm spec}\colon[n]\to[k]$, where
$\widehat{\mathsf{c}}_{\rm spec}(i)\coloneq\kappa_{i}^{\star}$.

###  4.5 Vector Autoregression per Cluster

The final step of our argument requires us to compute a high-probability
guarantee on the closeness of estimated vector autoregression parameters.

For each cluster $\iota\in[k]$, under Assumption 2, there is an index set,

$\\{j\in[n]:\mathsf{c}^{\star}(j)=\iota\\}=\Big{\\{}i_{1},\dots,i_{m}\Big{\\}}\subset[n],$

that specifies the indices of the time series that truly belong to cluster
$\iota\in[k]$. Define the following vectorized data:

$y_{\iota}^{(t)}=\begin{bmatrix}y_{i_{1}}^{(t)}\\\ \vdots\\\
y_{i_{m}}^{(t)}\end{bmatrix}\quad\mbox{and}\quad
x_{\iota}^{(t)}=\begin{bmatrix}x_{i_{1}}^{(t)}\\\ \vdots\\\
x_{i_{m}}^{(t)}\end{bmatrix}.$

Hence, $y_{\iota}^{(t)}\in\mathbf{R}^{m}$ and
$x_{\iota}^{(t)}\in\mathbf{R}^{md}$, for $\iota\in[k]$ and $t\in[T]$.

With this notation in hand, for each cluster $\iota\in[k]$, the observational
model stated in (1) is equivalent to

$y_{\iota}^{(t)}=(I_{m}\otimes(x_{\iota}^{(t)})^{\mathsf{T}})\Gamma_{\iota}^{\star}+\varepsilon_{\iota}^{(t)},\mbox{for
all}~{}t\in[T],~{}\iota\in[k].$

Here we define $\Gamma_{\iota}^{\star}\in\mathbf{R}^{m^{2}d}$ to be the
concatenation of the parameters in cluster $\iota\in[k]$. Specifically, in
view of our observational model (1), suppose that we define

$\widetilde{\gamma}_{i_{u}i_{v}}\coloneq\begin{cases}\theta_{i_{u}}&u=v\\\
\gamma_{i_{u}i_{v}}&u\neq v\end{cases},\quad\mbox{for}~{}u,v\in[m].$

Then

$\Gamma_{\iota}^{\star}=\Big{(}\widetilde{\gamma}_{i_{1}i_{1}},\widetilde{\gamma}_{i_{1},i_{2}},\ldots,\widetilde{\gamma}_{i_{2}i_{1}},\widetilde{\gamma}_{i_{2},i_{2}},\ldots,\gamma_{i_{m}i_{m}}\Big{)}\in\mathbf{R}^{m^{2}d}.$

This is our target of estimation.

#### Oracle VAR estimator.

The first estimator we study is what we refer to as the _oracle VAR estimate
for cluster $\iota$_:

$\Gamma^{\rm
VAR,or}_{\iota}\coloneq\underset{\Gamma\in\mathbf{R}^{m^{2}d}}{\mathop{\rm
arg\,min}}\bigg{\\{}\frac{1}{2T}\sum_{t=1}^{T}\|y_{\iota}^{(t)}-(I_{m}\otimes(x_{\iota}^{(t)})^{\mathsf{T}})\Gamma\|_{2}^{2}\bigg{\\}}.$

(4)

Crucially, it assumes knowledge of the set $\\{i_{1},\dots,i_{m}\\}$, which
are those time series which truly belong to cluster $\iota\in[k]$.

For this estimator, we have the following guarantee.

######  Lemma 2.

Suppose that Assumption 3 holds. Then for each $\iota\in[k]$, with probability
at least $1-\delta$ we have

$\big{\|}\Gamma^{\rm
VAR,or}_{\iota}-\Gamma^{\star}_{\iota}\big{\|}_{2}\leqslant\sqrt{2}\sqrt{\frac{\sigma^{2}m^{2}d}{T}}+\sqrt{3\frac{\sigma^{2}}{T}\log\frac{1}{\delta}}.$

Deferring the proof to the appendix, an outline is as follows. First we
demonstrate that under the stated assumptions, $\Gamma^{\rm
VAR,or}_{i}\sim\mathsf{N}(\Gamma^{\star},\tfrac{\sigma^{2}}{T}I_{m^{2}d})$.
Consequently,

$\operatorname*{\mathbf{E}}\ilimits@\Big{[}\|\Gamma^{\rm
VAR,or}_{i}-\Gamma_{\iota}^{\star}\|_{2}^{2}\Big{]}\leqslant\frac{\sigma^{2}m^{2}d}{T}.$

Establishing concentration around the mean is a consequence of tail-bounds for
$\chi^{2}$ random variates. We provide the full proof in Section B of the
appendix.

#### VAR estimator with estimated clustering.

We will show that the preceding result and Corollary 1, immediately imply our
main recovery guarantee, via a simple appeal to a union bound. To state this
result, we first need to describe how to construct the VAR estimate when the
ground truth clustering is unknown. To do this, we define the following index
set:

$\\{j\in[n]:\widehat{\mathsf{c}}(j)=\iota\\}=\\{\widehat{i}_{1},\dots,\widehat{i}_{m_{\iota}}\\}\subset[n].$

This is the set of time series which are estimated to belong to cluster
$\iota$. Define

$\widehat{y}_{\iota}^{(t)}=\begin{bmatrix}y_{\widehat{i}_{1}}^{(t)}\\\
\vdots\\\
y_{\widehat{i}_{m_{\iota}}}^{(t)}\end{bmatrix}\quad\mbox{and}\quad\widehat{x}_{\iota}^{(t)}=\begin{bmatrix}x_{\widehat{i}_{1}}^{(t)}\\\
\vdots\\\ x_{\widehat{i}_{m_{\iota}}}^{(t)}\end{bmatrix}.$

Hence, $\widehat{y}_{\iota}^{(t)}\in\mathbf{R}^{m_{\iota}}$ and
$\widehat{x}_{\iota}^{(t)}\in\mathbf{R}^{m_{\iota}d}$, for $\iota\in[k]$ and
$t\in[T]$. We now form our _VAR estimate for cluster $\iota$_:

$\Gamma^{\rm
VAR}_{\iota}\coloneq\underset{\Gamma\in\mathbf{R}^{m_{\iota}^{2}d}}{\mathop{\rm
arg\,min}}\bigg{\\{}\frac{1}{2T}\sum_{t=1}^{T}\|\widehat{y}_{\iota}^{(t)}-(I_{m}\otimes(\widehat{x}_{\iota}^{(t)})^{\mathsf{T}})\Gamma\|_{2}^{2}\bigg{\\}}.$

(5)

######  Theorem 2 (Recovery guarantee for mixed linear regression).

Suppose that Assumptions 1, 2, and 3 hold for the mixed linear regression
problem. Additionally assume that the separation conditions of Corollary 1
hold. Then, with probability at least $1-\delta-\exp(-0.08n)$, we have the
following guarantee:

$\max_{\iota\in[k]}\big{\|}\Gamma^{\rm
VAR}_{\iota}-\Gamma^{\star}_{\iota}\big{\|}_{2}\leqslant\sqrt{2}\sqrt{\frac{\sigma^{2}m^{2}d}{T}}+\sqrt{3\frac{\sigma^{2}}{T}\log\frac{k}{\delta}}.$

This result is interesting because it indicates that, up to logarithmic and
constant factors, the worst-case recovery error of $\Gamma^{\star}$ is
$\sqrt{\tfrac{\sigma^{2}m^{2}d}{T}}$, which is the standard rate for a
$(m^{2}d)$-dimensional parameter, with $T$ samples and noise variance
$\sigma^{2}$.

##  5 Empirical Results

\begin{overpic}[width=433.62pt]{figs/vs_clusters_k10.pdf} \end{overpic}

\begin{overpic}[width=433.62pt]{figs/vs_clusters_k20.pdf} \end{overpic}

\begin{overpic}[width=433.62pt]{figs/vs_clusters_k30.pdf} \end{overpic}

\begin{overpic}[width=433.62pt]{figs/vs_clusters_k10_p20.pdf} \end{overpic}

\begin{overpic}[width=433.62pt]{figs/vs_clusters_k20_p20.pdf} \end{overpic}

\begin{overpic}[width=433.62pt]{figs/vs_clusters_k30_p20.pdf} \end{overpic}

Figure 2: Here we demonstrate the performance of Cluster & Conquer on
simulated data with different configurations. We evaluate our method as a
function of the number of clusters and compare it to a Scalar AR model and an
ablation model that uses random clusters. We see that Cluster & Conquer is
always better than Scalar AR and the ablation model.
\begin{overpic}[width=186.45341pt]{figs/time_clusters.pdf} \end{overpic}
Figure 3: Training time of VAR models on the whole dataset decreases with an
increase in number of clusters. When a cluster is too big, VAR cannot scale as
the computation complexity increases super-quadratically w.r.t number of time
series in the VAR model.

In this section we demonstrate the performance of our approach on both
simulated and real-world datasets. Specifically, we show that our framework,
in combination with a simple linear function class (AR/VAR), is able to
outperform several recently proposed deep learning methods for time-series
forecasting.

#### Setup.

Given multivariate time series, we first use LinearSVR [18] for training local
AR models. Then we use the learned AR coefficients to cluster the time series.
To do so, we construct a $K$-nearest-neighbor graph ($K$=11) and then use the
graph-clustering algorithm METIS [24] to form partitions. METIS is extremely
scalable and has a sub-quadratic complexity in terms of the number of data
points ($n$), unlike direct $k$-means or spectral clustering. Finally, we
train global VAR models for each cluster using LinearSVR. We provide details
about our implementation, including hyper-parameter settings, in the appendix
in Section C. We refer to our algorithm as Cluster & Conquer.

#### Experiments on simulated data.

We perform a simulation study to empirically verify that our algorithm works
for non-i.i.d time series data given our generative assumptions. To do so, we
consider several different configurations, i.e., we vary the number of
clusters and the lag and period of the ground truth time series data. In each
case we construct dataset with $n=200$ time series. Details about the data
generation process are provided in the appendix in Section C.

Figure 2 shows results for the different configurations, along with 1-standard
deviation bands over 20 runs. The top row shows results for data with period
and lag of 20, whereas the bottom row shows results for data that have period
and lag of 10. The three columns correspond to ground truth datasets that have
10, 20 and 30 clusters, respectively, from left to right.

To evaluate our model, we consider a rolling validation task with forecast
horizon $24$ and with $7$ rolling windows. In each of the subplots, we plot
the performance as measured by WAPE (see (15) in the appendix) as a function
of the number of clusters ($k$) in our algorithm. We see that the performance
of the algorithm is optimal or near-optimal at the correct value of $k$ (i.e.,
when $k$ is equal to number of clusters in the ground truth data). We also
show that our method performs better compared to a Scalar AR model. Moreover,
the performance is also better compared to an ablation model that considers
random clusters followed by VAR within the clusters. We also note that the
algorithm is robust to $k$ provided $k$ is not too far away from the number of
clusters in the ground truth data.

In Figure 3, we plot the running time of the algorithm as $k$ increases. As
suspected when $k$ is small the running time is prohibitive as the VAR model
in each cluster has to learn $O((n/k)^{2})$ parameters. The shaded area
corresponds to the standard deviation across 20 runs.

| MAE | RMSE | MAPE
---|---|---|---
FC-LSTM [37] | 3.57 | 6.2 | 8.6
SFM [44] | 2.75 | 4.32 | 6.6
N-BEATS [33] | 3.41 | 5.52 | 7.65
DCRNN [27] | 2.25 | 4.04 | 5.30
LSTNet [26] | 2.34 | 4.26 | 5.41
ST-GCN [41] | 2.25 | 4.04 | 5.26
TCN [4] | 3.25 | 5.51 | 6.7
DeepState [34] | 3.95 | 6.49 | 7.9
DeepGLO [36] | 3.01 | 5.25 | 6.2
StemGNN [9] | 2.14 | 4.01 | 5.01
Scalar AR | 2.13$\pm$0.080 | 3.95$\pm$0.157 | 4.88$\pm$0.2
Cluster & Conquer |  1.99$\pm$0.030 |  3.68$\pm$0.069 |  4.52$\pm$0.1

Table 1: Rolling validation results for the PEMS07 dataset.

| MAE | RMSE | MAPE
---|---|---|---
FC-LSTM [37] | 0.32 | 0.54 | 31.0
SFM [44] | 0.17 | 0.58 | 11.9
N-BEATS [33] | 0.08 | 0.16 | 12.428
LSTNet [26] | 0.08 | 0.12 | 12.74
TCN [4] | 0.1 | 0.3 | 19.03
DeepState [34] | 0.09 | 0.76 | 19.21
GraphWaveNet [40] | 0.19 | 0.86 | 19.67
DeepGLO [36] | 0.09 | 0.15 | 12.45
StemGNN [9] | 0.05 | 0.07 | 10.58
Scalar AR | 0.024$\pm$0.000 | 0.041$\pm$0.000 | 5.67$\pm$0.1
Cluster & Conquer |  0.023$\pm$0.000 |  0.040$\pm$0.001 |  5.39$\pm$0.0

Table 2: Rolling validation results for the ECG dataset.

| WAPE | MAPE | SMAPE
---|---|---|---
DeepGLO [36] | 8.2 | 34.1 | 12.1
TCN (LeveledInit) [36] | 9.2 | 23.7 | 12.6
LSTM | 10.9 | 26.4 | 15.4
DeepAR [35] | 8.6 | 25.9 | 14.1
TCN [8] | 14.7 | 47.6 | 15.6
Prophet [16] | 19.7 | 39.3 | 22.1
TRMF [43] | 10.4 | 28.0 | 15.1
SVD+TCN [36] | 21.9 | 43.7 | 23.8
Scalar AR | 8.11$\pm$0.14 | 24.7$\pm$0.65 | 11.2$\pm$0.19
Cluster & Conquer |  7.97$\pm$0.09 |  20.48$\pm$0.58 |  10.8$\pm$0.08

Table 3: Rolling val. results for the Electricity dataset.

| WAPE | MAPE | SMAPE
---|---|---|---
DeepGLO [36] | 14.8 | 16.8 | 14.2
TCN (LeveledInit) [36] | 15.7 | 20.1 | 15.6
LSTM | 27.0 | 35.7 | 26.3
DeepAR [35] | 14.0 | 20.1 | 11.4
TCN [8] | 20.4 | 28.4 | 23.6
Prophet [16] | 30.3 | 55.9 | 42.0
TRMF [43] | 15.9 | 22.6 | 18.1
SVD+TCN [36] | 32.9 | 68.7 | 34.0
Scalar AR | 15.9$\pm$0.31 | 20.0$\pm$1.14 | 15.72$\pm$0.74
Cluster & Conquer |  14.0$\pm$0.16 |  17.18$\pm$0.49 | 13.68$\pm$0.33

Table 4: Rolling val. results for the Traffic dataset.

#### Experiments on real data.

Here, we follow the setups described by Cao et al. [9] and Sen et al. [36]. We
refer the reader to the appendix, particularly Section C, for details about
the baselines and our implementation, as well as statistics of the datasets.
For all datasets, we choose the number of clusters in Cluster & Conquer to be
$\left\lfloor n/10\right\rfloor$, where $n$ is the number of time series. The
number of clusters was manually tuned on the validation set between
$\\{\left\lfloor n/5\right\rfloor,\left\lfloor n/10\right\rfloor,\left\lfloor
n/20\right\rfloor\\}$ and in all cases we found $\left\lfloor
n/10\right\rfloor$ to be the best.

First we compare the Cluster & Conquer approach to baselines reported by [9].
We show results for the PEMS07 [10] dataset in Table 1, and the ECG dataset
[11] in Table 2. The Cluster & Conquer approach outperforms all other models
on all metrics, including state-of-the-art deep learning models. This is
remarkable, because we only instantiate our framework with linear models.

Next, we compare Cluster & Conquer against baselines reported by [36] for the
Electricity and Traffic datasets, using the same metrics used in [36] for a
fair comparison. Table 3 shows results for the Electricity dataset. Here, the
advantage of the Cluster & Conquer approach is pronounced, i.e., the gap
between our approach and other deep learning models is distinct on all three
metrics. On the other hand, the performance gap is less distinct for the
Traffic dataset in Table 4.

Given these results, a few remarks are in place. First, it should be noted
that the Cluster & Conquer approach always improves over the scalar AR model.
This shows that capturing intra-cluster dependencies is important for
forecasting. Second, the Cluster & Conquer approach is competitive to global
deep learning models (e.g. DeepGLO, DCRNN, LSTNet). One reason for that is
that the Cluster & Conquer approach suppresses noise from irrelevant time
series when predicting the future of a particular time series. Third, the
Cluster & Conquer approach does not use specially designed covariates. For
instance, DeepAR [35] uses specially designed time-based covariates such as
transformations of days of the week, and holidays.

##  6 Conclusions

We proposed a three-stage framework, based on a divide and conquer principle,
for large-scale multivariate time-series forecasting. Our Cluster & Conquer
framework is highly flexible, i.e., it allows for the use of any multivariate
time-series model coupled with a scalable clustering algorithm. In turn, this
framework is computationally efficient, embarrassingly parallel, and typically
more efficient than dense vector auto regression. We provide theoretical
guarantees for our framework in a mixed linear regression setting. Experiments
on both synthetic and multiple real-world datasets show that the Cluster &
Conquer approach achieves state-of-the-art results compared to baselines.

## References

  * Aghabozorgi et al. [2015] S. Aghabozorgi, A. Seyed Shirkhorshidi, and T. Ying Wah.  Time-series clustering - a decade review.  _Inf. Syst._ , 53(C):16–38, Oct. 2015.  ISSN 0306-4379.  doi: 10.1016/j.is.2015.04.007.  URL https://doi.org/10.1016/j.is.2015.04.007.
  * Azencot et al. [2020] O. Azencot, N. B. Erichson, V. Lin, and M. Mahoney.  Forecasting sequential data using consistent koopman autoencoders.  In _International Conference on Machine Learning_ , pages 475–485. PMLR, 2020.
  * Bai et al. [2020] L. Bai, L. Yao, C. Li, X. Wang, and C. Wang.  Adaptive graph convolutional recurrent network for traffic forecasting.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Bai et al. [2018] S. Bai, J. Z. Kolter, and V. Koltun.  An empirical evaluation of generic convolutional and recurrent networks for sequence modeling.  _CoRR_ , abs/1803.01271, 2018.  URL http://arxiv.org/abs/1803.01271.
  * Bandara et al. [2020] K. Bandara, C. Bergmeir, and S. Smyl.  Forecasting across time series databases using recurrent neural networks on groups of similar series: A clustering approach.  _Expert Syst. Appl._ , 140, 2020.  doi: 10.1016/j.eswa.2019.112896.  URL https://doi.org/10.1016/j.eswa.2019.112896.
  * Benidis et al. [2020] K. Benidis, S. S. Rangapuram, V. Flunkert, B. Wang, D. Maddix, C. Turkmen, J. Gasthaus, M. Bohlke-Schneider, D. Salinas, L. Stella, et al.  Neural forecasting: Introduction and literature overview.  _arXiv preprint arXiv:2004.10240_ , 2020.
  * Bhatia et al. [2010] N. Bhatia et al.  Survey of nearest neighbor techniques.  _arXiv preprint arXiv:1007.0085_ , 2010.
  * Borovykh et al. [2017] A. Borovykh, S. Bohte, and C. W. Oosterlee.  Conditional time series forecasting with convolutional neural networks.  _arXiv preprint arXiv:1703.04691_ , 2017.
  * Cao et al. [2020] D. Cao, Y. Wang, J. Duan, C. Zhang, X. Zhu, C. Huang, Y. Tong, B. Xu, J. Bai, J. Tong, and Q. Zhang.  Spectral temporal graph neural network for multivariate time-series forecasting.  In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/cdf6581cb7aca4b7e19ef136c6e601a5-Abstract.html.
  * Chen et al. [2001] C. Chen, K. Petty, A. Skabardonis, P. Varaiya, and Z. Jia.  Freeway performance measurement system: Mining loop detector data.  _Transportation Research Record_ , 1748(1):96–102, 2001.  doi: 10.3141/1748-12.  URL https://doi.org/10.3141/1748-12.
  * Dau et al. [2019] H. A. Dau, A. J. Bagnall, K. Kamgar, C. M. Yeh, Y. Zhu, S. Gharghabi, C. A. Ratanamahatana, and E. J. Keogh.  The UCR time series archive.  _IEEE CAA J. Autom. Sinica_ , 6(6):1293–1305, 2019.  doi: 10.1109/jas.2019.1911747.  URL https://doi.org/10.1109/jas.2019.1911747.
  * Ding et al. [2015] R. Ding, Q. Wang, Y. Dang, Q. Fu, H. Zhang, and D. Zhang.  Yading: fast clustering of large-scale time series data.  _Proceedings of the VLDB Endowment_ , 8(5):473–484, 2015.
  * Divine [2012] D. V. Divine.  M. mudelsee, , climate time series analysis: Classical statistical and bootstrap methods (2010) springer, dordrecht 978-90-481-9482-7.  _Comput. Geosci._ , 43:24, 2012.  doi: 10.1016/j.cageo.2012.02.030.  URL https://doi.org/10.1016/j.cageo.2012.02.030.
  * Erichson et al. [2019] N. B. Erichson, M. Muehlebach, and M. W. Mahoney.  Physics-informed autoencoders for lyapunov-stable fluid flow prediction.  _arXiv preprint arXiv:1905.10866_ , 2019.
  * Erichson et al. [2021] N. B. Erichson, O. Azencot, A. Queiruga, L. Hodgkinson, and M. W. Mahoney.  Lipschitz recurrent neural networks.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=-N7PBXqOUJZ.
  * Facebook [2020] Facebook.  Fb prophet, forecasting at scale.  https://github.com/facebook/prophet, 2020.
  * Faloutsos et al. [2020] C. Faloutsos, V. Flunkert, J. Gasthaus, T. Januschowski, and Y. Wang.  Forecasting big time series: Theory and practice.  In A. E. F. Seghrouchni, G. Sukthankar, T. Liu, and M. van Steen, editors, _Companion of The 2020 Web Conference 2020, Taipei, Taiwan, April 20-24, 2020_ , pages 320–321. ACM / IW3C2, 2020.  doi: 10.1145/3366424.3383118.  URL https://doi.org/10.1145/3366424.3383118.
  * Fan et al. [2008] R. Fan, K. Chang, C. Hsieh, X. Wang, and C. Lin.  LIBLINEAR: A library for large linear classification.  _J. Mach. Learn. Res._ , 9:1871–1874, 2008.  URL https://dl.acm.org/citation.cfm?id=1442794.
  * Goldenshluger and Zeevi [2001] A. Goldenshluger and A. Zeevi.  Nonasymptotic bounds for autoregressive time series modeling.  _Ann. Statist._ , 29(2):417–444, 2001.  ISSN 0090-5364.  doi: 10.1214/aos/1009210547.  URL https://doi.org/10.1214/aos/1009210547.
  * Guo et al. [2016] S. Guo, Y. Wang, and Q. Yao.  High-dimensional and banded vector autoregressions.  _Biometrika_ , 103(4):889–903, 2016.  ISSN 0006-3444.  doi: 10.1093/biomet/asw046.  URL https://doi.org/10.1093/biomet/asw046.
  * Hyndman and Athanasopoulos [2018] R. J. Hyndman and G. Athanasopoulos.  _Forecasting: principles and practice_.  OTexts, 2018.
  * Kalpakis et al. [2001] K. Kalpakis, D. Gada, and V. Puttagunta.  Distance measures for effective clustering of ARIMA time-series.  In N. Cercone, T. Y. Lin, and X. Wu, editors, _Proceedings of the 2001 IEEE International Conference on Data Mining, 29 November - 2 December 2001, San Jose, California, USA_ , pages 273–280. IEEE Computer Society, 2001.  doi: 10.1109/ICDM.2001.989529.  URL https://doi.org/10.1109/ICDM.2001.989529.
  * Kannan and Vempala [2008] R. Kannan and S. Vempala.  Spectral algorithms.  _Found. Trends Theor. Comput. Sci._ , 4(3-4):front matter, 157–288 (2009), 2008.  ISSN 1551-305X.  doi: 10.1561/0400000025.  URL https://doi.org/10.1561/0400000025.
  * Karypis and Kumar [1998] G. Karypis and V. Kumar.  _MeTiS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices_.  Department of Computer Science, University of Minnesota, version 4.0 edition, Sept. 1998.
  * Kuznetsov and Mohri [2015] V. Kuznetsov and M. Mohri.  Learning theory and algorithms for forecasting non-stationary time series.  In _NIPS_ , pages 541–549. Citeseer, 2015.
  * Lai et al. [2018] G. Lai, W. Chang, Y. Yang, and H. Liu.  Modeling long- and short-term temporal patterns with deep neural networks.  In K. Collins-Thompson, Q. Mei, B. D. Davison, Y. Liu, and E. Yilmaz, editors, _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval, SIGIR 2018, Ann Arbor, MI, USA, July 08-12, 2018_, pages 95–104. ACM, 2018.  doi: 10.1145/3209978.3210006.  URL https://doi.org/10.1145/3209978.3210006.
  * Li et al. [2017] Y. Li, R. Yu, C. Shahabi, and Y. Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  _arXiv preprint arXiv:1707.01926_ , 2017.
  * Li et al. [2018] Y. Li, R. Yu, C. Shahabi, and Y. Liu.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  In _International Conference on Learning Representations_ , 2018.
  * Lim et al. [2021] S. H. Lim, N. B. Erichson, L. Hodgkinson, and M. W. Mahoney.  Noisy recurrent neural networks.  _arXiv preprint arXiv:2102.04877_ , 2021.
  * Löffler et al. [2019] M. Löffler, A. Y. Zhang, and H. H. Zhou.  Optimality of spectral clustering for gaussian mixture model.  _CoRR_ , abs/1911.00538, 2019.  URL http://arxiv.org/abs/1911.00538.
  * Matsubara et al. [2014] Y. Matsubara, Y. Sakurai, W. G. van Panhuis, and C. Faloutsos.  FUNNEL: automatic mining of spatially coevolving epidemics.  In S. A. Macskassy, C. Perlich, J. Leskovec, W. Wang, and R. Ghani, editors, _The 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’14, New York, NY, USA - August 24 - 27, 2014_ , pages 105–114. ACM, 2014.  doi: 10.1145/2623330.2623624.  URL https://doi.org/10.1145/2623330.2623624.
  * McDonald et al. [2017] D. J. McDonald, C. R. Shalizi, and M. Schervish.  Nonparametric risk bounds for time-series forecasting.  _J. Mach. Learn. Res._ , 18:Paper No. 32, 40, 2017.  ISSN 1532-4435.
  * Oreshkin et al. [2020] B. N. Oreshkin, D. Carpov, N. Chapados, and Y. Bengio.  N-BEATS: neural basis expansion analysis for interpretable time series forecasting.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.  URL https://openreview.net/forum?id=r1ecqn4YwB.
  * Rangapuram et al. [2018] S. S. Rangapuram, M. W. Seeger, J. Gasthaus, L. Stella, Y. Wang, and T. Januschowski.  Deep state space models for time series forecasting.  In S. Bengio, H. M. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pages 7796–7805, 2018.  URL https://proceedings.neurips.cc/paper/2018/hash/5cf68969fb67aa6082363a6d4e6468e2-Abstract.html.
  * Salinas et al. [2019] D. Salinas, V. Flunkert, J. Gasthaus, and T. Januschowski.  Deepar: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ , 2019.
  * Sen et al. [2019] R. Sen, H. Yu, and I. S. Dhillon.  Think globally, act locally: A deep neural network approach to high-dimensional time series forecasting.  In H. M. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché-Buc, E. B. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 4838–4847, 2019.  URL https://proceedings.neurips.cc/paper/2019/hash/3a0844cee4fcf57de0c71e9ad3035478-Abstract.html.
  * Sutskever et al. [2014] I. Sutskever, O. Vinyals, and Q. V. Le.  Sequence to sequence learning with neural networks.  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada_ , pages 3104–3112, 2014.  URL https://proceedings.neurips.cc/paper/2014/hash/a14ac55a4f27472c5d894ec1c3c743d2-Abstract.html.
  * Venkataramana Kini and Chandra Sekhar [2013] B. Venkataramana Kini and C. Chandra Sekhar.  Bayesian mixture of AR models for time series clustering.  _PAA Pattern Anal. Appl._ , 16(2):179–200, 2013\.  ISSN 1433-7541.  doi: 10.1007/s10044-011-0247-5.  URL https://doi.org/10.1007/s10044-011-0247-5.
  * Warren Liao [2005] T. Warren Liao.  Clustering of time series data-a survey.  _Pattern Recogn._ , 38(11):1857–1874, Nov. 2005\.  ISSN 0031-3203.  doi: 10.1016/j.patcog.2005.01.025.  URL https://doi.org/10.1016/j.patcog.2005.01.025.
  * Wu et al. [2019] Z. Wu, S. Pan, G. Long, J. Jiang, and C. Zhang.  Graph wavenet for deep spatial-temporal graph modeling.  In _IJCAI_ , 2019.
  * Yu et al. [2017] B. Yu, H. Yin, and Z. Zhu.  Spatio-temporal graph convolutional networks: A deep learning framework for traffic forecasting.  _arXiv preprint arXiv:1709.04875_ , 2017.
  * Yu et al. [2018] B. Yu, H. Yin, and Z. Zhu.  Spatio-temporal graph convolutional networks: a deep learning framework for traffic forecasting.  In _Proceedings of the 27th International Joint Conference on Artificial Intelligence_ , pages 3634–3640, 2018.
  * Yu et al. [2016] H.-F. Yu, N. Rao, and I. S. Dhillon.  Temporal regularized matrix factorization for high-dimensional time series prediction.  In _Advances in neural information processing systems_ , pages 847–855, 2016.
  * Zhang et al. [2017] L. Zhang, C. C. Aggarwal, and G. Qi.  Stock price prediction via discovering multi-frequency trading patterns.  In _Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Halifax, NS, Canada, August 13 - 17, 2017_ , pages 2141–2149. ACM, 2017.  doi: 10.1145/3097983.3098117.  URL https://doi.org/10.1145/3097983.3098117.
  * Zhu and Shasha [2002] Y. Zhu and D. E. Shasha.  Statstream: Statistical monitoring of thousands of data streams in real time.  In _Proceedings of 28th International Conference on Very Large Data Bases, VLDB 2002, Hong Kong, August 20-23, 2002_ , pages 358–369. Morgan Kaufmann, 2002.  doi: 10.1016/B978-155860869-6/50039-1.  URL http://www.vldb.org/conf/2002/S10P04.pdf.

##  Appendix A Proof of result for near mixture of Gaussians

#### Notation:

Let us define the following loss function,
$\ell\colon[k]^{[n]}\times[k]^{[n]}\to\mathbf{R}_{+}$, given by

$\ell(\mathsf{c},\mathsf{c}^{\prime})\coloneq\min_{\varphi}\sum_{i=1}^{n}\mathbf{1}\Big{\\{}\mathsf{c}(i)\neq\varphi(\mathsf{c}^{\prime}(i))\Big{\\}},\quad\mbox{for
any}~{}\mathsf{c},\mathsf{c}^{\prime}\in[k]^{[n]},$

where the minimum ranges over bijections $\varphi\colon[k]\to[k]$.
Intuitively, this measures the number disagreements between the cluster
labellings $\mathsf{c},\mathsf{c}^{\prime}$. Let us set

$\Delta\coloneq\min_{i\neq
j\in[k]}\|\theta^{\star}_{i}-\theta^{\star}_{j}\|_{2}$

#### Proof:

Since Algorithms 4.4 and A.1 differ only by the application of the orthogonal
matrix $\widehat{U}$, it is striaghtforward to see that there is a bijection
$\phi\colon[k]\to[k]$ such that $\widehat{\mathsf{c}}_{\rm
appx}(i)=\phi(\widehat{\mathsf{c}}_{\rm spec}(i))$ for each $i\in[n]$.
Therefore, by definition of $\ell(\cdot,\cdot)$, we have that

$\ell(\widehat{\mathsf{c}}_{\rm
spec},\mathsf{c}^{\star})=\ell(\widehat{\mathsf{c}}_{\rm
appx},\mathsf{c}^{\star}).$

Suppose that the following good event holds:

$\mathcal{E}\coloneq\Big{\\{}\left\|G\right\|_{\rm
op}\leqslant\sqrt{2\gamma}(\sqrt{n}+\sqrt{d})\Big{\\}}.$

Consider the sets

$\mathcal{I}_{\rm
bad}\coloneq\Big{\\{}i\in[n]:\|\theta^{\star}_{\mathsf{c}^{\star}(i)}-\widehat{\theta}_{i}\|_{2}\geqslant\Delta/2\Big{\\}}\qquad\mbox{and}\qquad\mathcal{I}_{\rm
good}\coloneq[n]\setminus\mathcal{I}_{\rm bad}.$

Under $\mathcal{E}$, it is possible to show that

$|\mathcal{I}_{\rm bad}|\leqslant\frac{128k\left\|G\right\|_{\rm
op}^{2}}{\Delta^{2}}\leqslant\frac{512\gamma k(n+d)}{\Delta^{2}}.$

(The first inequality above is Proposition 1.) We can partition the good
indices into disjoint sets $I_{j}$ as follows:

$I_{j}\coloneq\Big{\\{}i\in[n]:\mathsf{c}^{\star}(i)=j,~{}\|\widehat{\theta}_{i}-\theta^{\star}_{\mathsf{c}^{\star}(i)}\|_{2}<\Delta/2\Big{\\}},\qquad
j=1,\ldots,k.$

Note that under the assumption on $\Delta$, $|\mathcal{I}_{\rm
bad}|\leqslant\tfrac{\beta n}{2k}$, which implies that $I_{j}$ are nonempty
for all $j\in[k]$. Additionally, $\widehat{\mathsf{c}}_{\rm appx}(I_{j})$ are
pairwise disjoint, as we demonstrate in Lemma 4, after which it follows that
there exists a permutation $\varphi$ such that

$\varphi(\widehat{\mathsf{c}}_{\rm appx}(i))=j,\quad\mbox{for each}~{}i\in
I_{j}.$

This implies that under $\mathcal{E}$

$\ell(\widehat{\mathsf{c}}_{\rm
appx},\mathsf{c}^{\star})\leqslant|\mathcal{I}_{\rm
bad}|\leqslant\frac{512\gamma k(n+d)}{\Delta^{2}}\leqslant\frac{1}{2},$

which implies that $\widehat{\mathsf{c}}_{\rm appx},\mathsf{c}^{\star}$ are
equal up to permutation. The result now follows by noting that
$\operatorname*{\mathbf{P}}\ilimits@(\mathcal{E})\geqslant
1-\mathrm{e}^{-0.08n}$, as we show in Lemma 3.

###  A.1 Auxillary results for Theorem 1

We assume that for all $i$, we have $\Sigma_{i}\preccurlyeq\gamma I_{d}$ for
some constant $\gamma>0$. It will be convenient for us to consider the sample
matrix $X\in\mathbf{R}^{d\times n}$, which has columns
$X_{1},\dots,X_{n}\in\mathbf{R}^{d}$. Then it is easy to see that

$X=\operatorname*{\mathbf{E}}\ilimits@X+(X-\operatorname*{\mathbf{E}}\ilimits@X)\eqcolon\Theta^{\star}+G,$

where $\Theta^{\star},G\in\mathbf{R}^{d\times n}$. In the decomposition above,
$\Theta^{\star}$ has columns
$\theta^{\star}_{\mathsf{c}^{\star}(1)},\dots,\theta^{\star}_{\mathsf{c}^{\star}(n)}$,
whereas $G$ has columns $g_{1},\dots,g_{n}\in\mathbf{R}^{d}$, where

$g_{i}\stackrel{{\scriptstyle\textnormal{ind.}}}{{\sim}}\mathsf{N}(0,\Sigma_{i}),\quad
i=1,\dots,n.$

######  Lemma 3.

Let $Z\in\mathbf{R}^{d\times n}$ be a random matrix with columns
$z_{i}\stackrel{{\scriptstyle\textnormal{ind.}}}{{\sim}}\mathsf{N}(0,\Sigma_{i})$,
for $i=1,\dots,n$. Let $\lambda_{\rm
max}\coloneq\max_{i=1,\dots,n}\left\|\Sigma_{i}^{1/2}\right\|_{\rm op}$. Then

$\operatorname*{\mathbf{P}}\ilimits@\Big{\\{}\left\|Z\right\|_{\rm
op}\geqslant\lambda_{\rm
max}\big{(}\sqrt{n}+\sqrt{d}+t\big{)}\Big{\\}}\leqslant\exp\Big{(}-\frac{t^{2}}{2}\Big{)}.$

###### Proof.

We note that $\left\|Z\right\|_{\rm op}\leqslant\lambda_{\rm
max}\left\|W\right\|_{\rm op}$ almost surely since we can write the columns of
$Z$ as $z_{i}=\Sigma_{i}^{1/2}w_{i}$ in distribution, where
$w_{i}\stackrel{{\scriptstyle\textnormal{ind.}}}{{\sim}}\mathsf{N}(0,I_{d})$.
Consequently,

$\operatorname*{\mathbf{P}}\ilimits@\Big{\\{}\left\|Z\right\|_{\rm
op}\geqslant\lambda_{\rm
max}\big{(}\sqrt{n}+\sqrt{d}+t\big{)}\Big{\\}}\leqslant\operatorname*{\mathbf{P}}\ilimits@\Big{\\{}\left\|W\right\|_{\rm
op}\leqslant\sqrt{n}+\sqrt{d}+t\Big{\\}}\leqslant e^{-t^{2}/2},$

by standard concentration bounds for Gaussian random matrices. ∎

We record the following corollary. To do so, we first consider the following
‘good’ event, when the operator norm of the noise matrix $G$ is small enough.

$\mathcal{E}\coloneq\Big{\\{}\left\|G\right\|_{\rm
op}\leqslant\sqrt{2\gamma}(\sqrt{n}+\sqrt{d})\Big{\\}}$ (6)

######  Corollary 2.

We have $\operatorname*{\mathbf{P}}\ilimits@(\mathcal{E})\geqslant
1-\mathrm{e}^{-0.08n}$.

###### Proof.

In Lemma 3, we may set $Z=G$, and take $\lambda_{\rm
max}\coloneq\sqrt{\gamma}$. The result follows by taking
$t=(\sqrt{2}-1)(\sqrt{d}+\sqrt{n})\geqslant\sqrt{0.4n}$. ∎

To finish the proof, we first introduce an alternative algorithm which is
simpler to analyze, which we refer to in the sequel as Algorithm (A.1).

  * Algorithm 3 [Clustering via low-rank approximation] _Algorithm for estimating clusters from samples._

input: sample matrix $X\in\mathbf{R}^{d\times n}$, number of clusters, $k$
---
|  1\. _Compute truncated singular value decomposition of sample matrix._
|
$X_{k}\coloneq\sum_{i=1}^{k}\widehat{\sigma}_{i}\widehat{u}_{i}\widehat{v}_{i}^{\mathsf{T}}\eqcolon\widehat{U}\widehat{\Sigma}\widehat{V}^{\mathsf{T}}\in\mathbf{R}^{d\times
n}$
|  2\. _Solve $k$-means on the columns of $Y$._ Compute a solution to the
optimization problem
|
$\mbox{minimize}\quad\frac{1}{2n}\sum_{i=1}^{n}\|X_{i}-\theta_{\kappa_{i}}\|_{2}^{2}$,
|  with variables $\theta_{1},\dots,\theta_{k}\in\mathbf{R}^{d}$ and
$\kappa_{1},\dots,\kappa_{n}\in[k]$.
output: estimates for labelling $\widehat{\mathsf{c}}_{\rm
appx}(\cdot)\colon[n]\to[k]$, where $\widehat{\mathsf{c}}_{\rm
appx}(i)\coloneq\kappa_{i}^{\star}$, for all $i$.

It turns out that there is a bijection $\phi\colon[k]\to[k]$ such that
$\widehat{\mathsf{c}}_{\rm appx}(i)=\phi(\widehat{\mathsf{c}}_{\rm spec}(i))$
for each $i\in[n]$ (for additional details, see Lemma 4.1 of [30]). Therefore,
by definition of $\ell(\cdot,\cdot)$, we have that

$\ell(\widehat{\mathsf{c}}_{\rm
spec},\mathsf{c}^{\star})=\ell(\widehat{\mathsf{c}}_{\rm
appx},\mathsf{c}^{\star}),$

and so it suffices to analyze the output of the algorithm that simply uses the
rank-$k$ approximant to $X$, algorithm A.1.

To do so, let us first introduce the following two sets of indices:

$\mathcal{I}_{\rm
bad}\coloneq\Big{\\{}i\in[n]:\|\theta^{\star}_{\mathsf{c}^{\star}(i)}-\widehat{\theta}_{i}\|_{2}\geqslant\Delta/2\Big{\\}}\qquad\mbox{and}\qquad\mathcal{I}_{\rm
good}\coloneq[n]\setminus\mathcal{I}_{\rm bad}.$

Above, $\widehat{\theta}_{i}$ is the optimal solution $\theta_{i}^{\star}$ in
step 2 of Algorithm A.1. We will also use the notation
$\widehat{\Theta}\in\mathbf{R}^{d\times n}$ to denote the matrix which has
columns $\widehat{\theta}_{i}$. Note that $\widehat{\Theta}$ by construction
is the projection of the rank-$k$ approximant $X_{k}$ onto the set of matrices
with at most $k$ distinct columns, in the Frobenius norm.

We then have the following proposition.

######  Proposition 1.

The cardinality of the set of bad indices enjoys the following upper bound:

$|\mathcal{I}_{\rm bad}|\leqslant\frac{128k\left\|G\right\|_{\rm
op}^{2}}{\Delta^{2}},$

almost surely.

###### Proof.

Note that

$\left\|\widehat{\Theta}-\Theta^{\star}\right\|_{\rm
F}^{2}=\sum_{i=1}^{n}\|\theta^{\star}_{\mathsf{c}^{\star}(i)}-\widehat{\theta}_{i}\|_{2}^{2}\geqslant|\mathcal{I}_{\rm
bad}|\frac{\Delta^{2}}{4}.$

Therefore, we see that

$|\mathcal{I}_{\rm
bad}|\leqslant\frac{4}{\Delta^{2}}\left\|\widehat{\Theta}-\Theta^{\star}\right\|_{\rm
F}^{2}$ (7)

Note that since $\widehat{\Theta},\Theta^{\star}$ are both (at most) rank-$k$
matrices we evidently have

$\left\|\widehat{\Theta}-\Theta^{\star}\right\|_{\rm
F}\leqslant\sqrt{2k}\left\|\widehat{\Theta}-\Theta^{\star}\right\|_{\rm
op}\leqslant\sqrt{8k}\left\|\Theta^{\star}-X_{k}\right\|_{\rm op},$

where in the last inequality we use the optimality of $\widehat{\Theta}$ along
with the triangle inequality. Adding and subtracting $X$, and then relying on
Eckart-Young, the display above implies

$\left\|\widehat{\Theta}-\Theta^{\star}\right\|_{\rm
F}\leqslant\sqrt{32k}\left\|X-\Theta^{\star}\right\|_{\rm
op}\leqslant\left\|G\right\|_{\rm op}\sqrt{32k}.$

Applying the display above in the cardinality bound (7) gives the result. ∎

Additionally, we define the following sets

$I_{j}\coloneq\Big{\\{}i\in[n]:\mathsf{c}^{\star}(i)=j,~{}\|\widehat{\theta}_{i}-\theta^{\star}_{\mathsf{c}^{\star}(i)}\|_{2}<\Delta/2\Big{\\}},\qquad
j=1,\ldots,k.$

Clearly, we have that $\\{I_{j}\\}$ form a disjoint partition of
$\mathcal{I}_{\rm good}$. However, as the following result shows, these sets
are also pairwise disjoint and nonempty under the estimated clustering,
$\widehat{\mathsf{c}}_{\rm appx}(\cdot)$.

######  Lemma 4.

The sets $\widehat{\mathsf{c}}_{\rm appx}(I_{j})$ are disjoint for each
$j\in[k]$.

###### Proof.

By contradiction, suppose that $\widehat{\mathsf{c}}_{\rm
appx}(I_{j_{1}})\cap\widehat{\mathsf{c}}_{\rm appx}(I_{j_{2}})\neq\emptyset$,
for some $j_{1}\neq j_{2}\in[k]$. This happens if and only if there exists
$i_{1},i_{2}\in\mathcal{I}_{\rm good}$ and $j_{0}\in[k]$ for which:

$\mathsf{c}^{\star}(i_{1})=j_{1},\quad\mathsf{c}^{\star}(i_{2})=j_{2},\quad\mbox{and}\quad
j_{0}=\widehat{\mathsf{c}}_{\rm appx}(i_{1})=\widehat{\mathsf{c}}_{\rm
appx}(i_{2}).$ (8)

However, then we note that

$\displaystyle\|\theta^{\star}_{j_{1}}-\theta^{\star}_{j_{2}}\|_{2}$
$\displaystyle\leqslant\|\theta^{\star}_{\mathsf{c}^{\star}(i_{1})}-\widehat{\theta}_{\widehat{\mathsf{c}}_{\rm
appx}(i_{1})}\|_{2}+\|\widehat{\theta}_{\widehat{\mathsf{c}}_{\rm
appx}(i_{1})}-\widehat{\theta}_{\widehat{\mathsf{c}}_{\rm
appx}(i_{2})}\|_{2}+\|\widehat{\theta}_{\widehat{\mathsf{c}}_{\rm
appx}(i_{2})}-\theta^{\star}_{\mathsf{c}^{\star}(i_{2})}\|_{2}$ (by triangle
inequality, conditions (8))
$\displaystyle=\|\theta^{\star}_{\mathsf{c}^{\star}(i_{1})}-\widehat{\theta}_{j_{0}}\|_{2}+\|\widehat{\theta}_{j_{0}}-\theta^{\star}_{\mathsf{c}^{\star}(i_{2})}\|_{2}$
(since $j_{0}=\widehat{\mathsf{c}}_{\rm appx}(i_{1})=\widehat{\mathsf{c}}_{\rm
appx}(i_{2})$) $\displaystyle<\Delta,$ (since $i_{1},i_{2}\in\mathcal{I}_{\rm
good}$)

which contradicts the definition of $\Delta$. ∎

##  Appendix B Other deferred proofs

###### Proof of Lemma 1.

Following the normal equations for the scalar autoregression problem (2),

$\widehat{\theta}^{\rm
AR}_{i}=\theta_{i}+\widehat{\Sigma}_{i}^{-1}\bigg{\\{}\frac{1}{T}\sum_{t=1}\mathsf{T}(\varepsilon_{i}^{(t)}+\Delta_{i}^{(t)})x_{i}^{(t)}\bigg{\\}}.$
(9)

Above we defined
$\Delta_{i}^{(t)}\coloneq\sum_{j:c_{j}=c_{i}}\gamma_{ij}^{\mathsf{T}}x_{j}^{(t)}$.
Note that

$\Delta_{i}^{(t)}\stackrel{{\scriptstyle\rm
d}}{{=}}\mathsf{N}\Big{(}0,(\tau\rho^{(t)}_{i})^{2}\Big{)},\qquad\mbox{where}\quad\rho_{i}^{(t)}\coloneq\bigg{(}\sum_{j:c_{j}=c_{i}}\|x_{j}^{(t)}\|_{2}^{2}\bigg{)}^{\\!\tfrac{1}{2}}.$

Therefore in view of (9), we see that $\widehat{\theta}^{\rm AR}_{i}$ is
indeed a sum of Gaussians with mean $\overline{\theta}_{c_{i}}$ and covariance

$\Lambda_{i}=\nu^{2}I_{d}+\frac{\sigma^{2}}{T}\widehat{\Sigma}_{i}^{-1}+\frac{\sigma^{2}}{T}\widehat{\Sigma}_{i}^{-2}\widehat{\Sigma}_{i}(\rho_{i}),$

which proves the result. ∎

###### Proof of Lemma 2.

To lighten notation, we define

$Z_{t}\coloneq I_{m}\otimes(x^{(t)})^{\mathsf{T}},\quad t\in[T].$

Note that $Z_{t}\in\mathbf{R}^{m\times m^{2}d}$. It follows from the
optimality condition for the VAR problem that we have the following
decomposition:

$\Gamma^{\rm
VAR}=\bigg{(}\frac{1}{T}\sum_{t=1}^{T}Z_{t}^{\mathsf{T}}Z_{t}\bigg{)}^{-1}\frac{1}{T}\sum_{t=1}^{T}Z_{t}^{\mathsf{T}}\varepsilon^{(t)}=\Gamma^{\star}+\frac{1}{T}\sum_{t=1}^{T}Z_{t}^{\mathsf{T}}\varepsilon^{(t)},$

where above we use the isotropy assumption. The decomposition above then
immediately implies that

$\Gamma^{\rm
VAR}\sim\mathsf{N}(\Gamma^{\star},\frac{\sigma^{2}}{T}I_{m^{2}d}).$

Note that $\tfrac{T}{\sigma^{2}}\|(\Gamma^{\rm VAR}-\Gamma^{\star})\|_{2}^{2}$
is equal in distribution to a $\chi^{2}$ random variable with $m^{2}d$-degrees
of freedom. The result now follows by using standard concentration bounds for
$\chi^{2}$-random variables. ∎

##  Appendix C Implementation Details

Simulation Study. In the simulation experiment in Figure 2, we generate a
synthetic dataset, where there are $n$ time series grouped into $k$ clusters
of equal size. In cluster $i$, the local AR parameters for every time series
is centered at $\theta_{c(i)}=z_{i}/C_{i}$ where $z_{i}$ is generated from
$\mathsf{N}(0,I_{d})$ distribution and $C_{i}$ is a constant. The noise around
$\theta_{c(i)}$ is set to have a standard deviation of 1e-2. The VAR
parameters for time series in a cluster is also generated from the standard
normal distribution. In other words $\gamma_{i,j}=z_{ij}/C_{ij}$ where
$z_{ij}\sim\mathsf{N}(0,I_{d})$ if $c(i)=c(j)$. The constants $C_{i}$’s and
$C_{ij}$’s are chosen such that the $\ell_{p}$ norm of all the lag
coefficients of any particular time series is one for $p=2.5$. The choice of
the normalization is made to promote stationarity. We also restrict the values
of the time series between $[-1,1]$. Thus the time-series pregression is
governed by the equation,

$\displaystyle
y_{i}^{(t)}=\min\left(1,\max\left(\theta_{i}^{\mathsf{T}}y^{(t-d:t-1)}_{i}+\sum_{j:c(i)=c(j)}\gamma_{ij}^{\mathsf{T}}y_{j}^{(t-d:t-1)}+\varepsilon_{i}^{(t)},-1\right)\right).$
(10)

Learning Algorithm. The scalar AR algorithm is implemented using one LinearSVR
[18] instance per time series. Similarly the VAR models are also implemented
using one LinearSVR instance per cluster. We tune the parameters of the AR
model on a randomly chosen time series from the whole datasets and then held
fixed for all other time series. The hyper-parameters are tuned via grid-
search among the grid {"max_iter": [1000], "C": np.logspace(-4, 0, 10),
"epsilon": [0, 1e-1, 1e-2, 1e-3]}. Similarly the VAR hyper-parameters are
tuned on a randomly chosen cluster among the same grid parameters as above,
and then held fixed.

Lag Indices. These are very important parameters for both the AR and VAR
model. For hourly datasets like Electricity and Traffic, we use the lag
indices np.concatenate([np.arange(1, 25), np.arange(7*24, 8*24),
np.arange(14*24, 15*24)]) i.e the last 24 hours, 24 hours on the same day last
week and 24 hrs on the same day two weeks back. For the datasets in 5min
intervals, like ECG and PEMS07 we set the lag indices as
np.concatenate([np.arange(1, 12), np.arange(12, 15), np.arange(24, 27)]) i.e
the last hour, the corresponding 15 min in the last hour and the corresponding
15 minutes two hours back.

Clustering. We cluster on the learnt AR parameters normalized to $\ell_{2}$
norm one. We fit a KD-Tree [7] in the AR paremeter set and form a $K$-nn graph
using 11 neighbors for each point. Then we use METIS graph clustering on this
graph. In general we set the number of clusters to be $k=\left\lfloor
n/10\right\rfloor$ for a dataset with $n$ time series. We tuned the $k$ among
$\\{\left\lfloor n/5\right\rfloor,\left\lfloor n/10\right\rfloor,\left\lfloor
n/20\right\rfloor\\}$, and $k=\left\lfloor n/10\right\rfloor$ was close to
optimal among the three on all 4 datasets.

Baselines. Now we can present details about the baselines. For the first two
experiments i.e the ones on ECG and PEMS07 all the baselines are in a setting
identical to the ones used in [9]. We merely recreate the values from the
table in [9]. It is a fair comparison as our algorithms are also evaluated on
the exact same task. Similary for the experiments onTraffic and Electricity
dataset the baseline numbers are taken from Table 2 in [36]. We make these
choices for the ease of reproducibility, moving forward.

Dataset | Num. Series ($n$) | Num. Training Points ($t$) | Forecast Horizon ($\tau$) | Rolling Windows ($n_{w}$) | Granularity
---|---|---|---|---|---
Electricity | 370 | 25,968 | 24 | 7 | Hourly
Traffic | 963 | 10,392 | 24 | 7 | Hourly
PEMS07 | 208 | 11,232 | 3 | 480 | 5min
ECG | 140 | 4,999 | 3 | 167 | 5min

Table 5: We present the salient features of the datasets. We use the last
$n_{w}\times\tau$ time-points to construct the test sets. Then, we perform
rolling validation with a forecast horizon of $\tau$ and number of rolling
windows as $n_{w}$.

##  Appendix D Metrics

The following well-known loss metrics are used in this paper. Here,
$\mathbf{Y}\in\mathbf{R}^{n^{\prime}\times t^{\prime}}$ represents the actual
values while $\widehat{\mathbf{Y}}\in\mathbf{R}^{n^{\prime}\times t^{\prime}}$
are the corresponding predictions.

(i) WAPE:  Weighted Absolute Percent Error is defined as follows,

$\displaystyle\mathcal{L}_{w}(\widehat{\mathbf{Y}},\mathbf{Y})=\frac{\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}|Y_{ij}-\widehat{Y}_{ij}|}{\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}|Y_{ij}|}.$
(11)

(ii) MAPE:  Mean Absolute Percent Error is defined as follows,

$\displaystyle\mathcal{L}_{m}(\widehat{\mathbf{Y}},\mathbf{Y})=\frac{1}{Z_{0}}\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}\frac{|Y_{ij}-\widehat{Y}_{ij}|}{|Y_{ij}|}\mathbf{1}\\{|Y_{ij}|>0\\},$
(12)

where
$Z_{0}=\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}\mathbf{1}\\{|Y_{ij}|>0\\}$.

(iii) SMAPE:  Symmetric Mean Absolute Percent Error is defined as follows,

$\displaystyle\mathcal{L}_{s}(\widehat{\mathbf{Y}},\mathbf{Y})=\frac{1}{Z_{0}}\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}\frac{2|Y_{ij}-\widehat{Y}_{ij}|}{|Y_{ij}|+|\widehat{Y}_{ij}|}\mathbf{1}\\{|Y_{ij}|+|\widehat{Y}_{ij}|>0\\},$
(13)

where
$Z_{0}=\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}\mathbf{1}\\{Y_{ij}|+|\widehat{Y}_{ij}|>0\\}$.

(iv) MAE:  Mean Absolute Error is defined as follows,

$\displaystyle\mathcal{L}_{mae}(\widehat{\mathbf{Y}},\mathbf{Y})=\frac{1}{n^{\prime}t^{\prime}}\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}|Y_{ij}-\widehat{Y}_{ij}|$
(14)

(iv) RMSE:  Root Mean Squared Error is defined as follows,

$\displaystyle\mathcal{L}_{rmse}(\widehat{\mathbf{Y}},\mathbf{Y})=\sqrt{\frac{1}{n^{\prime}t^{\prime}}\sum_{i=1}^{n^{\prime}}\sum_{j=1}^{t^{\prime}}(Y_{ij}-\widehat{Y}_{ij})^{2}}.$
(15)

###  D.1 Linear Generative Model and Analogy to Mixed Linear Regression

For the sake of completeness, we define a linear generative model for our
time-series setting under which we expect our algorithm to recover the true
parameters. In the second half of this section, we will relate this model to
the simplified mixed linear regression problem, under which we formally prove
our guarantees.

We will consider an auto-regressive lag of $d$. That is, any data point at a
given time can depend only on quantities no more than $d$ time steps in the
past. We consider a dataset with $n$ time series and $T$ time points. The
following process selects the local autoregressive parameters for each time
series $i\in[n]$:

$\displaystyle\mathsf{c}^{\star}(i)$
$\displaystyle\sim\mathsf{Categorical}(\pi_{1},\dots,\pi_{k}),\quad\mbox{and}\quad$
(16a) $\displaystyle\theta_{i}\mid\mathsf{c}^{\star}(i)$
$\displaystyle\sim\mathsf{N}_{d}(\overline{\theta}_{\mathsf{c}^{\star}(i)},\nu^{2}I_{d}),$
(16b)

where $\mathsf{Categorical}$ and $\mathsf{N}_{d}$ denote the Categorical and
$d$-dimensional Gaussian distributions respectively.

This means that the time series $i$ is first randomly assigned into one of $k$
clusters $\mathsf{c}^{\star}(i)$. Based on the cluster, the local
autoregressive parameters of the time series are drawn from the
$\mathsf{c}^{\star}(i)$-th component of a Gaussian Mixture model in
$\mathbf{R}^{d}$. The evolution of the time series is then given as follows:

$x_{i}^{(s)}=\theta_{i}^{\mathsf{T}}x_{i}^{(s-d:s-1)}+\sum_{\begin{subarray}{c}j\neq
i\\\
\mathsf{c}^{\star}(i)=\mathsf{c}^{\star}(j)\end{subarray}}\gamma_{ij}^{\mathsf{T}}x_{j}^{(s-d:s-1)}+\varepsilon_{i}^{(s)},$
(17)

for $s\in\\{d+1,\cdots,T\\}$. Here, $\gamma_{ij}\in\mathbf{R}^{d}$ are
parameters that govern how time series $i$ depends on the past of time series
$j$, only if both of them belong to the same cluster. This essentially implies
that apart from the dependence on its own past, a time series can depend on
the past of all other time series in the same cluster. Moreover the clustering
is defined by closeness with respect to the local autoregressive parameters.
From Equation (17), we see that the values of a time series do not depend on
those from clusters that the corresponding time series does not belong to.
This is a reasonable assumption in real-world multivariate time-series
settings, where each time series would only depend on a few related time
series.

Relation to MLR: The mixed linear regression problem is an i.i.d analogue of
the linear generative model presented above. Each time-series can be thought
of as a scalar regression problem, where each data-point in the MLR problem
corresponds to the target being the future time-point and the covariates being
the past $d$ time-points, i.e

$\displaystyle(x_{i}^{(j)},y_{i}^{(j)})\sim(x^{(s_{j}-d:s_{j}-1)}_{i},x_{i}^{(s_{j})}).$
(18)

where $s_{j}$ corresponds to the blocks separated by the mixing time, defined
in Section 3. Time-series that come from the same cluster would have
$\theta_{i}$’s from the same component of the mixture in MLR model defined in
Section 4. The time-series in the same cluster will have further VAR
dependencies on each other as captured by $\gamma$ coefficients in the MLR
model.