Toeplitz Inverse Covariance-based Clustering ~時系列クラスタリング~

　今回の記事では，time series clustering（時系列クラスタリング）に関する概要を述べ，その後に KDD 2017 paper である TICC (Toeplitz Inverse Covariance-based Clustering) という手法について解説する．
　また，説明の際に Markov Random Filed（マルコフ確率場, MRF），Expectation Maximumzation algorithm（EM algorithm）, Dynamic Programming（動的計画法, DP）, ADMM を用いる．MRF の説明は軽く行い，EM, DP, ADMM に関しては理解していなくても読み進められると思う．

# Time Series Clustering
　Time series clustering には，大きく３つのタスク（種類）がある．

* Whole time series clustering
    * 個々の時系列集合に対するクラスタリングのこと
    * (ex) 患者100人の脳波データを似ているデータごとにクラスタリング
* Subsequence time series clustering
    * 一本の時系列データのセグメント（時系列データをある間隔ごとに区切ること）に対するクラスタリングのこと
    * (ex) 100秒ごとに区切られた一人の患者の脳波データをクラスタリング
* Time point clustering
    * 一本の時系列データの各時点ごとのクラスタリング．ただし，どのクラスにも属さない（ノイズと判断される）時点もある
    * (ex) 一人の患者の脳波データの各時点の値に対するクラスタリング

　今回，説明していく TICC は time point clustering，subsequence time series clustering を同時に行える手法である[^1]．
[^1]: 正確には，一本の時系列データの segment と subsequence time series clustering を同時に行う手法である．


# TICC の手法の概要
　TICC は時系列データのセグメントとクラスタリングを同時に行う手法である．各クラスターごとに固有の**多層 MRF** を有すると仮定し，EM algorithm 的に最適解を探索して行く．

## Markov Random Field
　**Markov Random Filed (MRF)** とは，各特徴量をノードと捉え，特徴量間の関係をエッジとして捉えるグラフィカルモデルの一種である．特徴量$x_1,x_2$間に関係性がなければ，エッジの重みが0，すなわち$x_1,x_2$の間にエッジがないと考える．逆に，$x_1,x_2$間のエッジの重みが非0ならば，特徴量$x_1,x_2$間の関係性を重みで表す[^2]．

![mrf_graph.png](https://qiita-image-store.s3.amazonaws.com/0/218271/9d281315-0315-8082-661e-22bbcf603bf1.png)

[^2]: MRFの写真は， Kaikai Liu, Ting An, Hao Cai, Lirida Naviner, Jean-Francois Naviner, Herv´e Petit. A General Cost-effective Design Structure forProbabilistic-Based Noise-Tolerant Logic Functionsin Nanometer CMOS Technology. 2013. より拝借した．

　今回考える多層 MRF は次のように図式化される．各レイヤーが一つの時点における MRF を表しており，その時点における特徴量をノード，特徴量間の関係をエッジとして表している．また，時点を跨る特徴量間の関係を異なるレイヤー間のエッジで表している．

![スクリーンショット 2018-09-05 22.06.05.png](https://qiita-image-store.s3.amazonaws.com/0/218271/d23c5e6b-18d6-4b62-b176-156c612a3b1c.png)

## 学習概要
　最適化問題を解く際に，EM algorithm 的に E step, M step を収束するまで繰り返す形式となる．E step においてDP，M step において ADMM を使用していく[^3]．

[^3]: このように E step, M step を順次繰り返して最適化する手法を EM algorithm と呼ぶと思ってくれれば良い．


# 最適化問題の定式化
## 最適化問題のためのセットアップ
　長さ $T$ の時系列データ $\\{x_t\\}_{t=1}^T$ を考える．各時点の観測は $n$ 次元であるとし，これらを窓 $w$ 毎にまとめて横に並べたベクトル

$$X_t=(x_{t-w+1}^T,\cdots,x_t^T)\in\mathbb{R}^{nw}$$

を考える．構成としては，図でイメージするとわかりやすい．

![subclust.png](https://qiita-image-store.s3.amazonaws.com/0/218271/e8a946c8-f5d9-9673-f3cd-15901c7bb315.png)

　ここから最適化問題をセットアップしていく．クラスター数を $K$ とし，各クラスターは固有の精度行列 $\Theta_i\in\mathbb{R}^{nw\times nw}$ を有するとする．精度行列とは，ガウス型逆共分散行列のことであるが，MRF のネットワークの繋がりを記述したものである．例えば，$(j,k)$ 要素が0でなければ，ネットワークにおいては，変数 $x_j,x_k$ の間にエッジがあるということである．すなわち，クラスター毎にネットワーク構造が異なり，そのネットワーク構造の違いをもとにクラスタリングしようというのである．ここで，$\Theta_i$ はブロック Toeplitz 行列になっており，

```math
\Theta_i=\left[\begin{array}{cccccc}
A^{(0)}&(A^{(1)})^T&(A^{(2)})^T&\cdots&\cdots&(A^{(w-1)})^T\\
A^{(1)}&A^{(0)}&(A^{(1)})^T&\ddots&&\vdots\\
A^{(2)}&A^{(1)}&\ddots&\ddots&\ddots&\vdots\\
\vdots&\ddots&\ddots&\ddots&(A^{(1)})^T&(A^{(2)})^T\\
\vdots&&\ddots&A^{(1)}&A^{(0)}&(A^{(1)})^T\\
A^{(w-1)}&\cdots&\cdots&A^{(2)}&A^{(1)}&A^{(0)}
\end{array}\right]
```

と表せる．ただし，$A^{(j)}\in\mathbb{R}^{n\times n}$ である．ブロック Toeplitz 行列とは，このようにブロック対角が全て同じで，ブロック二重対角も（転置を除いて）同じといったように，対角からの外れ具合で同じ値を取るブロック行列のことを指す[^4]．

[^4]: Toeplitz 行列と言った場合は，対角成分，二重対角成分が全て同じ行列のことを言う．空間的・時間的距離に応じて値を一緒にしてパラメータ数を減らしたい場合などに用いられる．

![image.png](https://qiita-image-store.s3.amazonaws.com/0/218271/040368f6-c823-f3ba-5ed6-315916505f1e.png)

MRF をイメージするために再度，上の画像を見て欲しい．この場合は，window size $w=3$ であり，

```math
\Theta_i=\left[\begin{array}{ccc}
A^{(0)}&(A^{(1)})^T&(A^{(2)})^T\\
A^{(1)}&A^{(0)}&(A^{(1)})^T\\
A^{(2)}&A^{(1)}&A^{(0)}
\end{array}\right]
```

となる．$A^{(0)}$ はそれぞれのレイヤー（時点）におけるエッジを表している．例えば，画像の青いレイヤーに注目した時，どの変数同士が繋がっているかを表している．これが対角に並んでいるということは，紫や橙のレイヤーでも同じ変数同士の繋がりを持っているということを意味する．$A^{(1)}$ は隣のレイヤーの変数との繋がりを表している．例えば，画像の青・紫レイヤーに注目したときに，青レイヤーの変数と紫レイヤーの変数の繋がりを表している．$A^{(2)}$ 以降も同様に，２つ以上離れたレイヤーの変数間の繋がりを表現している．

　さて，最適化問題のセットアップに戻ろう．各クラスター $i$ に関して割り当て集合 $P_i$ を考える．例えば，$X_2,X_3,X_7$ がクラスター $i$ に属している場合，$P_i=\{2,3,7\}$ となる．集合用語で述べると，$\{1,\cdots,T\}$ を $P_1,\cdots,P_K$ に直和分解することになる．


## 定式化
　では，最適化問題を定式化しよう．最適化問題は次式のようになる：

```math
\mathrm{arg}\min_{\Theta\in\mathcal{T},P}\sum_{i=1}^K\left[\|\lambda\circ\Theta_i\|_1+\sum_{X_t\in P_i}\left(-ll(X_t,\Theta_i)+\beta\mathbf{1}\{X_{t-1} \notin P_i\}\right)\right].
```

これを **TICC 問題**と呼ぶことにする．第１項はスパースな MRF を得るための l1 norm である．第２項は，

$$ll(X_t,\Theta_t)=-\frac{1}{2}(X_t-\mu_i)^T\Theta_i(X_t-\mu_i)+\frac{1}{2}\log\det\Theta_i-\frac{n}{2}\log(2\pi)$$

で記述できる対数尤度関数である．正則化項がなければ，この項を最小化する問題となる．第３項は，temporal consistency と呼ばれる項であり，前の時点と別のクラスターに遷移するのに罰則 $\beta$ を課している．すなわち，クラスターが変わりすぎないようにするための項である．


# 部分問題
　TICC 問題は直接最適化するのが困難である．そこで，幾つかの部分問題に分割し，それらを順番に解くことで $P,\Theta$ の最適解を見つけていく．

## クラスター割り当て（E step）
　クラスター割り当てにおいては，部分問題

$$\mathrm{arg}\min_{P}\sum_{i=1}^K\sum_{X_t\in P_i}\left(-ll(X_t,\Theta_i)+\beta\mathbf{1}\{X_{t-1} \notin P_i\}\right)$$

を解く．この問題は，**動的計画法（DP）**を用いれば，大域最適解を見つけることができる．動的計画法とは，あるコストを最小にするパスを見つける際に用いられるアルゴリズムである．具体的には，途中時点では最もコストの小さいパスだけを保存し，最後にコスト最小パスをバックワードで見つけるものである．

![image.png](https://qiita-image-store.s3.amazonaws.com/0/218271/765395c9-d4b5-4a99-3c1e-c0f6b454a08a.png)

画像を見て欲しい． 例えば，$-ll(2,1)$ に行くには $K$ 通りのパスがある．$-ll(1,1)$ から来た場合のコストは $-ll(1,1)$，$-ll(1,2)$ から来た場合のコストは $-ll(1,2)+\beta$ といった具合になる．最終時点に到達したときに最小コストのパスを見つけたいので，途中時点においてもその時点に至るまでのパスの中で最小コストのものだけ保存しておけば十分だろう．


## 精度行列の更新（M step）
　精度行列の更新においては，次の部分問題（TGL; Toeplitz Graphical Lasso 問題）

```math
\mathrm{arg}\min_{\Theta_i\in\mathcal{T}}-\log\det\Theta_i+\mathrm{tr}(S_i\Theta_i)+\frac{1}{|P_i|}\|\lambda\circ\Theta_i\|_1
```

を各 $\Theta_i$ について解けば良い（各 $\Theta_i$ は並列に解くことができるので，以後簡単のため $\Theta$ と表記する）．ここで，$S_i$ はクラスター $i$ の標本共分散である．この形の問題に対しては，ADMM のような手法が有効であることが知られている．具体的には，$\Theta=Z$ なる $Z$ を導入してラグランジュ型の問題に落とし込んで行く：

```math
\mathrm{minimize\quad}-\log\det\Theta+\mathrm{tr}(S\Theta)+\|\lambda\circ Z\|_1,\\
\mathrm{subject\ to\quad}\Theta=Z,Z\in\mathcal{T}.
```

ここで，$\mathcal{T}$ は $nw\times nw$ Toeplitz 行列全体を示す．こうすることで，拡張ラグランジュ法が適用可能になる．拡張ラグランジュ法は，次のように制約条件の二乗ノルムを最適化問題に敢えて加えることで，最適化が容易になるテクニックだと思ってもらえれば良い（詳細は最後に述べる文献参照）．

```math
\mathcal{L}_p(\Theta,Z,U)=-\log\det\Theta+\mathrm{tr}(S\Theta)+\|\lambda\circ Z\|_1+\frac{\rho}{2}\|\Theta-Z+U\|_F^2
```

ここで，$Z$ は consensus 変数，$\rho>0$ は ADMM 罰則パラメータ，$U\in\mathbb{R}^{nw\times nw}$ は scaled 双対変数と呼ばれている．

ADMM はこのように拡張ラグランジュ法を使った手法の一つと思ってもらえればよく，この定式化によって，TGL 部分問題は次の３ステップを収束するまで繰り返せば良い．

```math
\mathrm{(a)}\quad\Theta^{k+1}:=\mathrm{arg}\min_\Theta\mathcal{L}_p(\Theta,Z^k,U^k)\\
\mathrm{(b)}\quad Z^{k+1}:=\mathrm{arg}\min_{Z\in\mathcal{T}}\mathcal{L}_p(\Theta^{k+1},Z,U^k)\\
\mathrm{(c)}\quad U^{k+1}:=U^k+(\Theta^{k+1}-Z^{k+1})
```

$\Theta, Z$ の更新も閉形式（cos, exp などのよく知られた関数で記述できる形式）で書くことができ，陽に各ステップ計算できる．また，TGL 部分問題は凸問題であることから，ADMM は大域解への収束が保証されている．


# 実験結果
## Macro F1-score
　**F1-score** とは，precision と recall の調和平均で定義される．クラスター $i$ の予測に関して，precision とは，クラスター $i$ に属すると予測したもののうち実際に属していた割合である．一方，recall とは，クラスター $i$ に属するもののうちクラスター $i$ と予測された割合である．一般に，precision と recall はトレードオフの関係であり，F1-score を用いることがある．この F1-score を各クラスターに対して計算し，その平均を取ったものが **Macro F1-score** である．Macro F1-score は高い方が良く，0以上1以下である．
　擬似生成したデータにおいて，各手法で Macro F1-score を比べた結果が次の表である．Temporal Sequence にある「1,2,1」とは，正解系列が1,2,1のクラスター遷移をしているということである．ここで，１つのクラスターに滞在している時点数は$100\times$（クラスター数）である．

![image.png](https://qiita-image-store.s3.amazonaws.com/0/218271/44d0e5c6-5f3e-12d4-0e56-24dec5013d6e.png)

ここで，GMM は混合ガウス法，EEV は正則化 GMM，DTW は動的時間伸縮法，GAK,Euclidean はそれぞれ global alignment kernel，Euclidean distance metric，Neural Gas は自己組織化マップ，K-means は Kmeans 法である．また，1クラスターに滞在するサンプル数を変化させた時の Macro F1-score の遷移は次のようになったと述べられている．

![image.png](https://qiita-image-store.s3.amazonaws.com/0/218271/d9bdabc0-6dac-790e-497f-1a73e4528ac7.png)

TICC が他の手法を上回っていることがわかる．

## Case Study
 Case Study として，車の走行状態のクラスタリングが示されている．7変数（ブレーキ，二次元加速度，方向，速度，RPM，ガソリン）を用いて，クラスター数5つとして TICC を適用したところ表のように，各クラスターの解釈が容易であった．また，GPS 位置に対して，クラスターを色で当てはめたところ，ブレーキを踏んでいる状態や曲がっている状態等のクラスタリングに成功している．

![image.png](https://qiita-image-store.s3.amazonaws.com/0/218271/c800043b-f7ee-442a-08b6-9ce17ee88787.png)

![image.png](https://qiita-image-store.s3.amazonaws.com/0/218271/92d30d79-b9bf-2897-dad8-d8d608d18045.png)

# 終わりに
　今回は TICC という時系列クラスタリング手法について述べてきた．時系列クラスタリングは，確立した手法が少なく，未だ未だ発展可能性を秘めている分野である．今回の手法は，https://github.com/davidhallac/TICC に公開されているので，ぜひ試してみてほしい．

# 参考文献
* D. Hallac, S. Vare, S. Boyd, J. Leskovec. Toeplitz Inverse Covariance-Based Clustering of Multivariate Time Series Data. in KDD, 2017. https://cs.stanford.edu/~jure/pubs/ticc-kdd17.pdf
    * 今回の記事の原論文である．

* S. Zolhavarieh, S. Aghabozorgi, and Y. W. Teh. A review of subsequence time series clustering. 2014. https://www.researchgate.net/publication/263470292_A_Review_of_Subsequence_Time_Series_Clustering.
    * 時系列クラスタリングの近年のレビュー論文である．

* N. Begum, L. Ulanova, J. Wang, and E. Keogh. Accelerating dynamic time warping clustering with a novel admissible pruning strategy. In KDD, 2015. http://www.cs.ucr.edu/~eamonn/Speeded%20Clustering%20Paper%20Camera%20Ready.pdf
    * DTW を用いた時系列クラスタリング手法が述べられている．

* C.M.ビジョップ著，元田浩・粟田多喜夫・樋口知之・松本裕治・村田昇監訳，パターン認識と機械学習 下 ベイズ理論による統計的予測，丸善出版，2008．
    * 言わずと知れたビジョップ本．MRF 等の勉強には良いかもしれない．

* 山下信雄，非線形計画法 (応用最適化シリーズ6)，朝倉書店，2015．
    * 拡張ラグランジュ法まで述べられており，最適化の入門書である．

* 大関真之，今日からできるスパースモデリング，http://www-adsys.sys.i.kyoto-u.ac.jp/mohzeki/Presentation/lecturenote20150902.pdf
    * ADMM 等がわかりやすく述べられている．