# TSUBASA: Climate Network Construction
on Historical and Real-Time Data

Yunlong Xu  ,  Jinshu Liu  University of Rochester<EMAIL_ADDRESS><EMAIL_ADDRESS>and  Fatemeh Nargesian  University of Rochester
<EMAIL_ADDRESS>
(2022)

###### Abstract.

A climate network represents the global climate system by the interactions of
a set of anomaly time-series. Network science has been applied on climate data
to study the dynamics of a climate network. The core task and first step to
enable interactive network science on climate data is the efficient
construction and update of a climate network on user-defined time-windows. We
present TSUBASA, an algorithm for the efficient construction of climate
networks based on the exact calculation of Pearson’s correlation of large
time-series. By pre-computing simple and low-overhead statistics, TSUBASA can
efficiently compute the exact pairwise correlation of time-series on arbitrary
time windows at query time. For real-time data, TSUBASA proposes a fast and
incremental way of updating a network at interactive speed. Our experiments
show that TSUBASA is faster than approximate solutions at least one order of
magnitude for both historical and real-time data and outperforms a baseline
for time-series correlation calculation up to two orders of magnitude.

time-series, climate data, climate networks, correlation matrix

††copyright: acmcopyright††journalyear: 2022††copyright:
acmcopyright††conference: Proceedings of the 2022 International Conference on
Management of Data; June 12–17, 2022; Philadelphia, PA, USA††booktitle:
Proceedings of the 2022 International Conference on Management of Data (SIGMOD
’22), June 12–17, 2022, Philadelphia, PA, USA††price: 15.00††doi:
10.1145/3514221.3526177††isbn: 978-1-4503-9249-5/22/06

##  1\. Introduction

To identify and analyze patterns in global climate, scientists and climate
risk analysts model climate data as complex networks – networks with non-
trivial topological properties (Abe and Suzuki, 2004; Kim et al., 2019;
Gozolchiani et al., 2008). The climate network architecture represents the
global climate system by a set of anomaly time-series (departure from the
usual behavior) of gridded climate data and their interactions (Tsonis and
Roebber, 2004). A climate data set includes remote and in-situ sensor
measurements (e.g. sea surface temperature and sea level pressure) covering a
grid (e.g. with a resolution of $2.5^{\circ}\times 2.5^{\circ}$). Nodes in a
climate network are geographical locations, characterized by time-series and
edges represent information flow between nodes. The edge weights indicate a
degree of correlation between the behaviors of time-series (e.g. Pearson’s
correlation). Note the geographical locality of nodes does not directly imply
the topology of a network.

Several studies have applied network science on climate data assuming dynamic
networks that are changing with real-time data (Berezin et al., 2012). Climate
networks have been shown to be powerful tools for gaining insights on
earthquakes (Abe and Suzuki, 2004), rainfalls (Kim et al., 2019), and global
climate events such as El Niño (Gozolchiani et al., 2008).

The common way for network dynamics analysis is to construct networks for each
hypothesized time-window and analyze them separately (Faghmous and Kumar,
2014). Figure 1 shows the steps of constructing a climate network. Given a
query window provided by a user, a correlation matrix is constructed by
computing the pairwise correlation of all time-series on the query window.
Pearson’s correlation is one of the most dominant measures for studying the
pairwise climatical correlation (Donges et al., 2009). The correlation matrix
enables visualization (Nocke et al., 2015), network dynamics analysis (Berezin
et al., 2012), as well as tasks such as community detection (Tantet and
Dijkstra, 2014). To analyze the topology of the network, a user-provided
correlation threshold can be applied on the matrix to find the significant
edges between nodes and obtain a boolean network matrix. From the mathematical
perspective, the analytical computation of the evolution of a complex system
(or even not so complex such as Ordinary Differential Equation systems),
depends on the robustness and correctness of the initial weights in the
complex network (Gozolchiani et al., 2008).

The core task in network construction is the problem of large-scale all-pair
time-series correlation calculation. The key challenges of interactive network
analysis include: 1) exact calculation of the complete correlation matrix, 2)
correlation calculation on time-windows of arbitrary size, and 3) efficiency
of network construction and update for historic and real-time data to achieve
interactivity.

Figure 1. Climate Network Construction.

The line of data management research that computes networks on time-series
(for example, for stock market data or climate data) apply pruning techniques
on the approximation of correlation measures. In particular, StatStream (Zhu
and Shasha, 2002) and MASS (Mueen et al., 2010) reduce the correlation of
time-series to the distance of their Discrete Fourier Transform (DFT)
coefficients and propose grid-based indexing (Zhu and Shasha, 2002) and
I/O-aware techniques (Mueen et al., 2010) for performing threshold-based
correlated time-series search. The accuracy of the network can be increased by
considering a very large number of DFT coefficients that are expensive to
compute. a query window size is restricted to be an integral multiple of the
basic window size, which limits the usability of these algorithms. In this
paper, we present TSUBASA a framework for efficient construction and update of
the correlation matrix for arbitrary query windows on historical and real-time
data. The differences of TSUBASA and existing work are three-fold. First,
unlike these algorithms that only identify edges with a correlation higher
than a threshold, TSUBASA computes and updates the complete correlation
matrix. This enables network dynamics analysis as well as choosing arbitrary
thresholds at query time. Second, TSUBASA computes the exact correlation
between time-series. Finally, unlike the existing work, TSUBASA allows
arbitrary query window sizes. In this paper, we make the following
contributions.

  * $\bullet$

We present the mathematical tools for the exact calculation of pairwise
Pearson’s correlation of time-series using the basic window model.

  * $\bullet$

TSUBASA relies on simple statistics of basic windows pre-computed by doing one
pass over the whole data. This provides a flexible and highly responsive
correlation calculation mechanism. Users can obtain a correlation matrix given
any query window without computing the correlation statistics repeatedly.

  * $\bullet$

We propose an incremental solution for real-time update of the correlation
matrix and climate network. Relying on the easy-to-compute statistics of basic
windows means the correlations can be updated quickly for frequently-updated
time-series.

  * $\bullet$

We enable queries with arbitrary time-window size and start and end point on
both historic and real-time by relaxing the restriction of the existing basic
window model on a query window size being an integral multiple of basic window
size.

Figure 2. Query Windows for Historical and Real-time Data.

##  2\. Preliminaries

We start by introducing the climate network construction problem, then, give
an overview of existing approximate techniques for calculating time-series
correlation.

###  2.1. Climate Network

We are given a collection $\mathcal{L}=\\{x^{1},\ldots,x^{n}\\}$ of geo-
labeled time-series, where $x^{i}$ denotes the time-stamped values of a
climatic variable collected at location $i$. A time-series $x^{i}$ is defined
as $[\textbf{x}^{i}_{1},\ldots,\textbf{x}^{i}_{m}]$, where
$\textbf{x}^{i}_{j}$ is the observed value at time $j$. We assume all time-
series in $\mathcal{L}$ are synchronized, i.e. each time-series has a value
available at every periodic time interval, namely time resolution.
Particularly, if the time resolution is $\gamma$ and the current timestamp is
$j$, every $x^{i}$ is $\mathcal{L}$ will have a value observed at time
$j+\gamma$. If an $x^{i}$ has missing value at $j$, a value is interpolated or
if multiple values appear between $j$ and $j+\gamma$, an aggregate value is
assigned. Table 1 shows a list of notations used throughout this paper.

Now, we turn our attention to user’s query. At query time, a user defines a
query time-window $w=(e,l)$ and a correlation threshold $\theta$. The query
window is defined with an end timestamp $e$ and a length $l$ that indicates a
sub-sequence of size $l$ in a time-series with a start timestamp $e-l+1$ and
an end timestamp $e$. For real-time data, the end timestamp can be the last
observed time, i.e. a user query on real-time data is $w=(``now",l)$, that
means the network is constructed on the last $l$ observed data points. We
consider the data points within $w$ for each time-series
$x=[\textbf{x}_{1},\ldots,\textbf{x}_{k}]$. For example,
$[\textbf{x}_{k-m+1},\ldots,\textbf{x}_{m}]$ is the sequence we consider for
$x$ on the query window $w=(k,m)$. When clear from the context, we call the
sequence of a time-series $x$, for a given query window simply query window or
time-series $x$. The climate network of $\mathcal{L}$ for a given query time-
window $w$ is a graph $\mathcal{N}=(G,V)$, where a node in $G$ corresponds to
a location $i$ and is represented by time-series $x^{i}$. An edge in $V$
between nodes $i$ and $l$ indicates that the correlation between time-series
$x^{i}$ and $x^{l}$ is above a user-defined threshold $\theta$. In this paper,
we focus on the most commonly used correlation measure i.e. Pearson’s
correlation coefficients. Given query windows
$x=[\textbf{x}_{1},\ldots,\textbf{x}_{m}]$ and
$y=[\textbf{y}_{1},\ldots,\textbf{y}_{m}]$, with means $\bar{x}$ and
$\bar{y}$, the Pearson’s correlation of $x$ and $y$ is calculated as follows.

(1) ${\displaystyle
Corr(x,y)={\frac{\sum_{i=1}^{m}(\textbf{x}_{i}-{\bar{x}})(\textbf{y}_{i}-{\bar{y}})}{{\sqrt{\sum_{i=1}^{m}(\textbf{x}_{i}-{\bar{x}})^{2}}}{\sqrt{\sum_{i=1}^{m}(\textbf{y}_{i}-{\bar{y}})^{2}}}}}}$

Constructing a network for a query window (e.g. the first six months of 2021)
requires computing the correlation for all pairs of time-series and pruning
links using threshold $\theta$. Existing techniques for fast all-pair
correlation calculation on large time series approximate pairwise correlation
by relying on the Fourier transform (Mueen et al., 2010; Zhu and Shasha, 2002;
Cole et al., 2005). In this paper, we present an efficient way of calculating
the exact correlation of time-series and constructing a network for historical
data and updating the network for real-time data. The existing work divide
time-series into cooperative and uncooperative to perform correlation
approximation. Although our core contribution is the exact calculation and the
update of correlations, we also present an approximate way of calculating
correlation for generic time-series (§ 3.2). We start by giving an overview of
correlation approximation.

###  2.2. Correlation Approximation Solutions

Computing the correlation of large-scale time-series is pervasively done using
the notion of basic windows (Zhu and Shasha, 2002; Mueen et al., 2010). Time-
series are processed in batches of size $B$, i.e. the stream
$[{\textbf{x}}_{1},\ldots,{\textbf{x}}_{n}]$ is equally divided into $n/B$
basic windows, where the $j$-th basic window contains data
$[{\textbf{x}}_{(j-1)*B},\ldots,{\textbf{x}}_{j*B}]$. Similarly, a query
window is a sequence of basic windows. A query window is assumed to be
divisible by the size of basic window. Figure 2 presents a visualization of
query and basic windows. The existing techniques approximate correlation using
the Discrete Fourier Transform (DFT) of basic windows. The DFT of a time-
series $x=[\textbf{x}_{1},\ldots,\textbf{x}_{k}]$ is a sequence
$X=[\textbf{X}_{1},\ldots,\textbf{X}_{k}]$ of complex numbers:

(2)
$\textbf{X}_{f}=\frac{1}{\sqrt{k}}\sum_{i=1}^{k}\textbf{x}_{i}e^{\frac{-j2\pi
fi}{k}},f=1,\ldots,k$

where $j=\sqrt{-1}$. Computing DFT coefficients has time complexity $O(n^{2})$
in the size of basic window. For normalized time series, DFT preserves the
Euclidean distance between two sequences, that is,
$Dist(\hat{x},\hat{y})=Dist(\hat{X},\hat{Y})$. The approximation techniques
consider the first few DFT coefficients to capture the shape and properties of
time-series. It has been shown that the correlation of two time-series can be
reduced to the Euclidean distance of the the DFT coefficients of their
normalized time-series (Rafiei, 1999; Zhu and Shasha, 2002). The normalization
of a basic window $x_{i}=[\textbf{x}_{1},\ldots,\textbf{x}_{B}]$ is
$\hat{x_{i}}=[\frac{\textbf{x}_{1}-\overline{x_{i}}}{\sigma_{i}},\ldots,\frac{\textbf{x}_{B}-\overline{x_{i}}}{\sigma_{i}}]$,
where $\overline{x}_{i}$ and $\sigma_{i}$ are the mean and standard deviation
of $x_{i}$. The correlation of two time-series can be obtained from the
Euclidean distance $d(.,.)$ of their normalized series.

(3) $c_{i}=1-\frac{1}{2}d^{2}(\hat{x_{i}},\hat{y_{i}})$

For more concise notation, we denote $d_{i}$ to be
$d(\hat{x_{i}},\hat{y_{i}})$. Suppose $\hat{X}_{i}$ and $\hat{Y}_{i}$ are the
DFT of normalized basic windows $\hat{x}_{i}$ and $\hat{y}_{i}$, and
$Dist_{n}(\hat{X}_{i},\hat{Y}_{i})$ is the Euclidean distance of the first $n$
DFT coefficients in $\hat{X}_{i}$ and $\hat{Y}_{i}$. Recall DFT preserves the
distance between coefficients and the original time-series. Therefore,
$d_{i}\simeq Dist_{n}(\hat{X}_{i},\hat{Y}_{i})$. The more coefficients are
used (the higher $n$), the more accurate the distance and correlation becomes.
So far, we have a way of computing the distance of basic windows. To compute
the distance of query windows, $Dist_{n}(x,y)$, the existing techniques assume
that the form and properties of time-series do not drastically change over a
query window, i.e. basic windows have similar statistics (mean and standard
deviation) to the query window (Zhu and Shasha, 2002; Mueen et al., 2010).
When the statistics do not change, $Dist_{n}(x,y)$ is the average of the
$d_{i}$ on all basic windows of $x$ and $y$. In § 3.2, we relax this
assumption and consider time series that change in form and properties over a
query window, i.e. the statistics of basic windows are not necessarily similar
to each other and the query window. Now, we apply Equation 3 on time windows,
to get $Corr(x,y)\simeq 1-\frac{1}{2}Dist_{n}(x,y)$. Again the higher $n$ we
use, the better approximation of correlation we obtain.

Now, we describe how $Dist_{n}(\hat{X},\hat{Y})$ is used to decide whether
$Corr(x,y)\geq\theta$. Zhu and Shasha show the following relationship between
correlation and the distance based on $n$ DFT coefficients of $\hat{X}$ and
$\hat{Y}$.

(4) $Corr(x,y)\geq 1-\epsilon^{2}\Rightarrow
Dist_{n}(\hat{X},\hat{Y})\leq\epsilon$

When using approximate techniques for climate network construction, to get the
pairs of time-series with $Corr(x,y)\geq\theta$, we can compute
$\epsilon=\sqrt{1-\theta}$. This allows us to prune pairs with condition
$Dist_{n}(\hat{X},\hat{Y})\leq\sqrt{1-\theta}$.

Using Equation 4, we get a superset of highly correlated time-series with no
false negatives. As we show in Figure 5(a), the false positives incur spurious
edges in the network and result in an inaccurate network. These false positive
edges can only be filtered at the cost of exact correlation calculation from
the raw data. To avoid false positives, TSUBASA calculates exact correlations
of time-series, even faster than approximation.

Applying Equation 4 requires normalizing time-series and calculating DFT
coefficients. When using DFT-based approximation, the accuracy of network
increases as more DFT coefficients are considered. Indeed, approximate
techniques consider very few coefficients (two in the case of StatStream (Zhu
and Shasha, 2002) for any basic window size). However, when dealing with
climate data sets, which are uncooperative time-series, the majority of
coefficients are needed to get near accurate results (Figure 5(a)). Statstream
proposes random projection for uncooperative time-series that similar to DFT
coefficient calculation approximates correlation and has high overhead. To
overcome modeling uncooperative time-series, Qiu et al. use Fourier transform
and neural network to embed time series into a low-dimensional Euclidean space
(Qiu et al., 2018). The search is done using a nearest neighbor search index
in the embedding space. In this paper, we choose a different approach and
compute exact correlation and network. Our solution, TSUBASA, relies on simple
statistics of basic windows, and as we show empirically is much faster than
approximation techniques. In particular, TSUBASA can update a network
constructed on real-time data more efficiently than approximation techniques.
This enables interactive network analysis on accurate networks for historical
and real-time climate data.

##  3\. Network Construction

Table 1. Table of Notations Symbol | Description
---|---
$x$ | a query window $[\textbf{x}_{1},\ldots,\textbf{x}_{m}]$ of stream x
$\textbf{x}_{i}$ | data value at time $i$ in a query window $x$
$\overline{x}$ | mean of $x$
$x_{j}$ | the $j$-th basic window of $x$
$\overline{x_{j}}$ | mean of the $j$-th basic window of $x$
$\overline{x_{i:j}}$ | mean of basic windows $x_{i},\ldots,x_{j}$ of $x$
$n_{s}$ | number of basic windows in a query window
$B$ | number of data points in a basic window
$\sigma_{xj}$ | standard deviation of the $j$-th basic window of $x$
$c_{j}$ | correlation of $x$ and $y$ on the $j$-th basic window

Before a deep dive into exact correlation calculation, we present a high-level
overview of TSUBASA’s end-to-end framework. Figure 3 illustrates the
components of TSUBASA for constructing and updating a climate network on
historical and real-time data. The data storage contains a collection of
frequently updated time-series accessible through locations. During the pre-
processing, every time-series is divided into basic windows. We sketch basic
windows of time-series, in one pass, and store the collected statistics. This
can be done at data ingestion time. At query time, the statistics of the basic
windows corresponding to a given query window of all time-series are retrieved
and all-pair correlation is calculated without the need to access the raw
data. For real-time data, the system constructs the initial network and
ingests the real-time raw data in chunks of size $B$. The sketching of the
newly ingested basic window is done on the fly and the correlations of time-
series are updated incrementally without computing the correlation from
scratch. In the following sections, we describe the details of data sketching
as well as the mathematical model for the exact and approximate calculation of
correlation.

Figure 3. Architecture of TSUBASA. update figure with parallel info.

###  3.1. Exact Pairwise Correlation

####  3.1.1. Historical Data

Algorithm 1 Preprocessing

1:streams $\mathcal{L}=\\{\textbf{x}^{1},\ldots,\textbf{x}^{n}\\}$; basic
window size $B$

2:statistics $S$

3:$n_{s}\leftarrow{\bf Len}(\textbf{x}^{1})/B$

4:$S\leftarrow\\{\\}$

5:for $x,y\in\mathcal{L}$ do

6: $x\leftarrow{\bf BasicWin}(\textbf{x},B)$; $y\leftarrow{\bf
BasicWin}(\textbf{y},B)$

7: for $j\in[1..n_{s}]$ do

8: $S_{x_{j}}\leftarrow{\bf Stats}(x_{j})$; $S_{y_{j}}\leftarrow{\bf
Stats}(y_{j})$

9: $c_{j}\leftarrow{\bf Corr}(x_{j},y_{j})$

10: $\hat{x}_{j}\leftarrow{\bf Normalize}(x_{j},S_{x_{j}})$;
$\hat{y}_{j}\leftarrow{\bf Normalize}(y_{j},S_{y_{j}})$

11: $\hat{X}_{j}\leftarrow{\bf DFT}(\hat{x}_{j})$; $\hat{Y}_{j}\leftarrow{\bf
DFT}(\hat{y}_{j})$

12: $d_{j}\leftarrow{\bf Dist_{n}}(\hat{X}_{j},\hat{Y}_{j})$

13: // 8-10 are performed for approximation method

14: $S\leftarrow$WriteStats($S_{x_{j}},S_{y_{j}},c_{j},d_{j}$)

15:return $S$

Algorithm 2 Network-Construct-Histo

1:streams $\mathcal{L}=\\{\textbf{x}^{1},\ldots,\textbf{x}^{n}\\}$; statistics
$S$; query $w$; basic window size $B$; threshold $\theta$

2:graph $(G,V)$

3:$G\leftarrow\\{1,\ldots,n\\}$; $V\leftarrow\\{\\}$

4:$b\leftarrow{\bf GetBasicWins}(w)$ // basic window ids in $w$

5:for $\textbf{x}\in\mathcal{L}$ and $\textbf{y}\in\mathcal{L}$ do

6: $S_{x}\leftarrow{\bf ReadStats}(S,b,x)$; $S_{y}\leftarrow{\bf
ReadStats}(S,b,y)$

7: $c\leftarrow{\bf Corr}(S_{x},S_{y})$ // use Lemma 1

8: if $c>\theta$ then

9: $V.{\bf Add}(x,y,c)$

10:return $(G,V)$

Our solution uses the basic window model to calculate the exact correlation of
times-series. Subdividing a series into basic windows allows us to process
data in smaller batches. Existing works for approximate correlation
calculation assume equal length across all basic windows (Zhu and Shasha,
2002; Mueen et al., 2010). Assuming an equal basic window length poses a
limitation on the length of the query window, that is the length of the query
window could only be an integral multiple of the length of the basic window.
Formally, we have $\mid x\mid=B\cdot n_{s}$, where $\mid x\mid$ is the length
of the basic window. TSUBASA relaxes this assumption and provides a way of
calculating the exact pairwise correlation of time-series for arbitrary query
window lengths.

###### Lemma 0.

Given query windows $x=[{\textbf{x}}_{1},\ldots,{\textbf{x}}_{m}]$ and
$y=[{\textbf{y}}_{1},\ldots,{\textbf{y}}_{m}]$ and the sizes of basic windows:
$\textbf{B}=[B_{1},B_{2},\ldots,B_{m}]$, where $B_{i}$ is the size of the
$i$th basic window size. The exact Pearson’s correlation of $x$ and $y$ is:

$\displaystyle Corr$
$\displaystyle(x,y)=\frac{\sum_{j=1}^{n_{s}}B_{j}(\sigma_{x_{j}}\sigma_{y_{j}}c_{j}+\delta_{x_{j}}\delta_{y_{j}})}{\sqrt{\sum_{i=1}^{n_{s}}B_{i}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}\sqrt{\sum_{i=1}^{n_{s}}B_{i}({\sigma_{y_{i}}}^{2}+{\delta_{y_{i}}^{2}})}}$
$\displaystyle\delta_{x_{i}}=\overline{x_{i}}-\frac{\sum_{k=1}^{n_{s}}\overline{x_{k}}}{n_{s}},~{}~{}\delta_{y_{i}}=\overline{y_{i}}-\frac{\sum_{k=1}^{n_{s}}\overline{y_{k}}}{n_{s}}$

where, $\sigma_{x_{i}}$ ($\sigma_{y_{i}}$) is the standard deviation of basic
window of $x_{i}$ ($y_{i}$), $c_{i}$ is the correlation of basic windows
$x_{i}$ and $y_{i}$, $\overline{x_{i}}$ ($\overline{y_{i}}$) is the mean of
basic window $x_{i}$ ($y_{i}$).

###### Proof.

This lemma has been provided as a possible general extension provided by
Dunlap (Dunlap, 1937), without proof. We provide a proof here. Let
$\Omega_{j}$ be the size of the tail of a time-series with $B_{1}$ to $B_{j}$
basic windows of arbitrary size.

$\displaystyle\Omega_{j}=\sum_{k=1}^{j}B_{k};\Omega_{0}=0$ $\displaystyle
Corr(x,y)=\frac{1}{T}\sum_{j=1}^{n_{s}}\sum_{i=\Omega_{j-1}+1}^{\Omega_{j}}(\frac{\textbf{x}_{i}-\overline{x}}{\sigma_{x}}).(\frac{\textbf{y}_{i}-\overline{y}}{\sigma_{y}})$
$\displaystyle=\frac{1}{T}\sum_{j=1}^{n_{s}}\sum_{i=\Omega_{j-1}+1}^{\Omega_{j}}\frac{\sigma_{x_{j}}\textbf{x}^{j,i}+\delta_{x_{j}}}{\sigma_{x}}.\frac{\sigma_{y_{j}}\textbf{y}^{j,i}+\delta_{y_{j}}}{\sigma_{y}}$
$\displaystyle=\frac{1}{T}\frac{\sum_{j=1}^{n_{s}}B_{j}(\sigma_{x_{j}}\sigma_{y_{j}}c_{j}+\delta_{x_{j}}\delta_{y_{j}})}{\sigma_{x}\sigma_{y}}$
$\displaystyle=\frac{\sum_{j=1}^{n_{s}}B_{j}(\sigma_{x_{j}}\sigma_{y_{j}}c_{j}+\delta_{x_{j}}\delta_{y_{j}})}{\sqrt{\sum_{i=1}^{n_{s}}B_{i}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}\sqrt{\sum_{i=1}^{n_{s}}B_{i}({\sigma_{y_{i}}}^{2}+{\delta_{y_{i}}^{2}})}}$

where $\textbf{x}^{j,i}(\textbf{y}^{j,i})$ is the
$\textbf{x}_{i}(\textbf{y}_{i})$ normalized in the $j$-th basic window. Now,
to show that
$\sigma_{x}=\sqrt{\frac{1}{T}{\sum_{1}^{n_{s}}B_{i}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}}$,
we evaluate:

$\displaystyle
T{\sigma_{x}}^{2}-{\sum_{1}^{n_{s}}B_{i}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}={\sum_{1}^{n_{s}}({\sigma_{x}}^{2}-{\sigma_{x_{i}}}^{2}-{\delta_{x_{i}}^{2}})}$
$\displaystyle=\sum_{k=1}^{n_{s}}\sum_{j=\Omega_{k-1}+1}^{\Omega_{k}}{\frac{1}{B_{k}}(\textbf{x}_{j}-\overline{x})^{2}-\frac{1}{B_{k}}(\textbf{x}_{j}-\overline{x_{k}})^{2}-\frac{1}{B_{k}}(\overline{x_{1}}-\overline{x_{k}})^{2}}$
$\displaystyle=\sum_{k=1}^{n_{s}}(\sum_{j=\Omega_{k-1}+1}^{\Omega_{k}}\frac{(-2\overline{x}{\textbf{x}_{j}}+2{\overline{x_{k}}}\,{\overline{x}})+(2{\textbf{x}_{j}}\overline{x_{k}}-2{\overline{x_{k}}}^{2})}{B_{k}})$
$\displaystyle=\sum_{k=1}^{n_{s}}\frac{1}{B_{k}}((-2B_{k}\overline{x}{\overline{x_{k}}}+2B_{k}{\overline{x_{k}}}\,{\overline{x}})+(2B_{k}{\overline{x_{k}}}^{2}-2B_{k}{\overline{x_{k}}}^{2}))=0$

∎

Using Lemma 1, we can pre-compute and store the statistics of basic windows of
once and compute the correlation of time-series for user-given time windows at
query time without performing a pass over time-series. Moreover, Lemma 1
allows us to support the arbitrary query window lengths. For instance, a user-
provided query window $x=[x_{i},\ldots,x_{j}]$ and $y=[y_{i},\ldots,y_{j}]$,
there exists a unique $\kappa\in\mathcal{N}$ such that $\kappa\cdot B\leq
i<(\kappa+1)\cdot B$, and there exists a unique $\chi\in\mathcal{N}$ such that
$\chi\cdot B<j\leq(\chi+1)\cdot B$. Let $B_{1}=(\kappa+1)\cdot B-i$,
$B_{n_{s}}=\chi\cdot B-j$, and $B_{k}$, for $k\in\\{2,\cdots,n_{s}-1\\}$. At
query time, we need to compute $\sigma_{x_{1}}(\sigma_{y_{1}}$),
$\sigma_{x_{n_{s}}}(\sigma_{y_{n_{s}}}$), $\delta_{x_{1}}(\delta_{y_{1}}$) and
$\delta_{x_{n_{s}}}(\delta_{y_{n_{s}}}$) from the raw data, and all the others
for the $B_{2},\cdots,B_{n_{s}-1}$ are pre-computed in the pre-processing.
Note that the case of equally subdividing time-series into basic windows of
size $B$ and a query window size being the integral multiple of the basic
window size is a special case of Lemma 1. For this special case, Algorithm 1
shows the steps of sketching basic windows and Algorithm 2 describes the steps
of constructing a network based on the exact correlation of time-series
calculated from the pre-computed statistics of basic windows.

####  3.1.2. Real-time Data

The correlation equation of Lemma 1 can be extended to deal with real-time
data. A user-defined query window on real-time data, $w=(``now",m)$, indicates
the sequence of the $m$ most recently observed data points of time-series.
That is, the size of the query window is fixed while the end timestamp is
changing as new data arrive. In our problem setting, since the data is
processed on the basis of basic windows, the algorithm waits until all new $B$
data points arrive. For time-series $x=[\textbf{x}_{1},\ldots,\textbf{x}_{m}]$
and $y=[\textbf{y}_{1},\ldots,\textbf{y}_{m}]$, and a query $w=(``now",m)$, we
can compute correlation at time $t$, namely $Corr_{t}(x,y)$, using Lemma 1.
This involves considering basic windows $[x_{1},\ldots,x_{n_{s}}]$ and
$[y_{1},\ldots,y_{n_{s}}]$, where $n_{s}=m/B$. At time $t+B$, the observed
time-series are $[\textbf{x}_{1},\ldots,\textbf{x}_{m+B}]$ and
$[\textbf{y}_{1},\ldots,\textbf{y}_{m+B}]$ and the basic windows are
$[x_{1},\ldots,x_{n_{s}+1}]$ and $[y_{1},\ldots,y_{n_{s}+1}]$. Based on query
$w=(``now",m)$, we need to consider $[x_{2},\ldots,x_{n_{s}+1}]$ and
$[y_{2},\ldots,y_{n_{s}+1}]$. According to Lemma 1, we can recalculate the
correlation at time $t+B$ from scratch. That is,

$\displaystyle
Corr_{t+B}(x,y)=\frac{\sum_{i=2}^{n_{s}+1}(\sigma_{x_{i}}\sigma_{y_{i}}c_{i}+\delta_{x_{i}}\delta_{y_{i}})}{\sqrt{\sum_{i=2}^{n_{s}+1}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}\sqrt{\sum_{i=2}^{n_{s}+1}({\sigma_{y_{i}}}^{2}+{\delta_{y_{i}}^{2}})}}$

Note that $\delta_{x_{j}}$’s and $\delta_{y_{j}}$’s have changed and needs to
be recalculated, since the means of the new query windows have probably
changed upon the arrival of new data. The following lemma allows us to compute
$Corr_{t+B}(x,y)$ by only using the statistics of the first and last basic
windows, without the need to calculate the statistics of the query window.

###### Lemma 0.

Given query windows $x=[\textbf{x}_{1},\ldots,\textbf{x}_{m}]$ and
$y=[\textbf{y}_{1},\ldots,\textbf{y}_{m}]$, basic windows
$[x_{1},\ldots,x_{n_{s}}]$ and $[y_{1},\ldots,y_{n_{s}}]$, and basic window
size $B=[B_{1},\ldots,B_{n_{s}}]$, where $T=\sum_{i=1}^{n_{s}}B_{i}$. Upon the
arrival of $B_{n_{s}+1}$ new data points, we have
$x=[\textbf{x}_{1},\ldots,\textbf{x}_{m+B_{n_{s}+1}}]$ and
$y=[\textbf{y}_{1},\ldots,\textbf{y}_{m+B_{n_{s}+1}}]$ and basic windows
$[x_{1},\ldots,x_{n_{s}+1}]$ and $[y_{1},\ldots,y_{n_{s}+1}]$. Let
$T^{\prime}=\sum_{i=2}^{n_{s}+1}B_{i}$. Considering a query window
$w=(``now",m)$, we can incrementally compute the Pearson’s correlation of $x$
and $y$ at time $t+B_{n_{s}+1}$ from their correlation at time $t$:

$\displaystyle Corr_{t+B_{n_{s}+1}}(x,y)=$
$\displaystyle\frac{1}{C~{}.~{}D}~{}\Big{(}T\sigma_{x}\sigma_{y}Corr_{t}(x,y)$
$\displaystyle+B_{n_{s}+1}(\sigma_{x_{n_{s}+1}}\sigma_{y_{n_{s}+1}}c_{n_{s}+1}+\delta_{x_{n_{s}+1}}\delta_{y_{n_{s}+1}})$
$\displaystyle-
B_{1}(\sigma_{x_{1}}\sigma_{y_{1}}c_{1}+\delta_{x_{1}}\delta_{y_{1}}){-}T^{\prime}\alpha_{x}\alpha_{y}\Big{)}$
$C=\sqrt{T\sigma_{x}^{2}+B_{n_{s}+1}({\sigma_{x_{n_{s}+1}}}^{2}+\delta_{x_{n_{s}+1}}^{2})-B_{1}(\sigma_{x_{1}}^{2}+\delta_{x_{1}}^{2})-T^{\prime}\alpha_{x}^{2}}$
$D=\sqrt{T\sigma_{y}^{2}+B_{n_{s}+1}({\sigma_{y_{n_{s}+1}}}^{2}+\delta_{y_{n_{s}+1}}^{2})-B_{1}(\sigma_{y_{1}}^{2}+\delta_{y_{1}}^{2})-T^{\prime}\alpha_{y}^{2}}$
$\alpha_{x}=\frac{B_{x_{n_{s}+1}}\delta_{n_{s}+1}-B_{1}\delta_{x_{1}}}{T}\quad
and\quad\alpha_{y}=\frac{B_{n_{s}+1}\delta_{y_{n_{s}+1}}-B_{1}\delta_{y_{1}}}{T}$
$\delta_{x_{n_{s}+1}}=\overline{x_{n_{s}+1}}-\overline{x_{1:n_{s}}}\quad
and\quad\delta_{y_{n_{s}+1}}=\overline{y_{n_{s}+1}}-\overline{y_{1:n_{s}}}$

where, $\sigma_{x}$ ($\sigma_{y}$) is the standard deviation of query window
$x$ ($y$) at time $t$, $\sigma_{x_{j}}$ ($\sigma_{y_{j}}$) is the standard
deviation of the basic window of $x_{j}$ ($y_{j}$), $c_{j}$ is the correlation
of the $j$-th basic windows of $x$ and $y$, $\overline{x_{j}}$
($\overline{y_{j}}$) is the mean of basic window $x_{j}$ ($y_{j}$), and
$\overline{x_{i:j}}$ ($\overline{y_{i:j}}$) is the mean of basic windows
$x_{i},\ldots,x_{j}$ ($y_{i},\ldots,y_{j}$).

###### Proof.

We denote $\sigma_{x}^{\prime}(\sigma_{y}^{\prime})$ to be the standard
deviation of the new query window $w$ after the arrival of the new basic
window, and $\delta_{x_{i}^{\prime}}(\delta_{y_{i}^{\prime}})$ to be
parameters of the new query window. We assume there is a linear transform
$Corr_{t+B_{n_{s}+1}}(x,y)=kCorr_{t}(x,y)+s$, where
$k=\frac{T\sigma_{x}\sigma_{y}}{T^{\prime}\sigma_{x}^{\prime}\sigma_{y}^{\prime}}$
and $s=\frac{1}{T^{\prime}\sigma_{x}^{\prime}\sigma_{y}^{\prime}}s^{\prime}$.
Looking at the change of the numerator, we have derived the following for
$s^{\prime}$:

$\displaystyle s^{\prime}$
$\displaystyle=T^{\prime}\sigma_{x}^{\prime}\sigma_{y}^{\prime}Corr_{t+B_{n_{s}+1}}(x,y)-T\sigma_{x}\sigma_{y}Corr_{t}(x,y)$
$\displaystyle=\sum_{i=2}^{n_{s}+1}B_{i}(\sigma_{x_{i}}\sigma_{y_{i}}c_{i}+\delta_{x_{i}^{\prime}}\delta_{y_{i}^{\prime}})-\sum_{i=1}^{n_{s}}B_{i}(\sigma_{x_{i}}\sigma_{y_{i}}c_{i}+\delta_{x_{i}}\delta_{y_{i}})$
$\displaystyle=B_{n_{s}+1}\sigma_{x_{n_{s}+1}}\sigma_{y_{n_{s}+1}}c_{n_{s}+1}-B_{1}\sigma_{x_{1}}\sigma_{y_{1}}c_{1}+\sum_{i=2}^{n_{s}+1}B_{i}\delta_{x_{i}^{\prime}}\delta_{y_{i}^{\prime}}$
$\displaystyle-\sum_{i=1}^{n_{s}}B_{i}\delta_{x_{i}}\delta_{y_{i}}$

Since
$\delta_{x_{i}}^{\prime}=\overline{x_{i}}-(\overline{x_{1:n_{s}}}+\alpha_{x})$
and
$\delta_{y_{i}}^{\prime}=\overline{y_{i}}-(\overline{y_{1:n_{s}}}+\alpha_{y})$,
we have

$\displaystyle\delta_{x_{i}}^{\prime}\delta_{y_{i}}^{\prime}$
$\displaystyle=(\overline{x_{i}}-(\overline{x_{1:n_{s}}}+\alpha_{x}))(\overline{y_{i}}-(\overline{y_{1:n_{s}}}+\alpha_{y}))$
$\displaystyle=(\delta_{x_{i}}-\alpha_{x})(\delta_{y_{i}}-\alpha_{y})$

Then, plugging it into the equation for $s^{\prime}$, we get

$\displaystyle s^{\prime}$
$\displaystyle=B_{n_{s}+1}(\sigma_{x_{n_{s}+1}}\sigma_{y_{n_{s}+1}}c_{n_{s}+1}+\delta_{x_{n_{s}+1}}\delta_{y_{n_{s}+1}})$
$\displaystyle-
B_{1}(\sigma_{x_{1}}\sigma_{y_{1}}c_{1}+\delta_{x_{1}}\delta_{y_{1}})-T^{\prime}\alpha_{x}\alpha_{y}$

To get $s$ and $k$, we also need the incremental equation for
$\sigma_{x}^{\prime}(\sigma_{y}^{\prime})$ if we look at the denominator. In
the proof of Lemma 1, we showed that
$\sigma_{x}=\sqrt{\frac{1}{T}{\sum_{1}^{n_{s}}B_{i}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}}$.
We have:

$\displaystyle T^{\prime}\sigma_{x}^{\prime 2}$
$\displaystyle=B_{2}{\sigma_{x_{2}}}^{2}+B_{2}{\delta_{x_{2}}^{\prime
2}}+\ldots+B_{n_{s}+1}\sigma_{{}_{n_{s}+1}}^{2}+B_{n_{s}+1}{\delta_{x_{n_{s}+1}}^{\prime
2}}$
$\displaystyle=\sum_{j=1}^{n_{s}}B_{j}{\sigma_{x_{j}}}^{2}+B_{n_{s}+1}{\sigma_{x_{n_{s}+1}}}^{2}-B_{1}{\sigma_{x_{1}}}^{2}+\sum_{j=2}^{n_{s}+1}B_{j}(\delta_{x_{j}}-\alpha_{x})^{2}$
$\displaystyle=\sum_{j=1}^{n_{s}}B_{j}{\sigma_{x_{j}}}^{2}+B_{n_{s}+1}{\sigma_{x_{n_{s}+1}}}^{2}-B_{1}{\sigma_{x_{1}}}^{2}$
$\displaystyle+\sum_{j=2}^{n_{s}+1}B_{j}({\delta_{x_{j}}}^{2}-2{\delta_{x_{j}}\alpha_{x}}+{\alpha_{x}}^{2})$
$\displaystyle=T\sigma_{x}^{2}+B_{n_{s}+1}({\sigma_{x_{n_{s}+1}}}^{2}+\delta_{x_{n_{s}+1}}^{2})-B_{1}(\sigma_{x_{1}}^{2}+\delta_{x_{1}}^{2})-T^{\prime}\alpha_{x}^{2}$

By replacing $k$ and $s$ in the above transform, the proof of the lemma
becomes complete. For real-time data, the mean of the query window changes as
new data arrives. For efficiency purposes, we do not want to compute that mean
for the calculation of parameters $\delta$. Lemma 1 allows us to compute the
standard deviation of a query window and correlation without computing the
mean. ∎

Algorithm 3 describes the steps of constructing a network for real-time data.

Algorithm 3 Network-Construct-RealTime

1:streams $\mathcal{L}=\\{\textbf{x}^{1},\ldots,\textbf{x}^{n}\\}$; statistics
$S$; query $w$; basic window size $B$; threshold $\theta$

2:graph $(G,V)$

3:$S\leftarrow{\sc Preprocessing}(\mathcal{L},B)$

4:$G,V\leftarrow{\sc Network-Construct-Histo}(\mathcal{L},S,w,B,\theta)$ //
create initial network

5:$b\leftarrow[]$ // most recent basic window

6:while do

7: $b\leftarrow{\bf IngestData}()$

8: if Len(b) == B then

9: $s\leftarrow{\bf Stats}(b)$

10: UpdateNetwork(G, V, s) // use Lemma 2

11: $b\leftarrow[]$

12:return

###  3.2. Approximate Pairwise Correlation

This considers non-aribtrary query window size. So far, we presented ways of
computing and updating the exact correlation of time-series. Now, we describe
how our model can be extended to approximate the correlation of time-series
over a query window for all time-series regardless of being cooperativeness or
uncooperative. Equation 3 shows how the DFT coefficients of two time series
can be reduced to the Euclidean distance of their normalized series, as
described in § 2.2 Note that, in our model, the necessary statistics for
normalization are collected during the sketch time.

####  3.2.1. Historical Data

Recall $d_{i}$ is the distance of the normalized $i$-th basic windows, namely
$\hat{x}_{i}$ and $\hat{y}_{i}$, $\hat{X}_{i}$ and $\hat{Y}_{i}$ are the DFT
of normalized basic windows $\hat{x}_{i}$ and $\hat{y}_{i}$, and
$Dist_{n}(\hat{X}_{i},\hat{Y}_{i})$ is the Euclidean distance of the first $n$
DFT coefficients in $\hat{X}_{i}$ and $\hat{Y}_{i}$. Since DFT preserves the
distance between coefficients and the original time-series, we have
$d_{i}\simeq Dist_{n}(\hat{X}_{i},\hat{Y}_{i})$. To compute the distance of
query windows, $Dist_{n}(x,y)$, from the distances of basic windows, without
any assumption about the form and properties of basic windows in a query
window, we can combine the equation of Lemma 1 and Equation 3 as follows.

$\displaystyle
1-\frac{1}{2}Dist_{n}(\hat{X},\hat{Y})^{2}\approx\frac{\sum_{i=1}^{n_{s}}(\sigma_{x_{i}}\sigma_{y_{i}}(1-\frac{d_{i}^{2}}{2})+\delta_{x_{i}}\delta_{y_{i}})}{\sqrt{\sum_{i=1}^{n_{s}}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}\sqrt{\sum_{i=1}^{n_{s}}({\sigma_{y_{i}}}^{2}+{\delta_{y_{i}}^{2}})}}$

We simplify the equation and obtain an approximation of the distance of two
query windows based on the distances of their basic windows.

(5) $Dist_{n}(\hat{X},\hat{Y})^{2}\approx
2+\frac{\sum_{i=1}^{n_{s}}\sigma_{x_{i}}\sigma_{y_{i}}d_{n}(\hat{X_{i}},\hat{Y_{i}})^{2}-2\sum_{i=1}^{n_{s}}(\sigma_{x_{i}}\sigma_{y_{i}}+\delta_{x_{i}}\delta_{y_{i}})}{\sqrt{\sum_{i=1}^{n_{s}}({\sigma_{x_{i}}}^{2}+{\delta_{x_{i}}^{2}})}\sqrt{\sum_{i=1}^{n_{s}}({\sigma_{y_{i}}}^{2}+{\delta_{y_{i}}^{2}})}}$

When all DFT coefficients are used, i.e. $n=B$, the $\approx$ becomes $=$,
turning into an exact calculation.

To perform all-pair correlation approximation in our framework, we can
normalize basic windows and compute their DFT coefficients, and pairwise
distances, during the sketch time (lines 8-10 of Algorithm 1). At query time,
we use Equation 5 to get $Dist_{n}(\hat{X},\hat{Y})$ and apply Equation 3 to
obtain the correlation. Algorithm 4 describes the steps of building a network
based on the approximation of correlation.

Algorithm 4 Network-Approximate

1:streams $\mathcal{L}=\\{\textbf{x}^{1},\ldots,\textbf{x}^{n}\\}$; statistics
$S$; query $w$; basic window size $B$; threshold $\theta$

2:graph $(G,V)$

3:$G\leftarrow\\{1,\ldots,n\\}$; $V\leftarrow\\{\\}$

4:$b\leftarrow{\bf GetBasicWins}(w)$ // basic window ids in $w$

5:for $x,y\in\mathcal{L}$ do

6: $S_{x}\leftarrow{\bf ReadStats}(S,b,x)$; $S_{y}\leftarrow{\bf
ReadStats}(S,b,y)$

7: $d_{1}\ldots d_{n_{s}}\leftarrow{\bf ReadStats}(S,x,y)$ // distance of
basic windows

8: if stats of basic windows $\simeq w$ then

9: $Dist\leftarrow{\bf Average}([d_{1},\ldots,d_{n_{s}}])$

10: else

11: $Dist\leftarrow{\bf Distance}(S_{x},S_{y},d_{1}..d_{n_{s}})$ // use
Equation 5

12: if $Dist\leq\sqrt{1-\theta}$ then

13: $V.{\bf Add}(x,y,Dist)$

14:return $(G,V)$

####  3.2.2. Real-time Data

Combining Equation 5 and Lemma 2, we can get the incremental update equation
for approximating pairwise correlation:

(6) $\begin{split}&2-Dist_{n}^{t+B}(\hat{X},\hat{Y})\underset{\text{= (when n
=b)}}{\approx}\\\
&\frac{1}{A~{}.~{}B}~{}\Big{(}n_{s}\sigma_{x}\sigma_{y}Dist_{n}^{t}(\hat{X},\hat{Y})+\sigma_{x_{n_{s}+1}}\sigma_{y_{n_{s}+1}}(1-\frac{d_{n_{s}+1})^{2}}{2})\\\
&-\sigma_{x_{1}}\sigma_{y_{1}}(1-\frac{d_{1}^{2}}{2})-\delta_{x_{1}}\delta_{y_{1}}{-}n_{s}\alpha_{x}\alpha_{y}{+}\delta_{x_{n_{s}+1}}\delta_{y_{n_{s}+1}}\Big{)}\end{split}$

Here, $Dist_{n}^{t+B}$ ($Dist_{n}^{t}(\hat{X},\hat{Y})$) is the DFT Distance
of the query window at time $t+B$ ($t$) using first $n$ coefficients in each
basic window. The new distance can be obtained by calculating the pairwise
distances for the last basic window $d_{n_{s}+1}$.

###  3.3. Complexity Analysis

In this section, we discuss the complexity analysis of query/sketch time and
space overhead of TSUBASA, the DFT-based algorithm, and the baseline algorithm
for non-arbitrary query windows. Next, we describe the synergies of time and
space with usability. Suppose $N$ is the number of time-series and each time-
series is in length $L$.

Space Complexity The space overhead of TSUBASA is
$\psi=\frac{L}{B}(2+\frac{N(N-1)}{2})$, where $B$ is the basic window size and
$\frac{L}{B}$ is the number of basic windows since we divide a time-series
evenly by default. For each basic window of a time-series, TSUBASA stores two
values for the mean and the standard deviation. In addition, for aligned basic
windows of all pairs of time-series, TSUBASA stores the correlation of each
pair of time-series. As a result, the space complexity of TSUBASA is
$O(\frac{LN^{2}}{B})$. The DFT-based approximate algorithm stores the mean and
the standard deviation for basic windows of each time-series and the distance
between the first few DFT coefficients of aligned basic windows of pairs of
time-series, thus, has the space complexity of $O(\frac{LN^{2}}{B})$. We
remark that this space overhead is in addition to the storage of raw time-
series for both algorithms if the raw time-series are not discarded after
sketching.

Time complexity The sketch time complexity of TSUBASA is independent of query
window size and is $O(L\cdot N^{2})$, since TSUBASA requires calculating
statistics over the aligned basic windows of all pairs of time-series. The
sketch time complexity of the approximate algorithm is worse than TSUBASA and
is $O(L^{2}\cdot N^{2})$, since the calculation of DFT coefficients for a
time-series of length $L$ is $O(L^{2})$ and coefficients are required for
calculating the distance of aligned basic windows in all pairs of time-series.
For a query window size $l^{*}=n_{s}\cdot B$, both TSUBASA and the approximate
algorithm scan all basic windows, therefore, the query time complexity of
TSUBASA and the approximate method are both $O(\frac{l^{*}}{B}\cdot N^{2})$.
However, the baseline algorithm scans the raw time-series and has the query
time complexity of $O(l^{*}\cdot N^{2})$.

The query time complexity of real-time TSUBASA is $O(B^{*}N^{2})$, where
$B^{*}$ is the size of the new coming basic window since TSUBASA needs to
compute statistics for the new window. The query time complexity of the real-
time approximate algorithm is $O(B^{*2}N^{2})$. The query time complexity of
real-time baseline algorithm is $O(L^{*}\cdot N^{2})$, where $L^{*}$ is the
size of the query window size.

Usability Discussion Let $M$ be the maximum space capacity available for the
storage of time-series sketches. Considering the above space analysis and
assuming equal-size basic windows, the minimum basic window size of TSUBASA
can be calculated by solving $\frac{L}{B}(2+\frac{N(N-1)}{2})\leq M$. That is,
with $M$ available storage the maximum basic window size handled by TSUBASA is
$\frac{L}{M}(2+\frac{N(N-1)}{2})$. Note that both time and space complexity
reduce as $B$ increases. Moreover, choosing a large $B$ means less space
capacity requirement. Therefore, should we just choose an extremely large $B$?
The answer is no. When an arbitrary query window is not supported, a large $B$
will reduce the flexibility of query windows, thus, usability. For the case of
the query window size being the integral multiple of the basic window size,
the chosen query window size by users becomes extremely limited. If we
consider the generic case of Lemma 1, we will observe a significant rise in
query time, since the start/end of a query window can fall anywhere in a basic
window, thus, when basic windows are large, the first and last basic windows
can be potentially large. Suppose the query window is in length $l^{*}$, where
$\exists n_{s}\in R$, such that $n_{s}\cdot B\leq l^{*}<(n_{s}+1)\cdot B$. The
time complexity is $O((\frac{l^{*}}{B}+B)\cdot N^{2})$. When $B>\sqrt{l}$,
$\frac{l^{*}}{B}+B$ is monotonically increasing. Since $B>\sqrt{l}$ at the
most meaningful queries, the query time increases when the $B$ increases for
the generic method.

###  3.4. Parallel and Disk-based TSUBASA

The disk-based TSUBASA stores sketches on the disk to be retrieved at query
time for correlation calculation. Moreover, despite the quadratic complexity
of the sketch time and query time, TSUBASA is embarrassingly parallelizable.
The set of all pairs of time-series can be partitioned into groups that are
processed in parallel. During sketching, workers are divided into a database
worker, that writes statistics to the database, and computation workers, that
perform sketch computation. Each worker sketches time-series pairs of a
partition and sends the sketches in batches to the database worker to write to
a disk-based database. During the query time, each worker is assigned a
partition, reads the sketches of time-series in batches directly from the
database and computes the pairwise correlations, and outputs a sub-matrix of
the correlation matrix.

To leverage data locality and minimize the number I/Os, for partitioning time-
series pairs, TSUBASA adopts an approach similar to the parallel block nested
loop join. Each partition contains a subset of time-series paired with all
time-series. i.e. each partition is a group of rows in a correlation matrix
and the processing is done row by row in batches. Batches of pairs are
assigned to a worker and once a worker is finished, it reads the statistics of
the next batch of pairs from the database. Since Pearson’s correlation is a
symmetric measure, TSUBASA needs to process $n(n-1)/2$ pairs to construct the
correlation matrix. For load balancing, TSUBASA assigns the same number of
pairs to each worker. Note that the same architecture can be used to make the
machinery described, in § 3.2, for correlation approximation.

###  3.5. Solutions to threshold correlation matrix

Suppose that we know $c_{xz}$ and $c_{yz}$, and we want to infer the range of
$c_{xz}$. Referring to (glass1970geometric), we have the upper bound and the
lower bound of $c_{xz}$:

(7) $\begin{split}c_{xz}c_{yz}-\sqrt{(1-c_{xz}^{2})(1-c_{yz}^{2})}\leq
c_{xy}\leq c_{xz}c_{yz}+\sqrt{(1-c_{xz}^{2})(1-c_{yz}^{2})}\end{split}$

The equation above provides us a way in predicting the whole correlation
matrix based on a small amount of correlations. For example, we are given $N$
time series and $c_{12},\cdots,c_{1N}$. We could make prediction on $c_{ij}$
for any $1\leq i,j\leq N$ from $c_{1i}$ and $c_{1j}$ based on Equation 7. For
example, given a positive threshold $\theta$, let $U_{xy:z}$ be the upper
bound of $c_{xy}$ estimated from $c_{xz}$ and $c_{yz}$, and $L_{xy:z}$ be the
lower bound of $c_{xy}$ estimated from $c_{xz}$ and $c_{yz}$. If
$L_{xy:z}\geq\theta$ or $U_{xy:z}\leq-\theta$, then we know that $m_{xy}=1$,
where $m_{ij}$ is the $i$ th row and $j$ th column of the output matrix. If
$L_{xy:z}\geq-\theta$ and $U_{xy:z}\leq\theta$, then we know that $m_{xy}=0$.

Algorithm 5 Pruning

1:streams $\mathcal{L}=\\{\textbf{x}^{1},\ldots,\textbf{x}^{N}\\}$;threshold
$\theta$

2:Matrix $M(m_{ij})$

3:$m_{ij}=-\infty$

4:for $i=1,2,\cdots,N$ do

5: if $\exists p,q\textbf{s.t.}{m_{pq}}<0$ then

6: $c_{i1},\cdots,c_{iN}\leftarrow\textbf{Computecorr}(\mathcal{L},i)$

7: for $j=1,2,\cdots,N$ do

8: for $k=1,2,\cdots,N$ do

9: $L_{jk:i},U_{jk:i}\leftarrow\textbf{Correct-Inference}(c_{ij},c_{ik})$

10: if $L_{jk:i}\geq\theta$ or $U_{jk:i}\leq-\theta$ then

11: $m_{jk}\leftarrow 1$

12: if $L_{jk:i}\geq-\theta$ and $U_{jk:i}\leq\theta$ then

13: $m_{jk}\leftarrow 0$

14:$M\leftarrow\textbf{Compute-Rest}$($M$)

15:return $M$

Figure 4. Correlation Inference

In the Figure 4, let the horizontal axis be the $c_{xz}$, and the vertical
axis be the $c_{yz}$. If $\theta$ is given to be $0.8$, then the colored
regions are the cases that we could know the $m_{xy}$ without computing
$c_{xy}$. In the blue regions, $m_{xy}=1$, and $m_{xy}=0$ in the red regions.
The white region is the uncertainty part, we would need to change a $z$
(anchor) to see whether we could infer/compute further. The algorithm 5
presents how we can use it to prune. We select the anchor randomly (or let
every time-series to be the anchor), then if we use the selected anchor to
scan the $M$ and see which $m_{jk}$ we could know by the bounds instead of
computing $c_{jk}$.

##  4\. Experiments

We have developed, in this paper, mathematical models and algorithms: Network-
Construct-Histo and Network-Construct-RealTime, for constructing and updating
correlation matrices to build exact networks on historical and real-time data.
Our empirical evaluation has two parts. First, we study these algorithms and
compare their query time and sketch time against a baseline, on historical and
real-time version of a climate data set. For these experiments, we use the in-
memory version of the algorithms, i.e. in-memory data structures are used for
storing raw data and sketches. Second, we evaluate the scalability and
efficiency of the disk-based and parallel TSUBASA and the approximate
algorithm as described in § 3.4. For all experiments, we assume equal basic
window sizes in time-series. All algorithms are implemented using Go language.
We use PostgreSQL for storing data sketches. All experiments are conducted on
a machine with 2 Intel® Xeon Gold 5218 @ 2.30GHz (64 cores), 512 GB DDR4
memory, a Samsung® SSD 983 DCT M.2 (2 TB).

NCEA Data
Set111https://www.ncei.noaa.gov/pub/data/uscrn/products/hourly02/2020/ is a
public data from the National Oceanic and Atmospheric Administration (NOAA).
The data is collected every hour, and uploaded publicly in 24-hour increments.
NOAA utilizes radiometric satellite collection, buoys, weather stations,
citizen scientists, and other methods for perpetual data gathering. The data
is collected from 157 nodes (time-series) across the US. Each node produces
approximately 8,760 points of data in a year. This data set is used for in-
memory experiments.

Berkeley Earth Data Set222http://berkeleyearth.org/data/ is a collection of
open-source data sets provided by an independent U.S. non-profit organization
(Berkeley Earth). We use NetCDF-format gridded data from this data set. The
climate data includes average temperature data on both lands and oceans. It
divides the earth by 1∘ $\times$ 1∘ latitude-longitude grid. We consider the
land time-series in this data set. The data set includes 18,638 nodes and each
nodes contains 3,652 data points. The time resolution is 24 hours. This data
set is used for scalability experiments.

###  4.1. Accuracy

We compared the accuracy of the climate network of NCEA data set, constructed
based on the correlation matrix computed by the DFT-based techniques (Zhu and
Shasha, 2002; Mueen et al., 2010) (as described in § LABEL:sec:bgapprox) and
exact calculation, followed by the application of a threshold. The approximate
technique (Zhu and Shasha, 2002; Cole et al., 2005) uses the first few DFT
coefficients for estimating the distance of aligned basic windows, then, basic
window distances are aggregated to obtain an approximation of the distance and
correlation of time-series on a query window (Zhu and Shasha, 2002). In our
experiments, we use the way, we believe, StatStream (Zhu and Shasha, 2002)
computes the distance (correlation) of query windows i.e. by averaging the
distance (correlation) of DFT coefficients over all basic windows.

We evaluate the impact of approximation on the accuracy of constructed
networks, using two measures: number of edges and the correlation similarity
ratio, inspired by (Napoletani and Sauer, 2008). A correlation matrix is an
$n\times n$ matrix, where $n$ is the number of time-series and a cell $c_{ij}$
is a binary value that indicates the correlation score of time-series $x^{i}$
and $x^{j}$ is higher than threshold $\theta$. The correlation similarity
ratio evaluates the percentage of identical edges in two networks. Formally,
given two complex networks represented by adjacency matrices $A:\\{a_{ij}\mid
0\leqslant i,j\leqslant n\\}$ and $B:\\{b_{ij}\mid 0\leqslant i,j\leqslant
n\\}$, the similarity ratio is defined as follows.

$D_{p}(A,B)=\frac{2\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}1-|a_{ij}-b_{ij}|}{n(n-1)}$

For instance, the correlation similarity ratio of networks with the adjacency
matrices $A$ and $B$ is $2/3$.

$A=\begin{pmatrix}1&1&0\\\ 1&1&1\\\ 0&1&1\\\ \end{pmatrix}\quad
B=\begin{pmatrix}1&0&0\\\ 0&1&1\\\ 0&1&1\\\ \end{pmatrix}\quad
D_{p}(A,B)=\frac{2}{3}$

(a)

(b)

(c)

(d)

Figure 5. In-memory: (a) Network Accuracy Comparison (b) Basic Window Size
Analysis (c) Query Window Size Analysis (d) Network Update Time.

For both techniques we consider the basic window size 200 and threshold 0.75,
while varying the number of DFT coefficients from 50 to 200 for the
approximate technique. Note that in the exact technique (basic window
correlation), the correlation of time-series are computed by aggregating the
correlation of basic windows as suggested by (Zhu and Shasha, 2002).
Therefore, the structure of this network (the solid red plot) is independent
of the number of DFT coefficients.

As shown in Figure 5(a), the number of edges in the network constructed by DFT
correlation calculation becomes equal to the number of edges in the network
constructed by exact calculation, only when all 200 coefficients are used.
This matches the theory, i.e. the approximation becomes identical to the exact
calculation, when all DFT coefficients are used. Note that the approximate
technique uses Equation 4 to find correlated time-series based on their DFT-
based distance. Following this rule, the DFT correlation calculation never
yields false negatives, but, creates false positive edges. This explains why
the number of edges in networks decrease as more coefficients are used.
Moreover, the similarity ratio of correlation matrices increases as the number
of considered coefficients increases and is at its highest value when all
coefficients are used to represent a basic window.

The main take-away is that constructing a network based on the approximation
of DFT-based distance can lead to inaccurate networks. For climate data sets,
near exact result is obtained only when a very large number of coefficients
are used for approximation. This means smaller basic windows are preferred for
approximation purposes which leads to higher number of basic windows,
therefore, higher correlation calculation time in addition to the high DFT
coefficient calculation time. These results highlight the necessity of
efficient algorithms for constructing and updating exact networks on large
collections of time-series.

###  4.2. Efficiency

We evaluate the efficiency of the in-memory version of correlation matrix
calculation algorithms with respect to query window size and basic window size
parameters.

Network Construction We compare the sketch time plus query time when using the
DFT-based approximation of StatStream with TSUBASA’s exact correlations. For
the approximation technique, we report on two scenarios: using all DFT
coefficients and using 75% of coefficients of a basic window. As shown in §
4.1, the former empirically yields a network similar to the network of exact
correlation calculation. During the sketch time, TSUBASA calculates the
statistics from Lemma 1 and the approximation algorithm calculates the
statistics necessary for Equation 5 for all basic windows of all time-series.
At query time, Lemma 1 and Equation 5 are used to combine sketched statistics
to get approximate and exact networks, respectively.

Figure 5(b) reports the run time when varying the size of basic window for a
query window of size 3,000. The sketch time of TSUBASA grows very gradually
with the basic window size, while the sketch time of the approximate algorithm
increases with the size of basic window. This is because of the $O(n^{2})$
complexity of DFT calculation. Our results show that TSUBASA outperforms the
approximation technique at sketch time and its query time is on par with the
approximate network construction technique.

Figure 5(c) shows the query time of TSUBASA, approximate calculation, and a
baseline, when varying the query window size, considering the constant basic
window size of 50. The baseline algorithm computes the Pearson’s correlation
of Equation 1 for all pairs of time-series directly from raw data at query
time without any sketching. In this experiment the approximate algorithm uses
75% of DFT coefficients of a basic window. Note that the distances of basic
windows ($d_{j}$’s in Equation 5) are calculated during the sketch time,
therefore, the query time of the approximate algorithm does not depend on the
number of considered DFT coefficients. TSUBASA is almost as fast as the
approximate algorithm for all query window sizes and outperforms the baseline
by two orders of magnitude. We remark that all algorithms have quadratic
complexity in the number of time-series. However, the exact and DFT-based
approximation are extremely efficient at computing correlation of each pair at
query time due to relying on statistics that are pre-calculated during the
sketch time.

Network Update We compare the network update time of TSUBASA with the DFT-
based approximation of § 3.2.2 for real-time NCEA data set. The initial
networks are constructed on the data set for a given query window. Then, after
the arrival of $B$ data points, both algorithms update the correlation and
network using the special case of Lemma 2 and Equation 6, respectively. Figure
5(d) shows the time taken by TSUBASA and approximate algorithm upon the
arrival of $B$ new data points for various basic window sizes for a query
window size of 3,000. The approximate algorithm uses 75% of DFT coefficients
in a basic window. For both algorithms, updates to the network depend on the
statistics of the first basic window of the current query window, which is
already calculated, and the most recently observed basic window, which needs
to be calculated on the fly. Indeed, the efficiency of update only depends on
the processing of the most recently observed basic window. Since the
approximation algorithm needs to calculate DFT coefficients we observe that it
is slower than TSUBASA at least one order of magnitude. The gap between the
two algorithms becomes more obvious for larger basic window sizes because of
the $O(n^{2})$ complexity of DFT calculation. In conclusion, TSUBASA can
compute exact correlation and networks for real-time much faster than the
approximation competitor.

(a)

(b)

(c)

(d)

Figure 6. Parallel and Disk-based: (a) Sketch Time Breakdown (b) Query Time
Breakdown (c) Impact of Parallelization (d) Space Overhead.

###  4.3. Scalability

We compare TSUBASA and the approximation algorithm in similar parallel and
disk-based configurations. To separate the impact of fine-tuning the database
on performance, in all experiments, we choose to use one database worker and
allocate the rest of workers for sketching and querying. For the scalability
experiments, we use subsets of time-series from the Berkeley Earth data set.
All experiments consider a basic window length of 120, a query window length
of 960, and 75% of DFT coefficients for correlation approximation.

Sketch Time Figure 6(a) shows the sketch time of TSUBASA and the approximate
algorithm for correlation matrix calculation for various number of time-series
on 63 partitions and 64 cores. The plot separates the write time from the
sketch calculation time. We observe that TSUBASA outperforms the approximate
algorithm in sketch time. This is due to quadratic complexity of DFT
calculation as opposed the linear complexity of computing TSUBASA sketches. We
observe that the majority of work by TSUBASA during sketching is spent on
writing sketches to a database, unlike matrix approximation which is on par
with the write time. Note that in this configuration the total sketch time of
TSUBASA and the approximate algorithm is bounded by database write time. The
total sketch time, sketch calculation, and write time of TSUBASA and the
approximate algorithm increase quadratically with the number of time-series.
However, due to parallelization, the growth is slower than what is expected
for a single-core configuration.

Query Time Figure 6(b) shows the query time of TSUBASA and the approximate
algorithm for correlation matrix calculation for various number of time-series
on 63 partitions and 64 cores. The plot separates the database read time from
the correlation matrix calculation time. Both TSUBASA and the approximate
algorithm have on par query time and take less than a minute for computing the
correlation matrix even for the largest number of time-series. We observe that
the read time during querying is negligible compared to matrix calculation.
The read time percentage is slightly higher for smaller networks due to the
database overhead compared to matrix calculation cost on small number of time-
series. The total query time, matrix calculation, and read time of TSUBASA and
the approximate algorithm increase quadratically with the number of time-
series. However, due to parallelization, the growth is slower than what is
expected for a single-core configuration.

Impact of Number of Partitions Figure 6(c) shows how TSUBASA scales with the
number of partitions. For these experiments, we use 2,000 time-series. Note
than we have 63 partitions/cores for sketch/matrix computation and we reserve
one core for database writes. Both sketch and matrix calculation times
decrease with the increase in the number of cores. We expect further
optimization of query and sketch time can be done by fine-tuning the database
and allocating further resources.

Space Overhead Figure 6(d) shows the size of databases used for storing
sketches of 2,000 time-series by TSUBASA and the approximate algorithm with
respect to basic window size. Both algorithms store sketches of the the same
size for each basic window and have the same space overhead. As the size of
basic window increase the number of basic windows decreases and total size of
sketches stored by both algorithms decrease.

##  5\. Related Work

Spatio-Temporal Databases represent the value of a climate variable with three
dimensions of geometry and time (i.e. latitude, longitude, and timestamp).
Systems such as Microsoft StreamInsight (Ali et al., 2010), GeoMesa (geo,
[n.d.]), and IBM PAIRS Geoscope (Klein et al., 2015) are designed for
processing streams of geospatial data, from sources such as satellites and IoT
sensors, benefiting from relational DBMSs, distributed column-oriented
databases, and scalable key-value data stores. The algorithms and mathematical
models we designed in this paper can be incorporated into geo-spatial systems.
We use the progressive and declarative processing of Trill (Chandramouli et
al., 2015) for storage and analysis. TSUBASA can be a stand-alone system with
other network analysis extensions for clustering and community detection.

Similarity Search on Time-Series To compute the similarity between time-
series, several measures have been proposed (Mirylenka et al., 2017;
Paparrizos et al., 2020). In TSUBASA, we consider the Pearson’s correlation
coefficient as it is the most commonly used measure for building climate
networks (Faghmous and Kumar, 2014; Tsonis and Roebber, 2004). There has been
extensive work from database community on similarity search of time-series
(Cai et al., 2021; Gogolou et al., 2020; Mirylenka et al., 2017; Zhu and
Shasha, 2002; Qiu et al., 2018). This line of work considers threshold queries
for similarity search in time-series databases and often involve calculating a
function over two or more streams and reporting when the threshold is crossed.
Time-series similarity search problems in the community have mostly focused on
identifying similar time-series to a query time-series. TSUBASA, however,
focuses on the construction of the complete and exact correlation matrix, a
task that requires all-pair correlation calculation.

Sketching and Data Reduction Alternative techniques to DFT for time-series
similarity approximation are Discrete Wavelet Transform (DWT) (Chaovalit et
al., 2011), Singular Value Decomposition (Jang et al., 2018), and Piecewise
Constant Approximation (Luo et al., 2015). Data reduction is based on the idea
of summarizing the population or sample data through smaller-sized matrices or
simple numbers (Khattree and Naik, 2000). As for time-series, it has been a
topic of interest that reduces data into low-dimensional data while preserving
its characteristics to a large extent (Zhao, 2006), which has IoT applications
(Papageorgiou et al., 2015).  TSUBASA sketches data into statistics that are
required for the efficient and exact calculation of correlation scores and
networks.

Data Streaming Systems Most data streaming systems, including S4 (Neumeyer et
al., 2010), Muppet (Lam et al., 2012), Spark (Zaharia et al., 2016), and Flink
(Carbone et al., 2015) are large-scale data processing systems driven by th
Reduce programming model. They can work on different kinds of data coming from
real-world sensors or IoT devices. The system we used in the experiments,
Trill, also has a distributed version called Quill (Chandramouli et al.,
2016). The distributed design will definitely reduce the latency of data
processing, and also increase the throughput in a given time.

##  6\. Discussion and Conclusion

We presented TSUBASA, an efficient and exact correlation matrix calculation
algorithm for climate network construction on historical and real-time data.
TSUBASA uses the basic-window model to subdivide a query time-window into
smaller windows. TSUBASA computes a cheap and simple sketch of basic windows
and reuses them at query time for building networks on arbitrary query
windows. We describe a way of approximating time-series correlation and
compared it with TSUBASA. Experiments show that TSUBASA can compute exact
correlation and network faster than DFT-based approximation techniques. The
techniques proposed in TSUBASA can be potentially applied to analyzing stock
market data (Memon and Yao, 2019) as well as biological data (Batushansky et
al., 2016). However, our preliminary investigation shows that approximate and
partial calculation of pairwise correlations suffices in scenarios such as
stock market data (Memon and Yao, 2019). For future work, we plan to extend
our problem definition to unaligned time-series, develop a pairwise
correlation pruning algorithm based on a threshold, and consider further
optimization of the parallel TSUBASA by fine-tuning.

## References

  * (1)
  * geo ([n.d.]) [n.d.].  GeoMesa.  https://github.com/locationtech/geomesa.
  * Abe and Suzuki (2004) Sumiyoshi Abe and Norikazu Suzuki. 2004.  Scale-free network of earthquakes.  _EPL (Europhysics Letters)_ 65, 4 (2004), 581\.
  * Ali et al. (2010) Mohamed H. Ali, Badrish Chandramouli, Balan Sethu Raman, and Ed Katibah. 2010.  Spatio-Temporal Stream Processing in Microsoft StreamInsight.  _IEEE Data Eng. Bull._ 33, 2 (2010), 69–74.
  * Batushansky et al. (2016) Albert Batushansky, David Toubiana, and Aaron Fait. 2016\.  Correlation-based network generation, visualization, and analysis as a powerful tool in biological studies: a case study in cancer cell metabolism.  _BioMed research international_ 2016 (2016).
  * Berezin et al. (2012) Y. Berezin, A. Gozolchiani, O. Guez, and S. Havlin. 2012\.  Stability of Climate Networks with Time.  _Scientific Reports_ 2 (2012).
  * Cai et al. (2021) Walter Cai, Philip A. Bernstein, Wentao Wu, and Badrish Chandramouli. 2021.  Optimization of Threshold Functions over Streams.  _PVLDB_ 14, 6 (2021), 878–889.
  * Carbone et al. (2015) Paris Carbone, Asterios Katsifodimos, Stephan Ewen, Volker Markl, Seif Haridi, and Kostas Tzoumas. 2015.  Apache Flink™: Stream and Batch Processing in a Single Engine.  _IEEE Data Eng. Bull._ 38, 4 (2015), 28–38.
  * Chandramouli et al. (2016) Badrish Chandramouli, Raul Castro Fernandez, Jonathan Goldstein, Ahmed Eldawy, and Abdul Quamar. 2016.  Quill: Efficient, Transferable, and Rich Analytics at Scale.  _PVLDB_ 9, 14 (2016), 1623–1634.
  * Chandramouli et al. (2015) B. Chandramouli, J. Goldstein, M. Barnett, and J. F. Terwilliger. 2015.  Trill: Engineering a Library for Diverse Analytics.  _IEEE Data Eng. Bull._ 38, 4 (2015), 51–60.
  * Chaovalit et al. (2011) Pimwadee Chaovalit, Aryya Gangopadhyay, George Karabatis, and Zhiyuan Chen. 2011.  Discrete wavelet transform-based time series analysis and mining.  _ACM Comput. Surv._ 43, 2 (2011), 6:1–6:37.
  * Cole et al. (2005) Richard Cole, Dennis E. Shasha, and Xiaojian Zhao. 2005\.  Fast Window Correlations over Uncooperative Time Series. In _SIGKDD_. 743–749.
  * Donges et al. (2009) Jonathan F Donges, Yong Zou, Norbert Marwan, and Jürgen Kurths. 2009.  Complex networks in climate dynamics.  _The European Physical Journal Special Topics_ 174, 1 (2009), 157–179.
  * Dunlap (1937) Jack W Dunlap. 1937\.  Combinative properties of correlation coefficients.  _The Journal of Experimental Education_ 5, 3 (1937), 286–288.
  * Faghmous and Kumar (2014) James H. Faghmous and Vipin Kumar. 2014.  A Big Data Guide to Understanding Climate Change: The Case for Theory-Guided Data Science.  _Big Data_ 2, 3 (2014), 155–163.
  * Gogolou et al. (2020) Anna Gogolou, Theophanis Tsandilas, Karima Echihabi, Anastasia Bezerianos, and Themis Palpanas. 2020\.  Data Series Progressive Similarity Search with Probabilistic Quality Guarantees. In _SIGMOD_. 1857–1873.
  * Gozolchiani et al. (2008) Avi Gozolchiani, Kazuko Yamasaki, Oz Gazit, and Shlomo Havlin. 2008.  Pattern of climate network blinking links follows El Niño events.  _EPL (Europhysics Letters)_ 83, 2 (2008), 28005\.
  * Jang et al. (2018) Jun-Gi Jang, Dongjin Choi, Jinhong Jung, and U Kang. 2018\.  Zoom-SVD: Fast and Memory Efficient Method for Extracting Key Patterns in an Arbitrary Time Range. In _CIKM_. 1083–1092.
  * Khattree and Naik (2000) Ravindra Khattree and Dayanand N Naik. 2000.  Multivariate data reduction and discrimination.  _SAS Institute, Cary, North Carolina_ (2000).
  * Kim et al. (2019) Kyunghun Kim, Hongjun Joo, Daegun Han, Soojun Kim, Taewoo Lee, and Hung Soo Kim. 2019\.  On complex network construction of rain gauge stations considering nonlinearity of observed daily rainfall data.  _Water_ 11, 8 (2019), 1578.
  * Klein et al. (2015) Levente J Klein, Fernando J Marianno, Conrad M Albrecht, Marcus Freitag, Siyuan Lu, Nigel Hinds, Xiaoyan Shao, Sergio Bermudez Rodriguez, and Hendrik F Hamann. 2015.  PAIRS: A scalable geo-spatial data analytics platform. In _Big Data_. 1290–1298.
  * Lam et al. (2012) Wang Lam, Lu Liu, STS Prasad, Anand Rajaraman, Zoheb Vacheri, and AnHai Doan. 2012\.  Muppet: MapReduce-Style Processing of Fast Data.  _Proc. VLDB Endow._ 5, 12 (2012), 1814–1825.
  * Luo et al. (2015) Ge Luo, Ke Yi, Siu-Wing Cheng, Zhenguo Li, Wei Fan, Cheng He, and Yadong Mu. 2015.  Piecewise linear approximation of streaming time series data with max-error guarantees. In _ICDE_. 173–184.
  * Memon and Yao (2019) Bilal Ahmed Memon and Hongxing Yao. 2019.  Structural change and dynamics of Pakistan stock market during crisis: A complex network perspective.  _Entropy_ 21, 3 (2019), 248.
  * Mirylenka et al. (2017) Katsiaryna Mirylenka, Michele Dallachiesa, and Themis Palpanas. 2017.  Data Series Similarity Using Correlation-Aware Measures. In _SSDBM_. 11:1–11:12.
  * Mueen et al. (2010) Abdullah Mueen, Suman Nath, and Jie Liu. 2010.  Fast approximate correlation for massive time-series data. In _SIGMOD_. 171–182.
  * Napoletani and Sauer (2008) Domenico Napoletani and Timothy D Sauer. 2008.  Reconstructing the topology of sparsely connected dynamical networks.  _Physical Review E_ 77, 2 (2008), 026103.
  * Neumeyer et al. (2010) Leonardo Neumeyer, Bruce Robbins, Anish Nair, and Anand Kesari. 2010.  S4: Distributed Stream Computing Platform. In _ICDMW 2010, The 10th IEEE International Conference on Data Mining Workshops, Sydney, Australia, 13 December 2010_. IEEE Computer Society, 170–177.
  * Nocke et al. (2015) Thomas Nocke, Stefan Buschmann, Jonathan Friedemann Donges, Norbert Marwan, H-J Schulz, and Christian Tominski. 2015.  visual analytics of climate networks.  _Nonlinear Processes in Geophysics_ 22, 5 (2015), 545–570.
  * Papageorgiou et al. (2015) Apostolos Papageorgiou, Bin Cheng, and Ernö Kovacs. 2015\.  Real-time data reduction at the network edge of Internet-of-Things systems. In _2015 11th international conference on network and service management (CNSM)_. IEEE, 284–291.
  * Paparrizos et al. (2020) John Paparrizos, Chunwei Liu, Aaron J. Elmore, and Michael J. Franklin. 2020.  Debunking Four Long-Standing Misconceptions of Time-Series Distance Measures. In _SIGMOD_. 1887–1905.
  * Qiu et al. (2018) Han Qiu, Hoang Thanh Lam, Francesco Fusco, and Mathieu Sinn. 2018\.  Learning Correlation Space for Time Series.  arXiv:1802.03628 [cs.LG]
  * Rafiei (1999) Davood Rafiei. 1999\.  On Similarity-Based Queries for Time Series Data. In _ICDE_. 410–417.
  * Tantet and Dijkstra (2014) Alexis Tantet and Henk A Dijkstra. 2014.  An interaction network perspective on the relation between patterns of sea surface temperature variability and global mean surface temperature.  _Earth System Dynamics_ 5, 1 (2014), 1–14.
  * Tsonis and Roebber (2004) A. A. Tsonis and P. J. Roebber. 2004.  The architecture of the climate network.  _Physica A_ 333 (Feb. 2004), 497–504.
  * Zaharia et al. (2016) Matei Zaharia, Reynold S. Xin, Patrick Wendell, Tathagata Das, Michael Armbrust, Ankur Dave, Xiangrui Meng, Josh Rosen, Shivaram Venkataraman, Michael J. Franklin, Ali Ghodsi, Joseph Gonzalez, Scott Shenker, and Ion Stoica. 2016\.  Apache Spark: a unified engine for big data processing.  _Commun. ACM_ 59, 11 (2016), 56–65.
  * Zhao (2006) Xiaojian Zhao. 2006\.  _High performance algorithms for multiple streaming time series_.  New York University.
  * Zhu and Shasha (2002) Yunyue Zhu and Dennis E. Shasha. 2002.  StatStream: Statistical Monitoring of Thousands of Data Streams in Real Time. In _VLDB_. 358–369.