機械学習チュートリアル (ポイントクラウド版)

まず利用するライブラリを読み込みます。Pythonの機械学習ライブラリ sklearn などを利用します。 sklearn については https://scikit-learn.org/stable/index.html を参考にしてください。その他のライブラリについても 勉強してください。

import os  # for makedirs
import homcloud.interface as hc  # HomCloud 
import homcloud.paraview_interface as pv # HomCloud <-> paraview interface
import numpy as np  # Numerical array library
from tqdm import tqdm_notebook as 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

%matplotlib inline
%load_ext autoreload
%autoreload 2

機械学習を始める前に

まずはデータがどんなものか確認していきます。チュートリアル用のデータは pc というディレクトリ以下にあります。

ls pc/

0000.txt  0026.txt  0052.txt  0078.txt  0104.txt  0130.txt  0156.txt  0182.txt
0001.txt  0027.txt  0053.txt  0079.txt  0105.txt  0131.txt  0157.txt  0183.txt
0002.txt  0028.txt  0054.txt  0080.txt  0106.txt  0132.txt  0158.txt  0184.txt
0003.txt  0029.txt  0055.txt  0081.txt  0107.txt  0133.txt  0159.txt  0185.txt
0004.txt  0030.txt  0056.txt  0082.txt  0108.txt  0134.txt  0160.txt  0186.txt
0005.txt  0031.txt  0057.txt  0083.txt  0109.txt  0135.txt  0161.txt  0187.txt
0006.txt  0032.txt  0058.txt  0084.txt  0110.txt  0136.txt  0162.txt  0188.txt
0007.txt  0033.txt  0059.txt  0085.txt  0111.txt  0137.txt  0163.txt  0189.txt
0008.txt  0034.txt  0060.txt  0086.txt  0112.txt  0138.txt  0164.txt  0190.txt
0009.txt  0035.txt  0061.txt  0087.txt  0113.txt  0139.txt  0165.txt  0191.txt
0010.txt  0036.txt  0062.txt  0088.txt  0114.txt  0140.txt  0166.txt  0192.txt
0011.txt  0037.txt  0063.txt  0089.txt  0115.txt  0141.txt  0167.txt  0193.txt
0012.txt  0038.txt  0064.txt  0090.txt  0116.txt  0142.txt  0168.txt  0194.txt
0013.txt  0039.txt  0065.txt  0091.txt  0117.txt  0143.txt  0169.txt  0195.txt
0014.txt  0040.txt  0066.txt  0092.txt  0118.txt  0144.txt  0170.txt  0196.txt
0015.txt  0041.txt  0067.txt  0093.txt  0119.txt  0145.txt  0171.txt  0197.txt
0016.txt  0042.txt  0068.txt  0094.txt  0120.txt  0146.txt  0172.txt  0198.txt
0017.txt  0043.txt  0069.txt  0095.txt  0121.txt  0147.txt  0173.txt  0199.txt
0018.txt  0044.txt  0070.txt  0096.txt  0122.txt  0148.txt  0174.txt  label.txt
0019.txt  0045.txt  0071.txt  0097.txt  0123.txt  0149.txt  0175.txt
0020.txt  0046.txt  0072.txt  0098.txt  0124.txt  0150.txt  0176.txt
0021.txt  0047.txt  0073.txt  0099.txt  0125.txt  0151.txt  0177.txt
0022.txt  0048.txt  0074.txt  0100.txt  0126.txt  0152.txt  0178.txt
0023.txt  0049.txt  0075.txt  0101.txt  0127.txt  0153.txt  0179.txt
0024.txt  0050.txt  0076.txt  0102.txt  0128.txt  0154.txt  0180.txt
0025.txt  0051.txt  0077.txt  0103.txt  0129.txt  0155.txt  0181.txt

0000.txt から 0199.txt まで 200個のテキストデータがあります。また、label.txtには 0 か 1 のラベルデータがあり、 これは対象のデータが2つのグループに分けられていることを示しています。

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

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.])

これは 200個の 0/1の配列です。では、どのデータが0/1でラベル付けされているのかを確認しましょう。 numpy の機能を使います。

np.nonzero(labels == 0), np.nonzero(labels == 1)

((array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
         17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
         34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
         51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
         68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84,
         85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99]),),
 (array([100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112,
         113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125,
         126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138,
         139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151,
         152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164,
         165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177,
         178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190,
         191, 192, 193, 194, 195, 196, 197, 198, 199]),))

0 から 99 までが 0 で、100から199までが 1 で、それぞれラベル付けられていることがわかります。 では、0番目のデータを100番目のデータを可視化してみましょう。homcloud.paraview_interfaceという3次元可視化のモジュールが HomCloudにはあるので使います。

pv.show([pv.PointCloud("pc/0000.txt", dim=3, delimiters=" ")], wait=False)

pv.show([pv.PointCloud("pc/0100.txt", dim=3, delimiters=" ")], wait=False)

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

ではファイルからポイントクラウドデータを読み込みます。

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

Widget Javascript not detected.  It may not be installed or enabled properly.

0番目のデータの「形」を見ると 1000x3 となっていることがわかります。これは3次元上の点が1000個あることを意味します。 どのポイントクラウドも同じです。パーシステント図による機械学習では点の個数は同じくらいでないとあまりうまくいかない (不可能ではないが注意点が多くなる)ことに注意しましょう。

pointclouds[0].shape

(1000, 3)

パーシステント図の計算

では、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))

Widget Javascript not detected.  It may not be installed or enabled properly.

このチュートリアルでは2次のパーシステント図を使います。つまり空隙構造に注目していきます。 もちろん1次や0次なども使えますが、実はこのデータは2次の所に一番特徴があるのでこうします。

より高度な解析法として、1次と2次を同時に使う方法もあります。

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

Widget Javascript not detected.  It may not be installed or enabled properly.

とりあえず計算したパーシステント図を表示してみます。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 Image (PI) でベクトル化

Persistence Image については https://arxiv.org/abs/1507.06217 http://www.jmlr.org/papers/volume18/16-337/16-337.pdf を参考にしてください。 ただし、ここで説明するやりかたはこの論文の手法そのままではなく、 https://arxiv.org/abs/1706.10082 https://link.springer.com/article/10.1007%2Fs41468-018-0013-5 で用いられている改変版です。

まずはベクトル化の仕様を決めます。 (0, 0.03), 128 で [0, 0.03]x[0, 0.03] という正方形内を128x128に分割して計算することを指定します。 sigma=0.002 はガウス分布の標準偏差を指定します。また、weight=("atan", 0.01, 3) で対角線からの距離による重みの変化を 指定します。この場合の重みを決めるルールは $\textrm{atan} (0.01 \cdot \ell^3)$ です。ただし$\ell$はlifetime(=birth-death)を意味 します。

vectorize_spec = hc.PIVectorizeSpec((0, 0.03), 128, sigma=0.002, weight=("atan", 0.01, 3))

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

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

ベクトルの要素の最大最小を見ると$10^{-10}$と非常に小さいです。これは機械学習に悪い影響を与えるので正規化しておきます。

pdvects.min(), pdvects.max()

(0.0, 2.099302720413663e-10)

pdvects = pdvects / pdvects.max()

主成分解析(PCA)

まずはPCAを使って解析します。PCAは高次元のデータをできるだけ情報を損わずに低次元に落とすための手法です。 あまり解析結果が破綻しないのでとりあえずパーシステント図との組み合わせでは試す価値があります。 パーシステント図との組み合わせにおいては、個々のパーシステント図を低い次元の点として表現でき、パーシステント図互いの 関係を理解しやすく表現します。

ここでは可視化しやすい2次元に落とします。

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

PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
    svd_solver='auto', tol=0.0, whiten=False)

さてどのように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 0x7f4307d642e8>

2次元平面上で赤と青がそれなりに分離できているようなので、ここでPIで計算したベクトルは2つのグループの間の違いをうまく表現できているようです。

ではこのX軸とY軸は何を表現しているのでしょうか？ここでは普通のPCAを使っているので、元の次元の所に ベクトルがあり、各点がそのベクトルの成分をどれだけ含んでいるかを表しています。数式で書くと、 X軸、Y軸に対応するベクトル$v_1, v_2$があり、PIによる高次元ベクトル$w_i$($i$はデータごとのインデックスです)は $$ w_i \approx \mu_{i1} v_1 + \mu_{i2} v_2 + v_0 $$ ただし$v_0$は$w_1,\ldots,$の平均、と近似されます(このあたりはPCAの理論を勉強してください)。 すると、$v_1, v_2$がどのようなベクトルであるかがわかればこの図の意味がわかるわけです。

sklearn の PCA では、$v_0$ が PCA.mean_, $v_1, v_2$ が PCA.components_ に保存されています。

pca.mean_

array([3.32287297e-05, 3.52588389e-05, 3.67793759e-05, ...,
       9.95550058e-08, 7.27036741e-08, 7.13024938e-08])

pca.components_

array([[ 2.59806997e-07,  2.70363482e-07,  2.77565448e-07, ...,
        -3.62878022e-09, -4.09395688e-09, -3.50500694e-09],
       [-6.08939378e-07, -6.42464500e-07, -6.65625771e-07, ...,
         4.65980431e-09, -5.76173853e-10,  1.41943997e-09]])

PIの一つの利点として、ベクトルをパーシステント図に変換できる点があります。これを使うとこれら平均や$v_1, v_2$をパーシステント図の形で 可視化できます。実際にやってみましょう。ベクトル化で使ったhc.PIVectorizerMeshのオブジェクトがこの変換(逆変換)のための情報を持っています。

vectorize_spec.histogram_from_vector(pca.mean_).plot()

vectorize_spec.histogram_from_vector(pca.components_[0, :]).plot(colorbar={"type":"linear-midpoint", "midpoint": 0})
vectorize_spec.histogram_from_vector(pca.components_[1, :]).plot(colorbar={"type":"linear-midpoint", "midpoint": 0})

PCAのベクトルの成分には±があり、それが各軸の±と対応しているので±が赤と黒に対応するよう表示させるとわかりやすいです。 それは、colorbar={"type":"linear-midpoint", "midpoint": 0}の部分で調整できます。ここからさらに この赤い部分、青い部分、に含まれるbirth-death pairを取り出して、逆解析で取り出してさらなる解析をする、ということもできます。 この逆解析の部分は次のロジスティック回帰の部分で実習しましょう。

ロジスティック回帰

次はロジスティック回帰です。ロジスティック回帰もPHとの組み合わせが比較的やりやすい手法です。 特に2クラス分類が解釈のしやすさもあってやりやすいです。

まず学習をする前にテスト用データ(test set)と学習用データ(training set)に分けます。

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

学習用データで学習します。sklearnのLogisitcRegressionというクラスを使います。正則化パラメータCはここでは あらかじめ適当なものを選んでいますが、LogisitcRegressionCVというクラスを使って交差検証を使って自動的にパラメータを 決めることもできます。

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

model.fit(pdvects_train, labels_train, )

LogisticRegression(C=0.01, class_weight=None, dual=False, fit_intercept=True,
                   intercept_scaling=1, l1_ratio=None, max_iter=100,
                   multi_class='warn', n_jobs=None, penalty='l2',
                   random_state=None, solver='lbfgs', tol=0.0001, verbose=0,
                   warm_start=False)

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

model.score(pdvects_test, labels_test)

0.82

なかなか良さそうな精度です。では、学習結果を可視化します。

vectorize_spec.histogram_from_vector(model.coef_.ravel()).plot(colorbar={"type": "linear-midpoint", "midpoint": 0.0})

この図の見方ですが、濃い赤の領域に沢山birth-death pairがあるデータはラベル1である可能性が高く、 濃い青の領域に沢山birth-death pairがあるデータはラベル0である可能性が高いことを意味します。

逆解析 (Optimal volume)

ではこの赤/青の領域を抽出しましょう。 値が0.005を基準として、0.005より大、-0.005より小、の部分を切り出します。 これには numpy の不等号演算が役立ちます。

coef = model.coef_.ravel()

red_area = vectorize_spec.mask_from_vector(coef > 0.005)

blue_area = vectorize_spec.mask_from_vector(coef < -0.005)

取り出した領域を可視化します。

blue_area.plot(title="blue", ax=plt.subplot(1, 2, 1))
red_area.plot(title="red", ax=plt.subplot(1, 2, 2))

取り出した領域に含まれているbirth-death pairを0番目のPD(ラベル0)100番目のPD(ラベル1)から取り出します。

pairs_in_blue_0 = blue_area.filter_pairs(pds[0].pairs())
pairs_in_red_0 = red_area.filter_pairs(pds[0].pairs())

pairs_in_blue_100 = blue_area.filter_pairs(pds[100].pairs())
pairs_in_red_100 = red_area.filter_pairs(pds[100].pairs())

個数を表示しましょう。

len(pairs_in_blue_0), len(pairs_in_red_0), len(pairs_in_blue_100), len(pairs_in_red_100)

(3, 2, 0, 7)

確かに0番目PDには青い領域にあるpairが多く、100番目PDには赤い領域にあるpairが多いようです。実際に見てみると以下のような pairのようです。

pairs_in_blue_0

[Pair(0.008484743298521377, 0.015315059365253281),
 Pair(0.0086407914437143, 0.016345383503399036),
 Pair(0.006687816186736091, 0.017153840879253375)]

さて、ここでoptimal volumeを計算して可視化してみます。 詳しくは説明しませんが、

• optimal volumeの理論: https://epubs.siam.org/doi/abs/10.1137/17M1159439
• optimal volume計算のインターフェース https://www.wpi-aimr.tohoku.ac.jp/hiraoka_labo/homcloud/python-api/interface.html#homcloud.interface.Pair.optimal_volume
• 可視化のインターフェース https://www.wpi-aimr.tohoku.ac.jp/hiraoka_labo/homcloud/python-api/paraview_interface.html

などを参考にしてください。これによってどのような空隙が分類に重要な役割を果たしているのか可視化できます。

optimal_volumes_blue_0 = [pair.optimal_volume(cutoff_radius=0.4) for pair in pairs_in_blue_0]
optimal_volumes_red_100 = [pair.optimal_volume(cutoff_radius=0.4) for pair in pairs_in_red_100]

pv.show([
    hc.OptimalVolume.to_paraview_node_for_volumes(
        optimal_volumes_blue_0,
        gui_name="Optimal volume: blue 0"
    ).threshold("isboundary", 1).set_color((0, 0, 1)).set_linewidth(4),
    pv.PointCloud("pc/0000.txt", dim=3, delimiters=" "),
], wait=False)

pv.show([
    hc.OptimalVolume.to_paraview_node_for_volumes(
        optimal_volumes_red_100,
        gui_name="Optimal volume: red 100"
    ).threshold("isboundary", 1).set_color((1, 0, 0)).set_linewidth(4),
    pv.PointCloud("pc/0100.txt", dim=3, delimiters=" "),
], wait=False)

ちょっとわかりにくい気もするので2つ重ねて表示してみます。

pv.show([
    hc.OptimalVolume.to_paraview_node_for_volumes(
        optimal_volumes_blue_0,
        gui_name="Optimal volume: blue 0"
    ).threshold("isboundary", 1).set_color((0, 0, 1)).set_linewidth(4),
    hc.OptimalVolume.to_paraview_node_for_volumes(
        optimal_volumes_red_100,
        gui_name="Optimal volume: red 100"
    ).threshold("isboundary", 1).set_color((1, 0, 0)).set_linewidth(4),
], wait=False)

青い空隙がラベル0に「典型的」な空隙で、赤い空隙がラベル1に「典型的」な空隙です。

青のほうが赤よりちょっと細い感じがすると思います。学習結果のパーシステント図を見ると、 青い領域は赤い領域よりY軸(death)で小さいほうにあります。実は2次のパーシステント図において、 deathはその空隙に収まる最大の球の半径の2乗に対応しています。つまりこの学習結果は ラベル1のほうがこのスケールの隙間が大きい、ことを意味しています。 赤い領域の中心のdeathが0.022、青い領域の中心が0.017くらいですので、 それぞれの平方根、つまり

np.sqrt(0.017), np.sqrt(0.022)

(0.130384048104053, 0.14832396974191325)

あたりがそれぞれのグループの特徴的空隙のスケールである、ということがわかります。 元々のパーシステント図を眺めると、もっと小さい空隙はどちらも共通に沢山ある、ということも 見てとれます。

以上でこのチュートリアルは終わりです。