機械学習チュートリアル（ポイントクラウド編）¶

まずは、必要なライブラリをインポートします。Python の機械学習ライブラリである scikit-learn などを使用します。 scikit-learn の詳細については、公式ドキュメント（https://scikit-learn.org/stable/index.html）を参照してください。 また、その他の使用ライブラリについても各自で学習してください。

In [2]:
import os  # for makedirs
import homcloud.interface as hc  # HomCloud 
import homcloud.plotly_3d as p3d  # 3D visualizaiton
import plotly.graph_objects as go  # 3D visualizaiton
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
import plotly.io
plotly.io.renderers.default = 'notebook'  # plotly Setting. Choose your favorite renderer.

データの確認¶

機械学習を始める前に、まずは使用するデータの内容を確認しましょう。 このチュートリアルで使用するデータは、pc ディレクトリ内に格納されています。

In [3]:
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 つのグループに分類されていることを示しています。

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

In [4]:
labels = np.loadtxt("pc/label.txt")
labels
Out[4]:
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.])

これは、0 または 1 の値を持つ 200 個のラベルからなる配列です。 それでは、どのデータが 0 または 1 にラベル付けされているのかを確認してみましょう。 ここでは NumPy の機能を使用します。

In [5]:
np.nonzero(labels == 0), np.nonzero(labels == 1)
Out[5]:
((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 番目のデータを可視化してみましょう。

In [6]:
go.Figure(data=[p3d.PointCloud(np.loadtxt("pc/0000.txt"))], layout=dict(scene=p3d.SimpleScene()))
In [7]:
go.Figure(data=[p3d.PointCloud(np.loadtxt("pc/0100.txt"))], layout=dict(scene=p3d.SimpleScene()))

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

それでは、ファイルからポイントクラウドデータを読み込んでみましょう。

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

100 個のデータはすべて 1000×3 の配列です。 これは、3 次元空間上に存在する 1000 個の点からなるポイントクラウドが 100 個あることを意味します。 以下のコードでその事実を確認します。

In [9]:
all(pointcloud.shape == (1000, 3) for pointcloud in pointclouds)
Out[9]:
True

パーシステンス図を用いた機械学習を行う場合、各ポイントクラウドの点の個数が完全に同じでなくても構いませんが、ある程度揃っていないと、うまく機能しないことがあります。

パーシステンス図の計算¶

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

保存されたデータは再利用可能なので、別の解析を行う際に毎回再計算する必要はありません。

In [10]:
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 次のパーシステンス図を使用します。つまり、空隙構造（ボイド）に注目して解析を行います。 もちろん、0 次や 1 次のパーシステンス図を使うことも可能ですが、このデータでは 2 次の特徴が最も顕著に現れるため、今回はそれに焦点を当てます。

今回は詳しくは説明しませんが、より高度な解析手法として、1 次と 2 次の情報を同時に活用する方法もあります。

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

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

In [12]:
pds[0].histogram(x_range=(0, 0.03)).plot(colorbar={"type": "log"})
pds[100].histogram(x_range=(0, 0.03)).plot(colorbar={"type": "log"})
No description has been provided for this image
No description has been provided for this image

Persistence Image（PI）によるベクトル化¶

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

  • http://www.jmlr.org/papers/volume18/16-337/16-337.pdf

ただし、ここで使用する手法はこの論文のオリジナルな方法ではなく、以下の文献で提案されている改変版に基づいています：

  • https://link.springer.com/article/10.1007%2Fs41468-018-0013-5

まずは、ベクトル化の設定を定義します。
(0, 0.03), 128 という指定は、[0, 0.03] × [0, 0.03] の正方形領域を 128×128 のグリッドに分割して計算することを意味します。
sigma=0.002 はガウス分布の標準偏差を表します。
また、weight=("atan", 0.01, 3) により、対角線からの距離に応じた重み付けが行われます。
この場合、重みは $\textrm{atan}(0.01 \cdot \ell^3)$ によって決まり、ここで $\ell$ は lifetime（= death − birth）を意味します。

In [13]:
vectorize_spec = hc.PIVectorizeSpec((0, 0.03), 128, sigma=0.002, weight=("atan", 0.01, 3))

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

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

ベクトルの要素の最大値，最小値を確認します。

In [15]:
pdvects.min(), pdvects.max()
Out[15]:
(np.float64(0.0), np.float64(2.099302720413663e-10))

最大値はおよそ $2 \times 10^{-10}$ と非常に小さい値です。
このままでは機械学習に悪影響を与える可能性があるため、正規化を行うことが望ましいでしょう。

正規化にはさまざまな方法がありますが、ここでは最も単純な方法として、配列の最大値が 1 になるようにスケーリングします。

In [16]:
pdvects = pdvects / pdvects.max()

主成分分析（PCA）¶

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

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

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

In [17]:
pca = PCA(n_components=2)
pca.fit(pdvects)
Out[17]:
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
n_components  2
copy  True
whiten  False
svd_solver  'auto'
tol  0.0
iterated_power  'auto'
n_oversamples  10
power_iteration_normalizer  'auto'
random_state  None

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

In [18]:
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のデータを青で描画
Out[18]:
<matplotlib.collections.PathCollection at 0x7276ecae3f50>
No description has been provided for this image

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

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

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, 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_ に格納されています。

In [19]:
pca.mean_
Out[19]:
array([3.32287297e-05, 3.52588389e-05, 3.67793759e-05, ...,
       9.95550058e-08, 7.27036741e-08, 7.13024938e-08], shape=(8256,))
In [20]:
pca.components_
Out[20]:
array([[ 2.59806997e-07,  2.70363482e-07,  2.77565448e-07, ...,
        -3.62878022e-09, -4.09395688e-09, -3.50500694e-09],
       [-6.08939377e-07, -6.42464501e-07, -6.65625771e-07, ...,
         4.65980431e-09, -5.76173853e-10,  1.41943997e-09]],
      shape=(2, 8256))

Persistence Imageの利点の一つは、ベクトルをパーシステンス図と同様の符号付きヒストグラムに変換できる点です。
この機能を使えば、平均ベクトルや主成分ベクトル（$v_1$, $v_2$）をパーシステンス図の形式で可視化することができます。

実際にやってみましょう。ベクトル化に使用した hc.PIVectorizerMesh オブジェクトは、この変換に必要な情報を保持しており、histogram_from_vector メソッドを使ってヒストグラムに変換できます。

In [21]:
vectorize_spec.histogram_from_vector(pca.mean_).plot()
No description has been provided for this image
In [22]:
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})
No description has been provided for this image
No description has been provided for this image

PCA のベクトルには正負の成分があり、それが各軸の正負方向に対応しています。
そのため、ヒートマップなどで表示する際には、正の値を赤、負の値を黒などに対応させると視覚的にわかりやすくなります。
これは colorbar={"type": "linear-midpoint", "midpoint": 0} の設定によって調整できます。

さらに、赤や青の領域に対応する birth-death ペアを抽出し、逆解析によってそれらの特徴を詳しく調べることも可能です。
この逆解析の方法については、次のロジスティック回帰のセクションで実習してみましょう。

ロジスティック回帰¶

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

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

In [23]:
pdvects_train, pdvects_test, labels_train, labels_test = train_test_split(pdvects, labels, test_size=0.25)

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

In [24]:
model = lm.LogisticRegression(C=0.01, solver="lbfgs")
In [25]:
model.fit(pdvects_train, labels_train, )
Out[25]:
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
penalty  'l2'
dual  False
tol  0.0001
C  0.01
fit_intercept  True
intercept_scaling  1
class_weight  None
random_state  None
solver  'lbfgs'
max_iter  100
multi_class  'deprecated'
verbose  0
warm_start  False
n_jobs  None
l1_ratio  None

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

In [26]:
model.score(pdvects_test, labels_test)
Out[26]:
0.92

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

In [27]:
vectorize_spec.histogram_from_vector(model.coef_.ravel()).plot(colorbar={"type": "linear-midpoint", "midpoint": 0.0})
No description has been provided for this image

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

逆解析¶

そこで、この赤と青の領域について詳しく調べてみましょう。まず、パーシステンス図から、それぞれの領域に含まれるペアを抽出します。

そのために、赤と青の領域を切り出します。基準値として0.005を用い、0.005より大きい部分と-0.005より小さい部分を抽出します。NumPyの不等号演算を使って真偽値ベクトルを作成し、それを vectorizer_spec の mask_from_vector メソッドに渡すことで、MaskHistogram オブジェクトを構築します。

このようにして構築されたオブジェクトは、切り出した領域の可視化や、領域内に含まれるペアの抽出に利用できます。

In [28]:
coef = model.coef_.ravel()
In [29]:
red_area = vectorize_spec.mask_from_vector(coef > 0.005)
In [30]:
blue_area = vectorize_spec.mask_from_vector(coef < -0.005)

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

In [31]:
blue_area.plot(title="blue", ax=plt.subplot(1, 2, 1))
red_area.plot(title="red", ax=plt.subplot(1, 2, 2))
No description has been provided for this image

切り出した領域に含まれる birth-death ペアを MaskHistogram の filter_pair メソッドを用いて0番目のパーシステンス図（ラベル0）と100番目のパーシステンス図（ラベル1）から抽出します。

In [32]:
pairs_in_blue_0 = blue_area.filter_pairs(pds[0].pairs())
pairs_in_red_0 = red_area.filter_pairs(pds[0].pairs())
In [33]:
pairs_in_blue_100 = blue_area.filter_pairs(pds[100].pairs())
pairs_in_red_100 = red_area.filter_pairs(pds[100].pairs())

個数を表示しましょう。

In [34]:
len(pairs_in_blue_0), len(pairs_in_red_0), len(pairs_in_blue_100), len(pairs_in_red_100)
Out[34]:
(3, 2, 1, 5)

確かに、0番目のパーシステンス図には青い領域に含まれるペアが多く、100番目のパーシステンス図には赤い領域に含まれるペアが多いようです。実際に確認してみると、以下のようなペアが含まれていることがわかります。

In [35]:
pairs_in_blue_0
Out[35]:
[Pair(0.008484743298521374, 0.015315059365253271),
 Pair(0.0086407914437143, 0.01634538350339904),
 Pair(0.006687816186736094, 0.01715384087925338)]

ここで、stable volume を計算し、可視化を行います。stable volume は、以前のポイントクラウドのチュートリアルでも使用した道具です。これにより、どのような空隙構造が分類において重要な役割を果たしているのかを視覚的に確認することができます。

In [36]:
stable_volumes_blue_0 = [pair.stable_volume(0.0003, cutoff_radius=0.4) for pair in pairs_in_blue_0]
stable_volumes_red_100 = [pair.stable_volume(0.0003, cutoff_radius=0.4) for pair in pairs_in_red_100]
In [37]:
go.Figure(data=[
    ov.to_plotly3d_mesh(color="blue") for ov in stable_volumes_blue_0
] + [
    p3d.PointCloud(pointclouds[0])
], layout=dict(scene=p3d.SimpleScene()))
In [38]:
go.Figure(data=[
    ov.to_plotly3d_mesh(color="red") for ov in stable_volumes_red_100
] + [
    p3d.PointCloud(pointclouds[100])
], layout=dict(scene=p3d.SimpleScene()))

違いがわかりにくいので重ねて描画してみましょう。

In [39]:
fig = go.Figure(data=[
    ov.to_plotly3d_mesh(color="blue") for ov in stable_volumes_blue_0
] + [
    ov.to_plotly3d_mesh(color="red") for ov in stable_volumes_red_100
], layout=dict(scene=p3d.SimpleScene()))
fig.update_traces(opacity=0.5, selector=dict(type="mesh3d"))

青い体積はラベル0に「典型的」な空隙を、赤い体積はラベル1に「典型的」な空隙を表しています。

青のほうが赤よりもやや細く、小さい印象を受けると思います。学習結果のパーシステンス図を見ると、青い領域は赤い領域よりもY軸（Death）の値が小さい位置にあります。実は、2次のパーシステンス図において death の値は、その空隙に収まる最大の球の半径の二乗に対応しています。つまり、この学習結果は「ラベル1の方が、このスケールで見たときに大きな空隙を持っている」ことを意味しています。

具体的には、赤い領域の中心の death 値は約 0.022、青い領域の中心は約 0.017 です。これらの平方根を取ると、それぞれの空隙に収まる最大の球の半径がわかります。つまり：

In [40]:
np.sqrt(0.017), np.sqrt(0.022)
Out[40]:
(np.float64(0.130384048104053), np.float64(0.14832396974191325))

あたりが、それぞれのグループにとって特徴的な空隙のスケールであることがわかります。

また、元のパーシステンス図を眺めると、より小さな空隙は両方のグループに共通して多数存在していることも確認できます。

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