Persistence Codebook を使った機械学習

このチュートリアルは Persistence Image を使った機械学習の一通りを済ませたことを仮定して進めます。 重複した内容は説明しません。

import os  # for makedirs
import homcloud.interface as hc  # HomCloud 
import homcloud.interface.codebook  # HomCloud's implementation of persistence codebook
import numpy as np  # Numerical array library
from tqdm.notebook import tqdm  # For progressbar
import matplotlib.pyplot as plt  # Plotting
import sklearn.linear_model as lm  # Machine learning
from sklearn.decomposition import PCA  # for PCA
from sklearn.model_selection import train_test_split

ラベルを読み込んで表示します。

labels = np.loadtxt("pc/label.txt")
labels

array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
       1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])

ポイントクラウドデータの読み込み

pointclouds = [np.loadtxt("pc/{:04d}.txt".format(i)) for i in tqdm(range(200))]

  0%|          | 0/200 [00:00<?, ?it/s]

パーシステンス図の計算

それでは、HomCloud を使ってパーシステンス図を計算しましょう。
計算されたファイルは pd ディレクトリに保存しておきます。

os.makedirs("pd", exist_ok=True)
for i in tqdm(range(200)):
    hc.PDList.from_alpha_filtration(pointclouds[i], save_boundary_map=True, 
                                    save_to="pd/{:04d}.pdgm".format(i))

  0%|          | 0/200 [00:00<?, ?it/s]

このチュートリアルでは、2 次のパーシステンス図を使用します。つまり、空隙構造（ボイド）に注目して解析を行います。

pds = [hc.PDList("pd/{:04d}.pdgm".format(i)).dth_diagram(2) for i in tqdm(range(200))]

  0%|          | 0/200 [00:00<?, ?it/s]

0番目(ラベル0)と100番目(ラベル1)のパーシステンス図を表示します。

pds[0].histogram(x_range=(0, 0.03)).plot(colorbar={"type": "log"})
pds[100].histogram(x_range=(0, 0.03)).plot(colorbar={"type": "log"})

Persistence Codebookによるベクトル化

Persistence Codebook については、以下の文献を参考にしてください：

• https://link.springer.com/article/10.1007/s10462-020-09897-4

この論文では多くのベクトル化法が提案されていますが、HomCloud では PBoW, Stable PBoW, PFV の3つのみ実装されています。 ここでは PBoW を使いましょう。 このチュートリアルの作者 (大林) の経験では、PBoWが一番手軽に利用できます。 パラメーター調整を適当にやってもそれなりに良い結果が得られます。

Persistence codebook では、ベクトル化の対象となるPDのbirth-death pairをすべて混ぜ合わせたものからランダムにサンプリングし、そこからpairのクラスターを計算します。 クラスタリングによってPDのpairの分布を低次元化することを目指します。 基本的にサンプリングはクラスタリングのコストを下げるために利用しますが、それだけでなく各pairを重み付けてサンプリングすることでpairの重要さを反映させる効果があります。

次のコードで、サンプリングするpairの個数(100000個)、クラスタ数 (20個)、重み関数 ($\arctan(0.01 (d - b)^3)$)を指定します。

def weight(b, d):
    return np.atan(0.01 * (d - b) ** 3)

vectorize_spec = hc.codebook.PBoWSpec(100000, 20, weight)

次のセルでクラスタを計算します。

vectorize_spec.fit(pds)

ではパーシステント図をベクトル化しましょう。

pdvects = np.vstack([vectorize_spec.vectorize(pd) for pd in pds])

主成分分析（PCA）

まずは PCA（主成分分析）を用いてデータを解析します。
PCA は、高次元のデータをできるだけ情報を失わずに低次元に圧縮するための手法です。
解析結果が大きく破綻することが少ないため、パーシステンス図と機械学習の組み合わせでは第一に試す価値があります。

パーシステンス図との組み合わせでは、各図を低次元の点として表現できるため、図同士の関係性を視覚的に理解しやすくなります。

ここでは、可視化しやすいように 2 次元に次元削減を行います。

pca = PCA(n_components=2)
pca.fit(pdvects)

PCA(n_components=2)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

どのように2次元上に表現されているか可視化します。 赤がラベル0、青がラベル1を意味します。

reduced = pca.transform(pdvects)  # すべてのデータを2次元に落とす
plt.gca().set_aspect('equal')  # 縦横のアスペクト比を揃える
plt.scatter(reduced[labels == 0, 0], reduced[labels == 0, 1], c="r")  # ラベル0のデータを赤で描画
plt.scatter(reduced[labels == 1, 0], reduced[labels == 1, 1], c="b")  # ラベル1のデータを青で描画

<matplotlib.collections.PathCollection at 0x7901e5904cb0>

2 次元平面上でX軸方向で赤と青の点がある程度分離されていることから、PBoWによって得られたベクトルは2つのグループの違いをうまく表現できていると考えられます。

では、この図の X 軸と Y 軸は何を表しているのでしょうか？

X 軸と Y 軸に対応するベクトルをそれぞれ $v_1$, $v_2$ とすると、PBoWによってベクトル $w_i$（$i$ はデータのインデックス）は次のように近似されます：

$$ w_i \approx \mu_{i1} v_1 + \mu_{i2} v_2 + v_0 $$

ここで $v_0$ は $w_1, w_2, \ldots$ の平均ベクトルで，$(\mu_{i1}, \mu_{i2})$が上の図の各点の座標です（このあたりの理論については PCA の基本を学んでください）。
つまり、$v_1$ と $v_2$ がどのようなベクトルであるかがわかれば、この図が何を意味しているかをより深く理解できます。

なお、scikit-learn の PCA では、$v_0$ は PCA.mean_、$v_1$, $v_2$ は PCA.components_ に格納されています。

pca.mean_

array([1382.53977497, 1113.12157663, 2316.20804098,  544.58056669,
        921.26733574, 1706.82253307, 1050.23883135, 1234.84060285,
        944.74688089, 2418.6258216 , 1530.33987473,  246.66255163,
        740.14725272,  349.92313533,  272.04492509,  355.24084419,
       1411.41960292, 1562.27499827,   62.13692049, 1643.1220043 ])

pca.components_

array([[ 0.39326077, -0.24055517,  0.21960124, -0.07153537, -0.20310608,
         0.48443186, -0.00377261, -0.06989213, -0.1630181 ,  0.18469233,
        -0.19420942, -0.05430333, -0.0748403 ,  0.00398742, -0.04671522,
        -0.04850072,  0.12773821,  0.09995429, -0.01214863,  0.56156237],
       [ 0.10362294,  0.23225009,  0.14514279,  0.09262878, -0.07614634,
        -0.41993717,  0.11772431,  0.16389753, -0.19611693,  0.21991956,
        -0.31901628, -0.00838165,  0.12423459,  0.09569389, -0.02517399,
         0.0405355 ,  0.176744  ,  0.65894105, -0.00131564, -0.04276406]])

各ベクトルの次元は20で、これはクラスター数と同じです。 Persistence codebook では、これらの係数はクラスタごとの重みとして解釈されます。 つまりクラスタリング結果を理解するためにはクラスタリング結果を調べるのが肝心です。

そこで、クラスタリングの中心をプロットしてみましょう。vectorize_spec.cluster_centers で各クラスターの中心座標が得られます。

vectorize_spec.cluster_centers

array([[0.01011662, 0.01670536],
       [0.00829549, 0.01951828],
       [0.0069341 , 0.0104699 ],
       [0.01272092, 0.02197064],
       [0.00880535, 0.02163759],
       [0.00746409, 0.01755715],
       [0.0148032 , 0.01946074],
       [0.01019561, 0.01915189],
       [0.00729412, 0.02075994],
       [0.00610172, 0.00811547],
       [0.00704536, 0.01882478],
       [0.00747822, 0.02268904],
       [0.01070353, 0.020874  ],
       [0.01520447, 0.02324038],
       [0.01036722, 0.0244698 ],
       [0.01825379, 0.0231114 ],
       [0.01239476, 0.01853661],
       [0.00885005, 0.01819233],
       [0.0079025 , 0.02517447],
       [0.00730205, 0.01608447]])

散布図でクラスター中心をプロットします。PDと比較しやすいよう対角線も描画します。

import matplotlib.lines

fig, ax = plt.subplots()
ax.scatter(vectorize_spec.cluster_centers[:, 0], vectorize_spec.cluster_centers[:, 1])
ax.set_aspect("equal")
ax.add_artist(matplotlib.lines.Line2D([0, 0.03], [0, 0.03], color="k"))
ax.set_xlim(0, 0.03)
ax.set_ylim(0, 0.03)

(0.0, 0.03)

death > 0.015 の領域に大半のクラスタ中心が分布していますが、death < 0.015の領域にも3箇所あります。

PDを重ね合わせたものをクラスタリングする際には対角線を割り引いてサンプリングしているので、対角線のそばには多くのbirth-death pairが分布しているにもかかわらずクラスタ中心が現われません。 この重みの調整はPBoWの運用において肝心な点です。

次にPCAの結果を可視化します。

PCAの結果で計算されるベクトル$v_i$の次元はクラスター数と同じで、このベクトルの各要素がクラスター毎の「重み」と見なせます。 詳しくはCodebookの論文を読んでください。 そこで、$v_i$の各要素をクラスター中心の点の色で表現しましょう。 その値が正なら青、負なら赤、0に近ければ白、となるように色付けを調整します。

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

for k in [0, 1]:
    ax = axes[k]
    sc = ax.scatter(
        vectorize_spec.cluster_centers[:, 0], 
        vectorize_spec.cluster_centers[:, 1],
        c=pca.components_[k, :],
        vmax=0.7,
        vmin=-0.7,
        cmap='RdBu'
    )
    
    cbar = fig.colorbar(sc, ax=ax)
    ax.grid(True)
    ax.set_aspect("equal")
    ax.add_artist(matplotlib.lines.Line2D([0, 0.03], [0, 0.03], color="k"))
    ax.set_xlim(0, 0.03)
    ax.set_ylim(0, 0.03)

# グラフの表示
plt.show()

この色の分布を見てPCAの結果を解釈します。 解釈法はロジスティック回帰と似ているので、そちらの説明を見てください。

ロジスティック回帰

次はロジスティック回帰です。ロジスティック回帰は、PCA に続いてパーシステントホモロジーとの組み合わせが比較的容易な手法です。
特に2クラス分類が解釈しやすく、扱いやすいという利点があります。

まずは学習を始める前にデータをテスト用（test set）と学習用（training set）に分割します。

pdvects_train, pdvects_test, labels_train, labels_test = train_test_split(pdvects, labels, test_size=0.25)

学習用データを使ってモデルを学習します。
ここでは、sklearn の LogisticRegression クラスを使用します。
正則化パラメータ C は、今回はあらかじめ適当な値を設定していますが、LogisticRegressionCV クラスを使えば、交差検証によって最適なパラメータを自動的に選ぶことも可能です。

model = lm.LogisticRegression(C=0.01, solver="lbfgs")

model.fit(pdvects_train, labels_train, )

LogisticRegression(C=0.01)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

まずテスト用データを使って性能を確認します。

model.score(pdvects_test, labels_test)

0.98

高い精度です。 次にロジスティック回帰の係数を以下で可視化します。

まずは係数の数値を見てみましょう。

model.coef_

array([[-0.01049108,  0.00159178, -0.00790289,  0.0028876 ,  0.00148073,
        -0.00490653,  0.00829817, -0.00032661,  0.00164977,  0.00560674,
         0.00312246,  0.00093899,  0.00170018,  0.00687491, -0.00083505,
         0.00625443,  0.00148933, -0.00139073,  0.00244847, -0.00490393]])

model.coef_.max(), model.coef_.min()

(np.float64(0.008298170615342543), np.float64(-0.010491077669949847))

fig, ax = plt.subplots()

sc = ax.scatter(
    vectorize_spec.cluster_centers[:, 0], 
    vectorize_spec.cluster_centers[:, 1],
    c=model.coef_,
    vmax=0.01,
    vmin=-0.01,
    cmap='RdBu'
)

cbar = fig.colorbar(sc, ax=ax)
ax.grid(True)
ax.set_aspect("equal")
ax.add_artist(matplotlib.lines.Line2D([0, 0.03], [0, 0.03], color="k"))
ax.set_xlim(0, 0.03)
ax.set_ylim(0, 0.03)

# グラフの表示
plt.show()

この図の見方ですが、濃い赤の領域に多くの birth-death ペアが存在するデータは、ラベル 0 である可能性が高く、 一方で、濃い青の領域に多くの birth-death ペアが存在するデータは、ラベル 1 である可能性が高いことを示しています。

この図ではdeathが0.018あたりに青と赤の境界があるように見えます。つまりdeath=0.018のラインより下にpairが多ければラベル0のデータ、上にpairが多ければラベル1のデータであることがわかります。

PDをベクトル化して得られたベクトルとこのベクトルのアダマール積を取ると、どのクラスターがどの程度予測に貢献しているかもわかります。

そこで、ラベル0とラベル1のデータを1つづつ(1番目と101番目)取ってきて計算し、それを上のように可視化しましょう。

pdvects[1, :] * model.coef_

array([[-3.21191502e+01,  0.00000000e+00, -1.68964610e+01,
         7.24862167e+00,  0.00000000e+00, -1.94181397e+01,
         3.86970391e+00, -1.30557095e+00,  0.00000000e+00,
         1.14009952e+01,  1.52847312e+01,  0.00000000e+00,
         0.00000000e+00,  0.00000000e+00, -0.00000000e+00,
         1.57046035e-06,  1.69002750e+00, -0.00000000e+00,
         0.00000000e+00, -1.59528157e+01]])

pdvects[101, :] * model.coef_

array([[-2.83204956e-01,  5.97278013e+00, -8.49660910e+00,
         0.00000000e+00,  3.28045154e+00, -8.53542659e+00,
         4.65867559e+00, -0.00000000e+00,  0.00000000e+00,
         5.79111333e+00,  6.55395527e+00,  0.00000000e+00,
         2.94293212e+00,  0.00000000e+00, -0.00000000e+00,
         1.35214975e-03,  7.65964220e-04, -3.05768193e+00,
         0.00000000e+00, -0.00000000e+00]])

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

sc = axes[0].scatter(
    vectorize_spec.cluster_centers[:, 0], 
    vectorize_spec.cluster_centers[:, 1],
    c=pdvects[1, :] * model.coef_,
    vmax=30,
    vmin=-30,
    cmap='RdBu'
)

cbar = fig.colorbar(sc, ax=axes[0])
axes[0].grid(True)
axes[0].set_aspect("equal")
axes[0].add_artist(matplotlib.lines.Line2D([0, 0.03], [0, 0.03], color="k"))
axes[0].set_xlim(0, 0.03)
axes[0].set_ylim(0, 0.03)
pds[1].histogram((0, 0.03), 128).plot(colorbar={"type": "log"}, ax=axes[1])
fig.tight_layout()
# グラフの表示
plt.show()

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

sc = axes[0].scatter(
    vectorize_spec.cluster_centers[:, 0], 
    vectorize_spec.cluster_centers[:, 1],
    c=pdvects[101, :] * model.coef_,
    vmax=10,
    vmin=-10,
    cmap='RdBu'
)

cbar = fig.colorbar(sc, ax=axes[0])
axes[0].grid(True)
axes[0].set_aspect("equal")
axes[0].add_artist(matplotlib.lines.Line2D([0, 0.03], [0, 0.03], color="k"))
axes[0].set_xlim(0, 0.03)
axes[0].set_ylim(0, 0.03)
pds[101].histogram((0, 0.03), 128).plot(colorbar={"type": "log"}, ax=axes[1])
fig.tight_layout()
# グラフの表示
plt.show()

左右を見比べると、どのあたりのpairが予測に効いているかがわかります。