位相的データ解析超入門 -Persistence Diagramを計算してみる-

16945 ワード
TDA 位相的データ解析 Python データ分析 力学系 Python テキストリンク

この記事では位相的データ解析(Topological Data Analysis)という技術について，中でも重要な要素であるpersistence diagram(パーシステント図)の概略と，パッケージを使った計算についてまとめてみたいと思います。

位相的データ解析とは？

位相的データ解析 ( Topological Data Analysis; TDA) は，Topologicalと名のつくようにデータの幾何的な性質を捉えることを目的とした手法です。例えば球体があったとき，その球体に空洞があるのかそうでないか，空洞があればその大きさはどのくらいかということを定量的に測ることができます。一見すると大したことが無いように思えますが，形状が複雑になればなるほど，この手法は威力を発揮します。

特にTDAの中で状態を定量的に捉える指標として，persistence diagramがあります。これは点群が幾何的に為すロバストな性質を捉える際に有用になります。

どんなところに使われる？

TDAは情報科学や生物学，物質科学の分野などの様々な分野で応用がなされています。
ここでは物質科学に適用された例として，液体とガラスの相転移点に関する話題を見てみましょう。

液体とガラスの分類

液体とガラスの境界は，古くから統計物理学の難問として取り組まれ続けています。というのも，このふたつの状態の境界は非常に曖昧で，一見同じ状態にしか見えないためです。

※ 例えば「位相的データ解析の基礎と応用」のp2を参照：https://www.wpi-aimr.tohoku.ac.jp/hiraoka_labo/introduction_j.pdf

ここで近年，東北大と統数研はpersistence diagramを応用して，${\rm SI O}_2$の相転移点を検出するということに取り組みました。物理的にはエンタルピー曲線を描き，微分が不連続となる点を推定することによってこの相転移点を求めるのですが，彼らの方法では，まず分子動力学法で計算して得られた点群に対して先のpersistence diagramを計算し，これをカーネル法を用いてベクトル表現として得た上で，カーネル変化点検出法と呼ばれる手法を用いて相転移点を検出しました。結果として，物理的な方法で推定される範囲内の点が相転移点として検出され，Topological Data Analysisによる手法が物質科学における解析でも有効であることを示しました。

Persistent HomologyとPersistence Diagram

じゃあそのpersistent diagramってどうやって計算するの？となりますね。
ここで計算の考え方について簡単に記述したいと思います。

まず次の集合$X_r$を考えます。

$$X = \{x_1,\cdot\cdot\cdot , x_n \}$$$$B(x_i;r) = \{x \in M | d_M(x_i, x) \leq r)\}$$$$X_r = \cup_{i=1}^n B(x_i;r)$$

この集合$X_r$は半径$r$を持ったボールからなる集合です。ここである半径$r$のときに，ボールによって囲まれた穴が生じたとします。このとき，この集合を代数的に捉えるためにまず鎖群として表現し，境界作用素と呼ばれる図形の連結状態を調べる作用素を施し計算することによって，代数的に連結成分の数や穴のあるなしを判定することができます。その結果として，穴や空洞に関する情報を持ったホモロジー群$H_q$を得ることができます($q$は境界作用素の次元)。

※ 鎖群とは頂点や辺，あるいは三角形からなる集合です。頂点からなる集合，辺からなる集合，三角形からなる集合の順に高次元になっていきます。

※ 鎖群に対して境界作用素$\partial_q$を作用させると，境界を取り出す操作になります。例えば頂点$|v_0|, |v_1|$からなる辺$|v_0v_1|$に対して境界作用素$\partial_1$を作用させると$\partial_1|v_0v_1|:=|v_0|+|v_1|$となり，辺の境界として頂点を取り出すことができます。また三角形がつくる面に対して境界作用素を施せば，境界として辺が3つ取り出すことができます。詳しくは下記を参照のこと。

ここで$r$を変化させればそれに対して位相状態は変化していき，異なるホモロジー群を得ることができます。そして$r$の増大の中で得られるホモロジー群$H_q$の列をpersistent homology(以下PH)と呼びます。

ところで球体の連結によって作り出される穴は，半径の増大の中で生成されたり消失したりしますね(図ではある半径で$\alpha_1$, $\alpha_2$が生成され，さらに大きな半径で$\alpha_2$が消失)。ここで穴が生まれる半径を生成半径(birth time), 消失する半径を消滅半径(death time)と呼びます。

こうして得られる各リングの生成半径と消滅半径を対にしてプロットすると図のようになり，これをpersistence diagram(以下PD)と呼びます。PDでは対角線の上側に点がプロットされることになりますが，このとき各点の対角線からの距離がすなわち半径の増大に対するリング等の寿命や頑健性を意味することになります。対角線に近い点ほど寿命が短くノイズ的な性質を持ち，遠い点ほど寿命が長くロバストな性質を持つことを意味します。

persistence diagramを計算してみる

ここまででPDに関する基本事項を確認してきました。ここからは，実際にライブラリを使って計算を行ってみましょう。

位相的データ解析を扱うライブラリはphatやdiphaなど多く存在しますが，今回は東北大学で開発されている，HomCloudというライブラリを利用してみます。HomCloudのインストール方法は下記に詳しいです。

まず点群を生成します。ここでは代表的な力学系であるLorenz方程式を用いて，点群を生成してみることにします。

Lorenz方程式

$$\frac{dx}{dt} = -px + py$$$$\frac{dy}{dt} = -xz + \gamma x -y$$$$\frac{dz}{dt} = xy - bz$$

Lorenz方程式はパラメータ$p, \gamma, b$を変化させることで異なる幾何状態が実現します。特に$p=10, r=28, b = 8/3$のケースは非常に有名です。

※ 半径$r$との混同を避けるために，一般にLorenz方程式が持つパラメータ$r$を$\gamma$としています。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Lorentz方程式のパラメータ
p, gamma, b = 10, 28, 8/3

# 時間パラメータ
dt = 0.01 # 刻み幅
t_0 = 0 #　初期時刻
t_n = 50 # 終了時刻

# 位置パラメータ
X_0 = np.array([1, 1, 1]) # 質点の初期位置


# ルンゲクッタ法に関する関数
def RungeKutta(t, X):
  k_1 = LorenzEquation(t, X)
  k_2 = LorenzEquation(t + dt/2, X + k_1*dt/2)
  k_3 = LorenzEquation(t + dt/2, X + k_2*dt/2)
  k_4 = LorenzEquation(t + dt, X + k_3*dt)

  return X + dt/6*(k_1 + 2*k_2 + 2*k_3 + k_4)

def LorenzEquation(t, X):
  x, y, z = X[0],  X[1],  X[2]

  return np.array([-p*x + p*y, -x*z + gamma*x - y, x*y - b*z])


if __name__ == '__main__':
    t = t_0
    X = X_0
    data = np.r_[X]

    while t < t_n:
      X = RungeKutta(t, X)
      t += dt
      data = np.c_[data, X]

    print(data)
    fig = plt.figure()
    ax = Axes3D(fig)

    ax.plot(data[0,:], data[1,:], data[2,:])
    plt.show()

    np.savetxt('/Users/shinoda/Downloads/homcloud-examples/pointcloud/input.txt', data, fmt='%.18e', delimiter=' ')

順にパラメータ$\gamma$を変化させていって見ましょう。

A. $p = 10, \gamma = 10, b = 8/3​$の時は

となり，またB. $p = 5, \gamma = 28, b = 4/3$のときは

で，さらにC. $p = 5, \gamma = 350, b = 4/3$のときは

となります。幾何的な特徴としては，螺旋構造によってAが数個のリングが重なったような形，Bはさらに複雑にリングが幾重にも重なったような形，Cはそれが崩壊してしまったような形になっています。

こうして生成された点群に対して，先ほどのHomCloudを用いてPDを計算してみます。下記サイトからhomcloud-examplesをダウンロードし，ディレクトリpointcloudの中に3次元データを記録したファイルをinput.txtとして保存をして，シェルスクリプトrun.shを

./run.sh

として実行すると，リングの情報を捉えた1次のPDをoutput.1.pngとして得ることができます。

Aのケースでは

またBのケースでは

さらにCのケースでは

となります。Aでは対角線の上側に，点が間隔を伸ばしながらプロットされていっているのがわかります。Bはよりロバストなリング構造を多数抽出しており，Cでは再びぱらぱらと細かなリングが抽出されているのがわかります。ただし，Aとは異なってリングが生じる間隔が乱れているのがわかります。

以上では用意されているシェルスクリプトにPDの生成まで任せましたが，さらに細かな解析をしたい場合はPDの情報がまとまっているpointcloud.idiagramというファイルを参照するとよいです。

※ シェルスクリプト実行時に一緒にディレクトリpointcloudの中に生成されています。

python3 -m homcloud.dump_diagram -S no -d 1 pointcloud.idiagram | tail

とすると

2948.0788842052393 2955.760470756709
2955.484855668927 2955.7618395347713
2892.180946856205 2955.7619854858694
2931.1911476159944 2955.7626685806863
2911.8643876078927 2955.762862465481
2910.2794456629676 2955.7630131635515
2955.1307368888242 2955.7643060783407
2955.6599549461243 2955.766776724913
2955.543271613065 2955.7669608942924
2955.7542468578076 2955.7698308547

のように中身を確認することができます。左側がbirth time, 右側がdeath timeになります。また先ほどのpython実行時の数字が1の場合は1次の，2の場合は2次のbirth death pairが出力されることになります。より詳しくはHomCloudのチュートリアルを参考にして下さい。

まとめ

今回は位相的データ解析について概略をまとめてみました。面白い概念・手法なので，今後の発展にさらに期待したいですね。ちなみにpersistence diagramの他の代表的な表現形式にbarcodeやpersistent landscapeなどがあるので，自分にあった形式を見つけて使って見てください。

参考