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 ディレクトリに保存しておきます。
if not os.path.isfile("pd/0000.pdgm"): # 計算結果のファイルがある場合は計算をスキップする
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))
このチュートリアルでは、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"})
テスト用データと訓練用データの分割¶
後で性能評価をするために、データをテスト用と訓練用に分割しておきましょう。
from sklearn.model_selection import train_test_split
pds_train, pds_test, pcs_train, pcs_test, labels_train, labels_test = train_test_split(
pds, pointclouds, labels, test_size=0.25, stratify=labels, random_state=424231
)
Persistence Codebookによるベクトル化¶
Persistence Codebook については、以下の文献を参考にしてください：
この論文では多くのベクトル化法が提案されていますが、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_train)
ではパーシステント図をベクトル化しましょう。
pdvects_train = np.vstack([vectorize_spec.vectorize(pd) for pd in pds_train])
pdvects_test = np.vstack([vectorize_spec.vectorize(pd) for pd in pds_test])
ベクトル化パラメータの保存法¶
Persistence codebook ではパーシステンス図の集合に対し適応的にベクトル化のパラメータを決めます。
これは性能向上などに有効ですが、一方以下のような欠点があります。
- ベクトル化パラメータの調整に乱数を用いるため、コードを実行するたびパラメータが変化してしまう (これは random_state を固定することで解決可能)
- あるパーシステンス図の集合で調整したパラメータを他のパーシステンス図に利用することが単純にはできない。
これを解決する一つの方法に、vectorize_spec.fitによってパラメーターが調整済みのベクトル化オブジェクト(hc.codebook.PBoWSpecなどのインスタンス)をjoblib.dumpで保存しておく、という手段があります。
https://scikit-learn.org/0.15/modules/model_persistence.html の文章が参考になるでしょう。
以下のセルはそれを実行しています。
import joblib
joblib.dump(vectorize_spec, "pbowspec.joblib") # "pbowspec.joblib" というファイルに保存
['pbowspec.joblib']
backup_vectorize_spec = joblib.load("pbowspec.joblib") # "pbowspec.joblib" から読み出す
本当にうまくいっているのか確認しましょう。まずは読み出したオブジェクトのクラスを確認します。
type(backup_vectorize_spec)
homcloud.interface.codebook.PBoWSpec
homcloud.interface.codebook.PBoWSpec になっているので問題なさそうです。では再生したオブジェクトでのベクトル化の結果が元のオブジェクトでの結果と一致しているか確認しましょう。
np.all(pdvects_train == np.vstack([backup_vectorize_spec.vectorize(pd) for pd in pds_train]))
np.True_
大丈夫そうです。
この方法に関する注意点を述べておきます。joblib.dumpの出力はPythonのバージョンや関連するライブラリ(joblib自体やNumPyなど)のバージョン、またHomCloud自体のバージョンの変化に影響されやすいため、バージョンアップの前に保存したパラメータ情報をバージョンアップ後に読み出すことができない場合があります。
そういう場合は最初からやりなおしてください。
主成分分析（PCA）¶
まずは PCA（主成分分析）を用いてデータを解析します。
PCA は、高次元のデータをできるだけ情報を失わずに低次元に圧縮するための手法です。
解析結果が大きく破綻することが少ないため、パーシステンス図と機械学習の組み合わせでは第一に試す価値があります。
パーシステンス図との組み合わせでは、各図を低次元の点として表現できるため、図同士の関係性を視覚的に理解しやすくなります。
ここでは、可視化しやすいように 2 次元に次元削減を行います。
pca = PCA(n_components=2)
pca.fit(pdvects_train)
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.
Parameters
どのように2次元上に表現されているか可視化します。 赤がラベル0、青がラベル1を意味します。
reduced = pca.transform(pdvects_train) # すべてのデータを2次元に落とす
plt.gca().set_aspect('equal') # 縦横のアスペクト比を揃える
plt.scatter(reduced[labels_train == 0, 0], reduced[labels_train == 0, 1], c="r") # ラベル0のデータを赤で描画
plt.scatter(reduced[labels_train == 1, 0], reduced[labels_train == 1, 1], c="b") # ラベル1のデータを青で描画
<matplotlib.collections.PathCollection at 0x7f92797c7650>
2 次元平面上でX軸方向で赤と青の点がある程度分離されていることから、PBoWによって得られたベクトルは2つのグループの違いをうまく表現できていると考えられます。
訓練データで計算された主成分ベクトルを用いてテストデータを低次元化するとどうなるか見てみましょう。
reduced_test = pca.transform(pdvects_test)
plt.gca().set_aspect('equal')
plt.scatter(reduced_test[labels_test == 0, 0], reduced_test[labels_test == 0, 1], c="r") # ラベル0のデータを赤で描画
plt.scatter(reduced_test[labels_test == 1, 0], reduced_test[labels_test == 1, 1], c="b") # ラベル1のデータを青で描画
<matplotlib.collections.PathCollection at 0x7f92787c4590>
学習(ここでは低次元特徴量の抽出方法の推定)結果をテスト用データに適用しても整合性がありそうなのがわかります。
では、この図の 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([1664.91352175, 2052.88524489, 1245.19091159, 522.47288541,
1151.1411731 , 284.66273806, 1467.01489227, 713.94170372,
1005.72547671, 1540.72881554, 1246.27042268, 953.07611447,
1492.50349112, 1418.41048829, 247.51906179, 833.88138663,
1587.26252753, 2048.74461734, 1354.886933 , 1470.5984355 ])
pca.components_
array([[ 0.49555303, 0.19325772, -0.25708223, 0.02364196, 0.00865811,
-0.04409085, -0.19945591, -0.11946515, -0.22842562, 0.47035825,
0.1008718 , -0.09814304, 0.12447496, 0.16074127, -0.04677165,
-0.08585197, 0.21264832, 0.1663711 , -0.10836854, 0.40694011],
[-0.35303629, 0.12793067, 0.18751048, 0.10830663, 0.12472173,
0.00845479, -0.28741045, 0.17162585, -0.12331415, -0.10610985,
0.1047861 , 0.22077473, 0.19986818, 0.10915412, -0.02639036,
-0.24017012, 0.59929631, 0.17617607, 0.31257349, 0.00989052]])
各ベクトルの次元は20で、これはクラスター数と同じです。 Persistence codebook では、これらの係数はクラスタごとの重みとして解釈されます。 つまりクラスタリング結果を理解するためにはクラスタリング結果を調べるのが肝心です。
そこで、クラスタリングの中心をプロットしてみましょう。vectorize_spec.cluster_centers で各クラスターの中心座標が得られます。
vectorize_spec.cluster_centers
array([[0.00736558, 0.01741746],
[0.00677702, 0.00980954],
[0.00826278, 0.01959199],
[0.01760754, 0.02365729],
[0.01470264, 0.01957506],
[0.01034828, 0.02398021],
[0.00705663, 0.01884239],
[0.01342636, 0.02251998],
[0.00879975, 0.02161568],
[0.00721275, 0.01602792],
[0.0054994 , 0.00678362],
[0.01095228, 0.02072767],
[0.01230069, 0.01833465],
[0.00718031, 0.01146237],
[0.00754263, 0.02331428],
[0.00721567, 0.02074812],
[0.00851198, 0.01818002],
[0.0062774 , 0.00844851],
[0.00991327, 0.01901712],
[0.00994233, 0.01675596]])
散布図でクラスター中心をプロットします。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クラス分類が解釈しやすく、扱いやすいという利点があります。
学習用データを使ってモデルを学習します。
ここでは、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.
Parameters
まずテスト用データを使って性能を確認します。
model.score(pdvects_test, labels_test)
1.0
高い精度です。 次にロジスティック回帰の係数を以下で可視化します。
まずは係数の数値を見てみましょう。
model.coef_
array([[-0.00424358, -0.00667812, 0.00199262, 0.00837489, 0.00074102,
0.0011112 , 0.00296459, 0.00450787, 0.00232443, -0.00695285,
0.00024307, 0.00457056, 0.00416643, -0.00367759, 0.00328515,
0.00359114, -0.00070893, 0.00276709, 0.00123688, -0.01018475]])
model.coef_.max(), model.coef_.min()
(np.float64(0.008374886430854937), np.float64(-0.010184751935497826))
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つづつ取ってきて計算し、それを上のように可視化しましょう。
labels_test
array([0., 0., 0., 1., 1., 1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0., 1.,
1., 1., 1., 0., 1., 1., 0., 1., 0., 0., 1., 0., 1., 0., 1., 1., 1.,
0., 0., 0., 1., 1., 0., 1., 0., 0., 1., 0., 1., 0., 0., 1., 1.])
index_label_zero_in_test = np.nonzero(labels_test == 0)[0][0]
index_label_one_in_test = np.nonzero(labels_test == 1)[0][0]
index_label_zero_in_test, index_label_one_in_test
(np.int64(0), np.int64(3))
テスト用データの0番目と3番目を使いましょう。
pdvects_test[index_label_zero_in_test, :] * model.coef_
array([[ -0. , -19.97386745, 0. , 5.56786462,
0.54949104, 0. , 0. , 0. ,
0. , -29.54896137, 0.49892172, 0. ,
4.6030616 , -9.11525355, 0. , 0. ,
-4.05909444, 8.58478956, 6.87359142, -59.43787589]])
pdvects_test[index_label_one_in_test, :] * model.coef_
array([[-0. , -6.06966184, 0. , 3.52621256, 0.53958006,
4.03067426, 0. , 0. , 5.77210346, -0. ,
0.15145563, 0. , 4.12648274, -2.7530174 , 0. ,
13.0618278 , -1.50767409, 3.33156074, 2.69511444, -0.07347319]])
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_train[index_label_zero_in_test, :] * 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_test[index_label_zero_in_test].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_test[index_label_one_in_test, :] * 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_test[index_label_one_in_test].histogram((0, 0.03), 128).plot(colorbar={"type": "log"}, ax=axes[1])
fig.tight_layout()
# グラフの表示
plt.show()
左右を見比べると、どのあたりのpairが予測に効いているかがわかります。