Dynamic Multi-Network Mining of Tensor Time Series

Fig. 1 shows the DMM results for clustering over Google Trends data, which consists of $10$ years of daily web search counts for six queries related to COVID-19 across $10$ countries, forming a $3^{rd}$-order tensor. Fig. 1 (a) shows the cluster assignments of the TTS, where each color represents a cluster. DMM splits the tensor into four segments and groups them into four clusters, each of which can be interpreted as a distinct phase corresponding to the evolving social response to COVID-19; thus, we name these phases “Before Covid,” “Outbreak,” “Vaccine,” and “Adaptation.” It is worth noting that this result is obtained with no prior knowledge.

Fig. 1 (b) presents the networks of each cluster, i.e., a country network, which has nodes plotted on the world map, reflects dependencies between different countries, and a query network for query dependencies. These networks, also known as a Markov Random Field (MRF) (Rue and Held, 2005), illustrate how the node affects the other nodes. The thickness and color of the edges in the network indicate the strength of the partial correlation between the nodes, which denotes a stronger relationship compared with a simple correlation. We learn the networks by estimating a Gaussian inverse covariance matrix. Then, by definition, if there is an edge between two nodes, the nodes are directly dependent on each other. Otherwise, they are conditionally independent, given the rest of the nodes. Moreover, we impose an $\ell_{1}$-norm penalty on the networks to promote sparsity, making it possible to obtain true networks and interpretability, as well as making the method noise-robust (Wytock and Kolter, 2013; Yuan and Lin, 2006). These networks provide visible and interpretable insights into the key relationships that characterize clusters.

We see that each of the four clusters exhibits unique networks that evolve with the different phases. In the “Before Covid” phase, the country network displays edges between English-speaking countries, indicating their interconnectedness. In the query network, the query “vaccine” correlates with “influenza.” However, during the “Outbreak” starting in $2020$, many countries respond to the COVID-19 pandemic, leading to various edges in the country network. In the query network of this phase, new edges related to “coronavirus” appear, and “coronavirus” and “virus” have a particularly strong connection. In the “Vaccine” phase, as people become more concerned about protection from COVID-19, the query “vaccine” forms an edge with “covid.” Moreover, since flu infects fewer people than in the past, “influenza” loses its edges. Lastly, during the “Adaptation” phase, as the world becomes accustomed to the situation, the country network reduces the number of edges, and the edges related to “influenza” reappear, reflecting a return to the networks observed in the “Before Covid” phase.

In summary, we propose DMM as a subsequence clustering method for TTS based on the MDL principle that enables each cluster to be characterized by multiple networks. The contributions of this paper can be summarized as follows.

• •

  Interpretable: DMM realizes the meaningful subsequence clustering of TTS, where each cluster is characterized by sparse dependency networks for each non-temporal mode, which facilitates the interpretation of the cluster from important relationships between variables.

• •

  Accurate: We define a criterion based on MDL to discover clusters with distinct networks. Thanks to the proposed criterion, any user-defined parameters can be determined, and DMM outperforms its state-of-the-art competitors in terms of clustering accuracy on synthetic data.

• •

  Scalable: The proposed clustering algorithm in DMM scales linearly as regards the input data size and is thus applicable to long-range and high-dimensional tensors.

We briefly introduce some definitions in tensor algebra from tensor related literature (Kolda and Bader, 2009; Cai et al., 2015).

Given a tensor $\mathcal{X}\in\mathbf{R}^{D^{(1)}_{1}\times\cdots\times D^{(1)}_{N}\times D^{(2)}_{1}\times\cdots\times D^{(G)}_{N}}$, we partition the modes into $G$, $P^{(g)}=\{gN+1,\cdots,g(N+1)\}$. The element is given by $re(\mathcal{X})^{(P^{(1)},\dots,P^{(G)})}_{i^{(1)},\dots,i^{(G)}}=\mathcal{X}_{d^{(1)}_{1},\dots,d^{(1)}_{N},d^{(2)}_{1},\dots,d^{(G)}_{N}}$, where $i^{(1)}=1+\sum_{g=1}^{N}(d^{(1)}_{g}-1)\prod_{n=1}^{g-1}D^{(1)}_{n}$.

Special cases of reordering are vectorization and matricization. Vectorization happens when $G=1$. $vec(\mathcal{X})=re(\mathcal{X})^{(\{-1\})}\in\mathbf{R}^{D}$, where $D=\prod_{n=1}^{N}D_{n}$ and $\{-1\}$ refers to the remaining unset modes. Mode-n matricization happens when $G=2$ and $P^{(1)}$ is a singleton. $mat(\mathcal{X})^{(n)}=re(\mathcal{X})^{(\{n\},\{-1\})}\in\mathbf{R}^{D_{n}\times D^{(\backslash n)}}$, where $D^{(\backslash n)}=\prod_{m=1(m\neq n)}^{N}D_{m}$.

We use graphical lasso as a part of our model. Given the mode-(N+1) matricization of the $(N+1)^{th}$-order TTS, $mat(\mathcal{X})^{(N+1)}\in\mathbb{R}^{T\times D}$, the graphical lasso (Friedman et al., 2008) estimates the sparse Gaussian inverse covariance matrix (i.e., network) $\theta\in\mathbb{R}^{D\times D}$, also known as the precision matrix, with which we can interpret pairwise conditional independencies among $D$ variables, e.g., if $\theta_{i,j}=0$ then variables $i$ and $j$ are conditionally independent given the values of all the other variables. The optimization problem is given as follows:

where $\theta$ must be a symmetric positive definite ($S^{p}_{++}$). $ll(x,\theta)$ is the log-likelihood and $\mu\in\mathbf{R}^{D}$ is the empirical mean of $mat(\mathcal{X})^{(N+1)}$. $\lambda\geq 0$ is a hyperparameter for determining the sparsity level of the network, and $\|\cdot\|_{od,1}$ indicates the off-diagonal $\ell_{1}$-norm. Since Eq. (1) is a convex optimization problem, its solution is guaranteed to converge to the global optimum with the alternating direction method of multipliers (ADMM) (Boyd et al., 2011) and can speed up the solution time.

A real-world complex $\mathcal{X}$ cannot be expressed by a single static network because it contains multiple sequence patterns, each of which has a distinct relationship/network. Moreover, we rarely know the optimal number of clusters and cluster assignments in advance. To address this issue, we want to provide an appropriate cost function and achieve subsequence clustering by minimizing the cost function. We now formulate the network-based TTS clustering problem. It assumes that $T$ time steps of $\mathcal{X}$ can be divided into $m$ time segments based on $K$ networks (i.e., clusters). Let $cp$ denote a starting point set of segments, i.e., $cp=\{cp_{1},cp_{2},\dots,cp_{m}\}$, the $i$-th segment of $\mathcal{X}$ is denoted as $\mathcal{X}_{cp_{i}:cp_{i+1}}$ where $cp_{m+1}=T+1$. We group each of the $T$ points into one of the $K$ clusters denoted by a cluster assignment set $\mathcal{F}=\{f_{1},f_{2},\dots,f_{K}\}$, where $f_{k}\subset\{1,2,\dots,T\}$, and we refer to all subsequences in the cluster $k$ as $\mathcal{X}[f_{k}]\subset\mathcal{X}$. Then, letting $\Theta$ be a model parameter set, i.e., $\Theta=\{\theta_{1},\theta_{2},\dots,\theta_{K}\}$, each $\theta_{k}\in\mathbb{R}^{D\times D}$ is a sparse Gaussian inverse covariance matrix that summarizes the relationships of variables in $\mathcal{X}[f_{k}]$. Therefore, the entire cluster parameter set is given by $\mathcal{M}=\{\mathcal{M}_{1},\mathcal{M}_{2},\dots,\mathcal{M}_{K}\}$, consisting of $\mathcal{M}_{k}=\{\theta_{k},f_{k}\}$. Overall, the problem that we want to solve is written as follows.

Assume $K,\mathcal{F}$ are given, here, we address how to define and infer the model $\theta_{k}$. The original graphical lasso allows $\theta_{k}$ to connect any pairs of variables in a tensor; however, it is too high-dimensional to reveal relationships separately in terms of the non-temporal modes. To avoid the over-representation, we aim to capture the multi-aspect relationships by separating $\theta_{k}$ into multimode to which we add a desired constraint for interpretability.

We assume that $\theta$ is derived from $N$ networks, $\{A^{(1)},\dots,A^{(N)}\}$, where $A^{(n)}\in\mathbf{R}^{D_{n}\times D_{n}}$ is the $n$-th network. For example, an element $a^{(n)}_{i,j}\in A^{(n)}$ refers to the relationship between the $i$-th and $j$-th variables of mode-n, In each network, the goal is to capture the dependencies between $D_{n}$ variables. We also assume that there are no relationships except among variables that differ only at mode-n. Thus, $\theta=\theta^{(N)}$ becomes an $N^{th}$ hierarchical matrix of shape $D\times D$. $\theta^{(n)}$ can be written as follows:

where $\theta^{(1)}=A^{(1)}$ and $C^{(n)}_{i,j}\in\mathbb{R}^{\prod_{m=1}^{n-1}D_{m}\times\prod_{m=1}^{n-1}D_{m}}$ is a diagonal matrix whose diagonal element is $a^{(n)}_{i,j}\in A^{(n)}$, i.e., $C^{(n)}_{i,j}=a^{(n)}_{i,j}\cdot\delta_{i.j}$ allows edges that differ only at mode-n, where $\delta_{i.j}$ is the Kronecker delta.

We extend graphical lasso to obtain $\theta$ by inferring a sparse $A^{(n)}$ from a TTS. The optimization problem is written as follows:

where $\mu_{d}\in\mathbb{R}^{D_{n}}$ is the empirical mean of the variable $re(\mathcal{X})_{:,d,:}\in\mathbb{R}^{T\times D_{n}}$. Eq. (8) is a convex optimization problem solved by ADMM. We divide the log-likelihood by $D^{(\backslash n)}$ to scale the sample size.

To determine the cluster assignment set $\mathcal{F}$ and the number of clusters $K$, we use the MDL principle (Grünwald, 2007), which follows the assumption that the more we compress the data, the more we generalize its underlying structures. The goodness of the model $\mathcal{M}$ can be described with the following total description cost:

We describe the four terms that appear in Eq. (10).

Our next goal is to find the best cluster parameter set $\mathcal{M}$ that minimizes the total description cost Eq. (10).

The first goal is to divide a given $\mathcal{X}$ into $m$ segments (i.e., patterns), but we assume that no information is known about them in advance. Therefore, to prevent a pattern explosion when searching for their optimal cut points, we introduce CutPointDetector based on the divide-and-conquer method (Keogh et al., 2001).

Specifically, it recursively merges a small segment set of $\mathcal{X}$ while reducing its total description cost, because neighboring subsequences typically exhibit the same pattern. We define $\mathbf{w}$ as a set of user-defined initial segment sizes, i.e., $\mathbf{w}=\{w_{i}\}_{i=1}^{m}$, such as the number of days in each month or any small constant. An example illustration is shown in Fig. 2. Let $\theta_{i:i+1}$ be a model of $\mathcal{X}\{cp_{i}:cp_{i+1}\}$ at the $i^{th}$ segment. Given the three subsequent segments illustrated in Fig. 2 (a), we evaluate whether to merge the middle segment with either of the side segments (Fig. 2 (b)(c)). The total description cost for Fig. 2 (a) is given by $Cost_{T}(\mathcal{X};\{\theta_{i:i+1},\theta_{i+1:i+2},\theta_{i+2:i+3}\})$, where we omit the cluster assignment (e.g., $\{j\}_{j=cp_{i}}^{cp_{i+1}-1}\}$) from the cost for clarity. If the cost for the original three segments is reduced by merging, it eliminates the unnecessary cut point and employs a new model $\theta$ for the merged segment. By repeating this procedure for each segment, $m$ decreases monotonically until convergence. See Appendix B.1 for the detailed procedure.

DMM searches for the best number of clusters by increasing $K=1,2,\dots,m$, while the total description cost $Cost_{T}(\mathcal{X};\mathcal{M})$ is decreasing. To compute the cost, however, we must solve two problems, namely obtain the cluster assignment set $\mathcal{F}$ and the model parameter set $\Theta$, either of which affects the optimization of the other. Therefore, we design ClusterDetector with the expectation and maximization (EM) algorithm. In the E-step, it determines $\mathcal{F}$ to minimize the data coding cost, $Cost_{C}(\mathcal{X}|\mathcal{M})$, which is achieved by solving:

for the $i$-th segment, and then inserts time points from $cp_{i}$ to $cp_{i+1}$ (i.e., $\{j\}_{j=cp_{i}}^{cp_{i+1}-1}$) to the best $k$-th cluster $f_{k}\in\mathcal{F}$. In the M-step, for $1\leq k\leq K$ the algorithm infers $A^{(n)}_{k}(1\leq n\leq N)$ according to Eq. (8) to obtain $\theta_{k}\in\Theta$ for a given $\mathcal{X}[f_{k}]$. Note that ClusterDetector starts by randomly initializing $\Theta$.