【论文翻译】Toeplitz Inverse Covariance-Based Clustering of Multivariate Time Series Data

摘要

多元时间序列的子序列聚类是发现时态数据中重复模式的有用工具。 一旦这些模式被发现，看似复杂的数据集就可以解释为仅有少量状态或集群的时间序列。 例如，来自健身跟踪应用程序的原始传感器数据可以表示为选定的几个动作（即步行、坐着、跑步)的时间线。 然而，发现这些模式是具有挑战性的，因为它需要同时分割和聚类时间序列。 此外，解释由此产生的簇是困难的，特别是当数据是高维的。 在这里，我们提出了一种新的基于模型的聚类方法，我们称之为TICC。在TICC方法中，每个集群都是由相关网络或马尔可夫随机场(MRF)定义的，它描述了该集群的典型子序列中不同观测值之间的相互依赖关系。 基于这种图形表示，TICC同时分割和聚类时间序列数据。 利用期望最大化(EM)算法的变形——交替最小化来解决TICC问题。我们分别通过动态规划和交替方向乘子法(ADMM)，导出了封闭形式的解以可伸缩的方式有效地求解这两个子问题。 我们通过在一系列合成实验中将TICC与几个最先进的基线进行比较来验证我们的方法，然后在汽车传感器数据集上演示 TICC 如何用于学习现实世界场景中的可解释集群。

1.介绍

许多应用，从汽车[32]到金融市场[35]和可穿戴传感器[34]，产生大量的时间序列数据。 在大多数情况下，这些数据是多元的，其中每个时间戳观测由来自多个实体或传感器的读数组成。 这些长时间序列通常可以分解成一个状态序列，每个状态由一个简单的“模式”定义，其中状态可以重复发生多次。例如，来自健身跟踪设备的原始传感器数据可以解释为动作[38]（比如走10分钟，跑30分钟，坐一个小时，再跑45分钟）的时间序列。同样，使用汽车传感器数据，单个驾驶会话可以表示为几个关键状态的顺序时间线：转弯、加速、减速、直行、红灯停车等。 这种表示可以用来发现重复的模式，了解趋势，检测异常，更广泛地说，更好地解释大型和高维数据集。

图1. 我们的TICC方法将一个时间序列分割成一系列状态，或“簇”(即a、b或c)。 每个集群的特征是一个相关网络，或MRF，定义在一个大小为W的短窗口上。这个MRF控制属于该集群的段内任何窗口的（时间不变）偏相关结构。 在这里，TICC学习集群MRFS和时间序列分割。

为了实现这种表示，需要同时分割和聚类时间序列。 这个问题比标准时间序列分割[17,20]更困难，因为多个片段可以属于同一个集群。 然而，它也比子序列聚类困难[3,43]因为每个数据点不能独立聚类（因为相邻点最好属于同一个簇）。 此外，即使一个人能够同时分割和聚类数据，问题仍然是如何解释不同的集群。 这些集群很少事先知道，因此最好通过数据来学习。 然而，如果没有先验知识，就很难理解每个集群所指的是什么。传统的聚类方法并不特别适合于在数据中发现可解释的结构。 这是因为它们通常依赖于基于距离的度量，例如动态时间扭曲[4]。这些方法侧重于匹配原始值，而不是在数据中寻找更细微的结构相似性，例如，汽车中的不同传感器在不同时间之间是如何相互关联的。在本文中，我们提出了一种新的多元时间序列聚类方法，我们称之为托普利兹逆协方差聚类(TICC)。 在我们的方法中，我们将每个集群定义为一个依赖网络，在一个短时间内子序列中显示不同传感器之间的关系（图1）。例如，在一个与汽车“转弯”对应的簇中，这个被称为马尔可夫随机场(MRF)的网络可以显示通用时刻t的制动脚板如何影响 t + 1 时刻的方向盘角度。在这里，不同集群的MRF，比如“slow down”，将在这两个传感器之间有一个非常不同的依赖结构。在这些 MRFs中，一条边表示两个变量之间的偏相关[26,41,46]。重要的是要注意 MRFs表示的关系比简单的相关性强得多;偏相关是用来控制其他混杂变量的影响，因此MRF中边的存在意味着数据中两个变量之间存在直接的依赖关系。因此，MRF 提供了可解释的洞察力，以精确地描述每个集群的关键因素和关系。
在我们的TICC方法中，我们通过估计一个稀疏高斯逆协方差矩阵来学习每个簇的MRF[14,48]。逆协方差Θ,如果Θi, j = 0,那么根据定义,在Θ中元素 i 和 j 是条件独立的(给定所有其他变量的值)。因此,Θ定义了MRF依赖网络的邻接矩阵[1,45]。这个网络有多个层，在一个层内和不同层之间都有边。在这里，层的数量对应于我们定义的MRF的短子序列的窗口大小。例如，图1中对应于簇A和簇B的MRFs都有三层。该多层网络代表了属于该簇的某个段内观测的任意窗口的非时变相关结构。我们通过解决一个约束逆协方差估计问题来学习每个聚类的这种结构，我们称之为Toeplitz graphical lasso。我们所施加的约束确保所得到的多层网络具有块Toeplitz结构[16] (即l 层和l + 1层之间的任何边也存在于l + 1层和l + 2层之间)。这个Toeplitz约束保证了我们的聚类定义是非时变的，所以聚类分配不依赖于子序列的确切起始位置。相反，我们只根据时间序列当前的结构状态对这些短子序列进行聚类。
为了解决TICC问题，我们使用期望最大化(EM)方法，基于交替最小化，迭代地聚类数据，然后更新聚类参数。尽管TICC涉及解决一个高度非凸的极大似然问题，但我们的方法在实践中能够非常有效地找到(局部)最优解。在将数据分配到聚类时，我们还有一个时间一致性的目标，即鼓励时间序列中相邻数据属于同一聚类。然而，这产生了一个组合优化问题。我们使用动态规划开发了一个可扩展的解决方案，它允许我们有效地学习最优分配，只需O(KT)时间来将T点分配到K个簇。然后，为了求解聚类参数，我们开发了一个求解Toeplitz graphical lasso问题的算法。由于从数据中学习这个图形结构是一个计算开销很大的半定规划问题[6,22]，因此我们开发了一种专门的基于交替方向乘子[5]方法的消息传递算法。在Toeplitz图形拉索算法中，我们推导了每个ADMM子问题的封闭形式更新，以显著加快求解时间。
然后，我们实现了我们的TICC方法，并将其应用于真实和合成数据集。我们从几个已知的地面实况聚类的综合例子开始评估性能。我们将 TICC 与几种最先进的时间序列聚类方法进行了比较，在聚类分配精度方面至少比它们都高出41%。我们还对每种方法精确的集群恢复所需的数据量进行了量化，我们发现，要实现类似的性能，TICC需要的观察量比第二优的方法少三倍。此外，我们发现我们的方法能够准确重建底层的MRF依赖网络，在我们的实验中网络恢复得分F1是0.79-0.90。然后我们分析一个汽车传感器数据集，以了解如何使用TICC从现实数据中学习可解释的见解。应用我们的方法，我们发现汽车数据集有五个真实的聚类，每个聚类对应汽车经常处于的“状态”。然后，我们通过检查驾驶会话的纬度/经度位置以及结果的聚类分配来验证我们的结果，以展示 TICC如何成为从多元时间序列进行无监督学习的有用工具。

相关工作

这项工作与最近在时间序列聚类和凸优化方面的进展有关。时序数据的子序列聚类是一个比较成熟的领域。方法包括动态时间扭曲的几种变化[3,23,25,39]、符号表示[29,30]和基于规则的motif发现[11,28]。对于时间序列数据的同时聚类和分割也有研究，称为时间点聚类[15,49]。然而，这些方法通常依赖于基于距离的度量，在某些情况下，已经显示出不可靠的结果[24]。相反，我们的TICC方法是一种基于模型的聚类方法，类似于基于ARMA [47]， 高斯混合[13]或隐马尔科夫模型[43]的聚类方法。据我们所知，我们的方法是第一个基于每个子序列的图形依赖性结构来进行时间序列聚类的方法。这为我们的集群提供了可解释性，防止了过度拟合，并且，正如我们在第6节和第7节中所展示的，允许我们发现其他方法无法发现的模式类型。为此，我们提出了一个结构化协方差反估计问题，我们称之为Toeplitz图形套索。这个问题是著名的图形lasso问题[14]的变体，在这个问题中，我们在解上强制执行一个块Toeplitz结构。虽然存在许多算法来解决标准图形lasso[1,21,22]，但我们不知道有任何方法专门适用于Toeplitz块。本文提出了一种基于ADMM算法的优化问题求解方法，该方法可以将全局优化问题分解为对 ADMM 友好的子问题，在子问题中我们可以推导出闭型近邻运算子[5]来快速求解优化问题。

2.构建问题

考虑一个T序列观测的时间序列，

式中，x_i为第i个多元观测。我们的目标是把这T个观测值聚到K个聚类中。然而，我们不是孤立地聚类每一个观测，而是在时间序列的前几个点的上下文中处理每一个点。因此，不仅仅是看x_t，我们反而聚集了一组尺寸为w的时间序列，这包括观察x_t−w+1，…, x_t，我们把它连接成一个nw维的向量，我们称之为X_t。我们将这个从X₁到X_T的新序列称为X。注意，在每个点x_T与其子序列X_T之间存在一个双射，或者说双向一对一映射.(x_orig的第一个w观测结果映射到一个更短的子序列，因为时间序列直到x₁才开始）。这些子序列是一个有用的工具，可以为每个T观测提供适当的上下文。例如,在汽车行驶过程中、一个观察可能显示汽车的当前状态(例如,开车直行15英里/小时),但一个短窗口,甚至只有几分之一秒,允许更完整的数据理解(例如,汽车是否加速或者减速)。因此，我们的方法不是直接对观测结果进行聚类，而是对这些子序列X1，…, XT。我们这样做是为了鼓励相邻的子序列属于同一个集群，这个目标我们称之为时间一致性。因此，我们的方法可以看作是[49]时间点聚类的一种形式，同时对时间序列进行分段和聚类。
Toeplitz Inverse Covariance-Based Clustering (TICC). 我们用一个逆高斯协方差Θ_i∈R ^nw×nw定义每个集群。回想一下,逆协方差之间的显示条件独立性的结构变量[26],所以Θ_i定义了一个马尔可夫随机场来编码聚类i的结构表示。除了提供可解释的结果，稀疏图形表示是防止[27]过拟合的有用方法。因此,我们的目标是解决这些K逆协方差Θ={Θ1,…,ΘK}, 每个集群一个,和由此产生的结果集P = {P1,…, PK}，其中Pi⊂{1,2，…T}。在这里，每个T点被精确地分配给一个簇。我们的总体优化问题是
这就是TICC问题。$\tau$ 为Toeplitz对称块 nw×nw矩阵的集合,${\Vert \lambda \circ {\Theta }_i\Vert }_1$是对Hadamard乘积的一范式惩罚，用来激活稀疏逆协方差（$\lambda$是正则参数）。此外，$ll\left(X_t,{\Theta }_i\right)$是X_t来自簇 i 的对数似然值。

其中，${\mu }_i$是群集 i 的经验均值。$\beta$是一个强制时间一致性的参数， I { X t − 1 ∉ P i } \mathbb{I}\left \{ X_{t-1}\notin P_{i} \right \} I{Xt−1​∈/​Pi​}是一个指示函数，用来检查相邻的点是否被分配到同一个簇中。

Toeplitz Matrices 注意,我们限制逆协方差$\Theta$是Toeplitz块。因此，每个nw×nw矩阵可以表示为

其中，$A^{(0)},A^{(1)},...,A^{(w-1)}\in {\mathbb{R}}^{n\times n}$。子块$A^{(0)}$表示时间内的偏相关性，因此$A_{i,j}^{(0)}$表示传感器i与j当前值的关系。在该簇对应的MRF中,$A^{(0)}$定义每一层内边的邻接矩阵。另一方面，非对角的子块指的是“跨时间”边。例如,$A_{i,j}^{(1)}$表示t时刻的传感器i与t + 1时刻的传感器j之间是如何相关的,$A^{(2)}$示时间t和时间t + 2之间的边结构。逆协方差的Toeplitz块结构意味着我们在这个长度w窗口上做了一个定常假设(我们通常希望这个窗口的大小比平均段长度小得多)。因此，例如在图1中，层1和层2之间的边也必须存在于层2和层3之间。我们使用这个假设是因为我们正在寻找一种独特的结构模式来识别每个集群。我们将每个集群视为特定的“状态”。当时间序列处于这种状态时，不管窗口的起始点是什么，它都会保留一个特定的(定常的)结构，这个结构会在整个片段中持续存在。通过在逆协方差上执行一个托普利兹结构,这次我们能够模拟不变性，将${\Theta }_i$添加到我们的估计中。

Regularization Parameters 我们的TICC优化问题有两个正则化参数:$\lambda$确定描述每个簇的MRFs中的稀疏度;$\beta$是鼓励相邻子序列被分配到同一个集群的平滑性惩罚。注意，即使$\lambda$是nwxnw维矩阵，我们通常将它的所有值设置为一个常量，从而将搜索空间减少到只有一个参数。在有适当稀疏性或时间一致性先验知识的应用中，可手动选择$\lambda$和$\beta$。通常，参数值的选择也可以采用更有原则的方法，如贝叶斯信息准则(BIC)[18]或交叉验证。
Window Size 回想一下，我们不是孤立地聚类每一个点xt，而是聚类一个短窗口，或者说是子序列，从时间t - w + 1到t，我们将其连接成一个nw维的向量，我们称之为Xt。Toeplitz约束假设每个集群都有一个非时变结构，但这个窗口大小仍然是一个相关参数。特别地，它允许我们学习交叉时间相关性(例如，在t时刻的传感器i影响在t +1时刻的传感器j)。窗口越大，这些跨时间边界的距离就越远。然而，我们不希望我们的窗口太大，因为它可能会在分割边界上正确分类点存在困难，在那里我们的非时变假设可能不成立。因此，我们一般保持w的值相对较小。但是，它的准确值通常应该根据应用程序、观察的粒度和平均预期段长度来选择。它也可以通过BIC或交叉验证来选择，尽管我们在第6节中发现，我们的TICC算法对于这个窗口大小参数的选择是相对稳健的。
Selecting the Number of Clusters 与许多聚类算法一样，TICC中聚类的数量K是一个重要的参数。有不同的方法可以做到这一点。如果有一些标记的地面实况数据，我们可以在测试集上使用交叉验证或标准化互信息[8]来评估性能。如果我们没有这样的数据，我们可以使用BIC或者剪影评分[40]来选择这个参数。然而，集群的确切数量通常取决于应用程序本身，特别是在我们还在寻找精确性之外的可解释性的情况下。

3.交替最小化

问题(1)是一个混合组合连续优化问题。有两组变量,集群结果集P和逆协方差Θ,两者都是耦合在一起的高度非凸问题。因此，没有可控制的方法来求解全局最优解。相反，我们使用期望最大化(EM)算法的一种变化，在分配点给集群和更新集群参数之间交替进行。虽然这种方法不一定达到全局最优，但类似类型的方法已被证明在相关[13]问题上表现良好。在这里，我们定义了组成我们方法的两个步骤的子题。然后,在第4节，我们推导了快速算法来解决子问题，并正式描述了我们的整体算法来解决TICC问题。

3.1给集群分配点

我们通过整理Θ的值来分配点到各个聚类，解决以下P组合优化问题 P = { P 1 , . . . , P K } P=\left \{ P_{1},...,P_{K} \right \} P={P1​,...,PK​},

该问题将每个T子序列分配到K个聚类中的一个，以联合最大化对数似然性和时间一致性，并通过正则化参数函数$\beta$调节两个目标之间的权衡。当$\beta$= 0时，表示子序列X1，…，XT可以独立分配，因为没有惩罚来鼓励相邻的子序列属于同一簇。这可以通过简单地将每个点分配到概率最大的集群来解决。随着$\beta$变大，相邻的子序列越来越有可能被分配到同一个集群中。当$\beta \to \infty$，切换代价变得非常大，以至于时间序列中的所有点都被归为一个簇。尽管问题3是组合问题，在4.1节我们可以看到可以使用动态规划来有效地寻找这个TICC子问题的全局最优解。

3.2 Toeplitz Graphical Lasso

给定结果集P,我们的下一个任务是通过解决问题(1)来更新集群参数${\Theta }_1,...,{\Theta }_K$,同时保持P不变。我们可以并行求出每个Θi 。为此,我们可以针对每个Θ重写问题(2)中的负对数似然。这种似然可以表示为：

其中，$\left|P_i\right|$是分配给簇i的点数，Si是这些点的经验协方差,和C是一个不取决于Θi的常数。因此，我们的EM算法的m步为：

问题(4)是一个凸优化问题，我们称之为托普利兹图形套索。这是著名的图形lasso问题[14]的一个变形，在这个问题中，我们在逆协方差上添加了一个Toeplitz块约束。原图形套索定义了两个目标之间的权衡,受参数约束:减少负对数似然,确保Θi是稀疏的。当Si是可逆的,Θi会靠$S_i^{-1}$.我们需要添加额外的约束确保Θi是Toeplitz块。λ是一个nw×nw矩阵,因此它可以用来规范Θi的每一子块不同。可以通过简单的尺度化$\lambda$来将$\frac{1}{|P_i|}$纳入正则化。因此，为了简化符号，我们通常在写问题4的时候略去这一项(相应地缩放$\lambda$)。

4.TICC算法

在这里，我们将描述我们的算法，将X1，…，将XT聚类为K个簇。我们的方法，在章节4.3中详细描述，依赖于两个关键子例程:在AssignPointsToClusters中，我们使用一个动态规划算法来将每个Xt分配到一个集群中;在UpdateClusterParameters中，我们通过使用基于交替方向乘数法(ADMM)的算法来解决Toeplitz图形套索问题来更新集群参数。注意，这类似于期望最大化(EM)，两个子例程分别对应于E步和M步。

4.1 聚类任务

给定K个聚类的模型参数(即逆协方差)，求解问题(3)，分配T个子序列X1，…，XT到K个聚类中，同时最小化集群分配在时间序列中变化的次数。鉴于T点的K个可能的聚类分配，这个组合优化问题有$K^T$个可能的解。然而，我们只能在O(KT)操作中求出全局最优解。我们通过算法（1）中描述的动态规划方法来做到这一点，该方法相当于为长度t序列寻找最小代价维特比路径[44]，如图2所示。

问题(3)等价于寻找从时间戳1到时间戳T的最小代价路径，其中节点代价是该点被分配到一个给定集群的负对数似然概率，而每当集群分配切换时，边代价是$\beta$.

4.2 求解Toeplitz图形套索

一旦我们完成聚类分配，EM算法的m步就是更新给定分配给每个聚类的点的逆协方差。这里，我们解决问题(4)中定义的Toeplitz图形套索。对于较小的协方差，这个半定规划问题可以用标准内点法求解[6,36]。然而，为了解决整体的TICC问题，我们需要在算法的每次迭代中为每个集群解决一个单独的Toeplitz图形套索。在TICC收敛之前，我们可能需要解决问题(4)数百次，因此有必要快速解决该问题。我们通过乘数交替方向法(ADMM)来实现这一点，这是一种分布式凸优化方法，已被证明在大规模优化任务中表现良好[5,37]。利用ADMM将问题分解为两个子问题，并使用消息传递算法迭代收敛到全局最优解。当ADMM子问题可以找到封闭形式的解决方案时，ADMM具有可扩展性，我们可以推导出Toeplitz图形套索。
为了使问题(4)具有ADMM友好的形式，我们引入一个一致变量Z，并将问题(4)重写为其等价问题，

增广拉格朗日[19]可以表示为

其中，$\rho >0$是ADMM的惩罚参数，U∈R ^nw×nw为尺度对偶变量[5，§3.1.1]，Z∈$\tau$。
ADMM由以下三个步骤组成，重复直到收敛:

其中k为迭代次数。这里，我们通过Θ和Z交替优化问题(5),在每次迭代之后我们更新了双重变量U。由于Toeplitz图解套索问题是凸的，ADMM保证收敛到全局最优。我们使用了基于原始残值和对偶残值接近零的停止准则;参见[5]。
Θ-Update Θ- update可以写成

该优化问题存在一个已知的解析解[10]，

其中，$QD{Q}^T$是$\frac{Z^k-U^k}{\rho }-S$的特征分解。
Z-Update Z-update可以写成：

该近邻算子可以对每个子块$A^{\left(0\right)},...,A^{\left(w-1\right)}$并行求解。此外，在每个子块中,子块中的第(i, j)个元素也可以并行求解。在$A^{\left(0\right)}$中，有$\frac{n\left(n+1\right)}{2}$个独立问题（它是对称的），其他w-1个块中，这个数字是$n^2$.因此,问题(7)可以分解为$\left(w-1\right)n^2+\frac{n\left(n+1\right)}{2}$个独立的问题。每一个问题都有相同的形式，可以用相同的方式解决。

我们用R表示每个元素出现的次数（对于子块A^m，这等于2(w−m),除了A⁽⁰⁾的对角线,发生w次）。我们对R出现次数进行排序，我们让$B_{ij,l}^{\left(m\right)}$表示A^m的第(i, j)个元素第l次出现在Z中的索引，l = 1, . . . , R。因此$B_{ij,l}^{\left(m\right)}$返回Z的nwxnw矩阵中的索引（x,y).用这个符号，我们可以用相同的方法解出每一个$\left(w-1\right)n^2+\frac{n\left(n+1\right)}{2}$子问题的Zupdate近邻算子。求出Z对应于$B_{ij}^{\left(m\right)}$的元素，我们让这些元素都等于：

让$Q={\sum }_{l=1}^R\left|{\lambda }_{B_{ij,l}^{(m)}Z}\right|$,$S_l={\left({\Theta }^{k-1}+U^k\right)}_{B_{ij,l}^{(m)}}$,问题8可以简化为

这只是一个软阈值近端算子[37]，它有以下封闭形式的解，

我们在$Z^{k+1}$中填入R个元素(对应于$B_{ij}^{\left(m\right)}$)。我们对所有$\left(w-1\right)n^2+\frac{n\left(n+1\right)}{2}$个子问题做同样的处理，每一个都可以并行求解，我们就得到了整个Z更新的最终结果。

4.3 TICC聚类算法

为了解决TICC问题，我们将4.1节中的动态规划算法和4.2节中的ADMM方法结合成一个迭代EM算法。我们从随机初始化集群开始。从那里，我们交替使用E步和m步，直到集群分配平稳(也就是说，问题已经收敛)。算法(2)概述了TICC的总体方法。

5 完成

我们已经建立了一个自定义的Python求解器来运行TICC算法。我们的求解器以原始多元时间序列和问题参数作为输入。然后返回时间序列中每个点的聚类分配，以及每个聚类的结构MRF表示。

6 实验

我们在几个综合算例上测试了TICC方法。我们这样做是因为有已知的“地面实况”聚类来评估我们的方法的准确性。
生成数据集 我们在R⁵中随机生成合成的多元数据。每个K集群平均值是0,聚类结果完全基于数据的结构。对于每个聚类，我们生成一个随机的地面实况Toeplitz逆协方差[33]:

然后通过构建聚类段的时间序列生成整体时间序列(例如,在三个聚类段中，每一个都有序列“1、2、1”的200个样本,来自两个逆协方差Θ1和Θ2)。然后根据之前w - 1个样本的值，一次抽取一个样本。请注意，当我们刚刚切换到一个新集群时，我们正在部分基于前一个集群生成的数据绘制一个新示例。
我们在四种不同的时间序列上进行实验:“1、2、1”,“1、2、3、2、1”,“1、2、3、4、1、2、3、4”,“1、2、2、1、3、3、3、1”。每个例子中的每个片段在R⁵中都有100K个观察值，其中K为该实验的簇数(分别为2,3,4,3)。这些例子被选择来表达不同类型的时间序列在不同的时间长度。
性能标准 我们通过聚类时间序列中的每个点并与地面实况聚类进行比较来评估性能。由于TICC和基线方法都使用非常相似的方法来选择适当的集群数量，所以我们将K固定为对TICC和所有基线来说的“真实”集群数量。这就产生了一个简单的多分类问题，它允许我们通过度量macro-f1得分来评估聚类的准确性。对于每个聚类，F1得分是我们估计的精度和召回率的调和平均值。然后，macro-F1得分是所有集群的F1得分的平均值。我们使用这个分数来比较我们的TICC方法与几个著名的时间序列聚类基线。
Baseline Methods. 们使用多模型和基于距离的聚类方法作为我们的基线。我们使用的方法是:
TICC，$\beta =0$——这是我们的没有时间一致性约束的TICC方法。在这里，每个子序列被分配给一个独立于其在时间序列中的位置的集群。
GMM——采用高斯混合模型[2]进行聚类。
EEV——正则化的GMM，对高斯协方差矩阵进行形状和体积约束[13]。
DTW，GAK——基于动态时间翘曲(DTW)的聚类，使用全局对齐内核[9,42]。
DTW，Euclidean——DTW使用欧几里德距离度量[42]。
Neural Gas——基于自组织映射的人工神经网络聚类方法[12,31]。
K-means——使用欧氏距离的标准的K-means聚类算法。
聚类精度 我们在表1中测量了四种不同时间序列的宏观f1得分。在这里，所有八种方法都使用完全相同的合成数据，以分离每种方法对性能的影响。如图所示，TICC显著优于基线。我们的方法获得了介于0.90和0.98之间的宏观f1得分，四个例子的平均得分为0.95。这比第二好的方法GMM高出41%（不算TICC, β = 0），平均宏观f1得分只有0.67分。我们还使用micro-F1评分来运行我们的实验，它使用一个加权平均值来对拥有更多样本的集群进行更重的衡量，我们得到了非常相似的结果(在macro-F1评分的1-2%范围内)。请注意，我们示例中的K个集群始终是零均值，它们只根据数据的结构有所不同。在这些实验中，基于距离的方法比基于模型的方法得分低。

Effect of the Total Number of Samples. 接下来我们关注每种方法需要多少样本才能准确地对时间序列进行聚类。我们取表1中“1,2,3,4,1,2,3,4”的例子，改变样本数。我们在图3中绘制了8种方法的宏f1得分与每个片段的样本数量的对比图。可以看出，当样本数为100时，没有一种方法能够准确地聚类数据。但是，当我们观察到更多的样本时，无论是TICC还是TICC，$\beta =0$都能快速提高性能。当有200个样本时，TICC的宏观f1得分已经在0.9以上。即使是在有限的数据量，我们的TICC方法仍然能够准确地聚类数据。此外，我们注意到，由$\beta$定义的时间一致性约束在这个区域只有很小的影响，因为TICC和TICC，$\beta =0$达到类似的结果。因此，由于我们在TICC方法中施加的稀疏块Toeplitz约束，最可能得到准确的结果。然而，随着样本数量的增加，这两个图开始出现分歧，随着TICC趋于1.0，而TICC, $\beta =0$在0.9左右徘徊。这意味着，一旦我们有足够的样本，最终的性能改进是由于时间一致性惩罚，鼓励邻近的样本被分配到同一簇。
Network Recovery Accuracy. 我们的TICC方法还有一个额外的好处，那就是它学习的集群是可解释的。TICC将每个集群建模为多层马尔可夫随机场,边对应逆协方差矩阵中的非零元素Θi.我们可以将我们估计的网络与“真实的”MRF网络进行比较，并测量我们估计的所有集群的平均宏f1得分。我们查看与表2相同的四个示例，并将结果绘制在表2中。我们恢复了F1值在0.79到0.90之间的网络底层边缘结构。这说明TICC既能准确地对数据进行聚类，又能恢复底层聚类的网络结构。请注意，我们的方法是第一种能够显式地重建这个网络的方法，而其他基线方法无法做到这一点。
Window Size Robustness. 接下来，我们将研究窗口大小w的选择如何影响结果。我们运行相同的“1,2,3,4,1,2,3,4”示例，除了现在我们改变窗口大小w.(回忆一下“true”窗口大小为5。)根据经验，我们发现，任何窗口大小在4到15之间的宏观f1聚类精度得分在0.95到0.98之间。类似地，对于窗口大小在5到14之间的情况，我们的网络恢复宏f1得分保持在0.87到0.89之间。只有当窗口大小降到4以下或15以上时，结果才开始变得更糟。因此，我们的TICC方法对w的选择似乎是相对稳健的。

Scalability of TICC ICC算法的一次迭代包括运行动态规划算法，然后求解每个簇的Toeplitz图形拉索问题。重复这些步骤，直到收敛。迭代的总次数取决于数据，但通常不超过几十次迭代。由于T通常比K和n都大得多，我们可以预期在分配阶段会出现最大的瓶颈，其中T可能有数百万。为了评估算法的可伸缩性，我们改变T并计算算法在一次迭代中的运行时间。我们观察R⁵⁰中的样本，估计窗口大小为3的5个集群，并在几个数量级上变化。我们在图4中用对数-对数尺度绘制结果。注意，我们的ADMM求解器在4秒内推断每个150×150逆协方差(因为nw = 50×3 = 150)，但这个运行时间与T无关，所以ADMM对图中的常数偏移量有贡献。如图所示，当T值较大时，我们的算法与点的数量呈线性缩放。我们的TICC求解器可以聚类1000万个R⁵⁰中点，每个迭代运行时间大约为25分钟。

7 案例学习

在这里，我们将我们的TICC方法应用到一个真实的例子中，以演示如何使用这种方法在无监督的情况下从时间序列数据中找到有意义的见解。
我们分析由一家大型汽车公司提供的数据集，其中包含来自真实驾驶会话的传感器数据。这一环节持续1小时，并发生在一个欧洲大城市郊区的真实道路上。我们每0.1秒观察7个传感器:

因此，在这一小时的会议中，我们获得了36000次x7维的时间序列。我们窗口大小为1秒(或10个样本)。我们使用BIC来选择簇的数量，最终确定K=5。

我们通过对TICC输出的5个集群进行分析，来理解和解释它们各自所指的“驾驶状态”。每个集群都有一个定义其结构的多层MRF网络，为了分析结果，我们使用网络分析来确定簇网络中每个店的相对的“重要性”。我们将每个节点的介数中心性评分[7]绘制在表3中。我们可以看到，5个集群中的每一个都有一个独特的“签名”，不同的传感器在每个集群中具有不同的介数中心性分数。例如，y -加速度传感器在5个集群中只有2个非零分:#2和#5。因此，我们期望这两个集群指的是汽车转弯的状态，而其他三个集群指的是汽车直线行驶的时间间隔。类似地，集群#1是惟一跟油门踏板无关，也是刹车踏板得分最大的集群。因此，我们希望无论什么时候汽车减速时，这个状态就是集群分配。我们也看到,集群3和4有相同的非零传感器、速度和油门踏板,所以我们可能期望他们指州连续驾驶汽车时不减速,最“重要”的传感器在两个集群是集群4中的速度。因此，我们可以使用这些中间值来以一种有意义的方式解释这些集群。例如，从表3中，一个合理的假设可能是，集群指的是1)减速，2)转弯，3)加速，4)在直道上巡航，5)在弯曲路段行驶。

驾驶过程中的两个真实世界的转弯。指示针的颜色代表从我们的TICC算法中分配的集群(绿色=走直线，白色=减速，红色=转弯，蓝色=加速)。由于我们基于结构而不是距离进行聚类，因此在TICC下左转和右转看起来非常相似。

Plotting the Resulting Clusters 为了验证我们的假设，我们可以绘制驱动器的纬度/经度位置，以及生成的集群分配。通过分析这些数据，我们实证地发现，五个集群中的每一个都有一个清晰的现实世界的解释，这与我们基于表3中中介分数的估计非常接近。此外，我们注意到在这一小时中出现了许多一致和重复的模式。例如，每当司机接近转弯时，他或她都会遵循同样的顺序:直行、减速、转弯、加速，然后再次直行。我们在数据集中绘制了两个典型的回合，根据它们的集群分配给时间戳着色，如图5所示。重要的是要注意，在这里，同样的模式出现在左转弯和右转弯。虽然基于距离的方法处理这两种情况会非常不同(因为几个传感器的值完全相反)，但TICC却根据结构相似性对时间序列进行聚类。这样，左右分配的TICC就变成了同一个底层集群。

8 结论和展望

文定义了一种多元时间序列子序列的聚类方法。我们的方法，Toeplitz逆协变聚类(TICC)，是一种新型的基于模型的聚类，能够在数据中找到准确的和可解释的结构。我们的TICC算法同时对数据进行分段和聚类，将高维时间序列分解成清晰的顺序时间线。我们根据子序列的结构相关性对其进行聚类，通过多层MRF定义每个集群，使我们的结果具有高度的可解释性。为了发现这些簇，TICC通过动态规划和ADMM来更新簇的MRFs，以时间一致的方式分配点给簇。
TICC在合成数据和真实世界数据上的有前景的结果为未来的研究带来了许多潜在的方向。例如，我们的方法可以扩展到任意异质指数族MRF参数化的依赖网络。这将允许更广泛的数据集(如布尔值或分类读数)被纳入到现有的TICC框架中，为这项工作提供新的潜在应用。