はじめに
人工知能（AI）の分野の中で深層学習（Deep Learning）の技術が進化、そして画期的な生成AI登場により、世の中猫も杓子もAI祭りである。社会人である私も、上司からなんでもいいからAIを使うのじゃと命令されておるのであります。でもひねくれものの私は、どうにも深層学習のブラックボックス的な感じがもやもやするわけです。
下記に示すMNISTデータは、機械学習でよく使われる手書き数字画像のデータセットです。すべての数字画像のサイズは28x28ピクセルで、本当に手書きっぽい（少し読みにくい）データが並んでいます。
By Suvanjanprasai - Own work, CC BY-SA 4.0, Link
人間っぽい分類はできないだろうか？
数字ひとつに対して、MNISTの28x28ピクセル合計784個のDataをわたしたちは、例えば右から3番目、上から10番目のピクセルが白か黒かはという情報によって数字を認識していませんよね（と思っているだけかもしれませんが）。数字の形が丸いとか、穴が空いているとか、線になっているとか、形を把握した情報をもとに、数字という文字を判断していると思っています。なので読み間違えは当然あるわけです。上記のサンプルの数字は、1と7は区別がつかないものも多くありますよね。
深層学習のような複雑なピクセルごと相互関係、重み付けにより特徴値が自動抽出され、解釈が難しいとされるブラックボックスとは異なる、もっと直感的な数字画像の形の情報をもとに、自らの思いついた特徴値を使ってMNISTのData分析ができないかと思ったのが、この記事を書くきっかけでです。
トポロジー
ここで登場するのが、「位相幾何学」とも呼ばれるトポロジーです。この学問は幾何学の分野の1つであり、図形を構成する点の連続的位置関係のみに着目してその性質を研究するものです。
トポロジーはいろいろわかりやすい書物に書かれていますし、書物にもかかれている様に、何らかの形を連続変化、何らかの形を連続変形（伸ばしたり曲げたりすることはするが切ったり貼ったりはしないこと）しても保たれる性質（位相不変量と言う）に焦点を当てたものです。よく引き合いに出されるのがコーヒーカップとドーナツですね。これらは連続に変形すればトポロジーとしては位相不変量が同一というわけです（下記のWikiの図参照）。
Lucas Vieira - https://commons.wikimedia.org/w/index.php?curid=1236079
トポロジーの考えを数字に当てはめてみましょう。8は穴が2つ、0と6と9は穴が1つの画像を期待したいですが、実際のMNISTのDataをみると、手書きの影響により、期待通りの穴の数となっていない場合があります。これはしょうがないか 
ホモロジーからパーシステントホモロジーへ
位相不変量にはいくつか種類がありますが、よく使われるのがホモロジーです。詳細は専門書に任せるとして、
- 0次のホモロジーは、構成要素の連結部分の数
- 1次のホモロジーは、構成要素の中にある穴の数
- 2次のホモロジーは、3次元にて構成要素にある空洞
（トムとジェリーに出てくるチーズみたいな感じ）
MNISTのDataに対してこれら0次と1次のホモロジーを求めることで、数字を判別することができるのではないかと思ったわけです。数字の8は1次のホモロジーが2であることを期待しますが、しかしながら、手書きでは下図のよう右上が切れている場合もあります。これだと1次のホモロジーは、1となってしまいます。
ここで登場するのが、パーシステントホモロジーです。この8の字のひとつひとつのピクセルを太らせるとどうでしょうか？8の数字の右上が切れているのが塞がって穴をつくり、さらにピクセルを太らせると最終的は穴が塞がってしまいます。パーシステントホモロジーでは、ピクセルを痩せさせたり、太らせたりして、ホモロジーの数だけでなくホモロジーの生死（1次のホモロジーの場合、穴が生成されてなくなるまで）がどうなるかという情報を入手することができるのです。
HomCloud
世の中の先人たちのお陰で、このパーシステントホモロジーを扱うライブラリがすでに用意されています。下記のサイトを参照すれば、簡単にpipを使ってインストールすることができます。
HomCloud
https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2166-23.pdf
ではさっそくこのツールを入手し、トポロジーをつかったData分析をしていきましょう。期待したいのは、8の正答率は高くあってほしいな。
まずはそのまま分析開始
Dataの準備
MMINSTのDataをシンプルにするため、グレースケールを白黒のみ変更し、閾値は階調255のうち、100にしておきました。より膨らんだ状態がいいか、痩せたほうがいいかは本来はチューニングの要素かと思います。Data数は10000にて。
from tensorflow.keras.datasets import mnist
from sklearn.model_selection import train_test_split
# データ読み込み
(X_full, y_full), (X_test_full, y_test_full) = mnist.load_data()
# Data関連 Data数の定義と、2値化（白黒はっきりさせる、グレーな部分は白に寄せる）
data_num = 10000
X_full = X_full[:data_num] > 100
y_full = y_full[:data_num]
X_train, X_val, y_train, y_val = train_test_split(X_full, y_full, test_size=0.2, random_state=42)
実装
計算結果を比較するために、ピクセルベースの学習も同時に計算するプログラムを作成しました。それぞれ標準化を実施しています。分類はサポートベクトルマシン分類(SVC）を使いました。
import numpy as np
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
# ピクセルの場合
X_train_pixels = X_train.reshape(len(X_train), -1).astype(np.float32) / 255.0 # 784次元 標準化
X_val_pixels = X_val.reshape(len(X_val), -1).astype(np.float32) / 255.0
## 標準化
scaler_pixels = StandardScaler()
X_train_pixels = scaler_pixels.fit_transform(X_train_pixels)
X_val_pixels = scaler_pixels.transform(X_val_pixels)
## 学習・評価
clf_pixels = SVC(kernel='rbf', random_state=42, C=10.0, class_weight='balanced')
clf_pixels.fit(X_train_pixels, y_train)
score_pixels = clf_pixels.score(X_val_pixels, y_val)
# パーシステントホモロジーの場合
X_train_ph = np.array([compute_ph_features(x) for x in X_train])
X_val_ph = np.array([compute_ph_features(x) for x in X_val])
## 標準化
scaler_ph = StandardScaler()
X_train_ph = scaler_ph.fit_transform(X_train_ph)
X_val_ph = scaler_ph.transform(X_val_ph)
## 学習・評価
clf_ph = SVC(kernel='rbf', random_state=42, C=10.0, class_weight='balanced')
clf_pixels.fit(X_train_pixels, y_train)
score_pixels = clf_pixels.score(X_val_pixels, y_val)
print("Validation accuracy with raw pixel features SVC:", score_pixels)
print("Validation accuracy with PH summary features SVC:", score_ph)
Homcloudをつかってパーシステントホモロジーを求め特徴値を定義
さきほどのコードにあるcompute_ph_features(x)にて（詳細コードは）、パーシステントホモロジーを求めます。ここからわかるのは、0次と1次のホモロジーそれぞれの発生時間と消滅時間、総数であり、発生時間から消滅時間を差し引いた寿命を求めることができます。例えば、0次のホモロジーの場合は、構成要素の連結部分の数ですので、シンプルな形に比べて、濃淡があってくねくねしている形は数が多く、1次のホモロジーの場合は、穴の存在なので、大きな穴であればあるほど寿命が長いことになります。この関数により、なんとなく思いつきそうな特徴値を下記のように定義します。ここはもっといい方法がきっとあるかと思います。今後の課題ですね。
- 特徴値 feat0 は、0次のホモロジーの数、寿命の平均、寿命の最大値、寿命の最小値
- 特徴値 feat1 は、1次のホモロジーの数、寿命の平均、寿命の最大値、寿命の最小値
つまり784のピクセルのDataをパーシステントホモロジーによって合計8つにして計算したことになります。少しわかりにくい0次のホモロジーの特性を理解するには、この記事が非常に有益でした。
そのままの計算結果
下記のとおりとなりました。
Validation accuracy with ピクセル SVC: 0.932
Validation accuracy with パーシステントホモロジー SVC: 0.437
ピクセル分析に、ダブルスコアをつけられてしまいました。完敗です。個別の数字を見て分析すると下記のマトリクスにあるように、正答率が50%を切っているのは、2,3,4,5,6です。トポロジー的に穴があるなしで言えば、4という数字は悩ましく、MNISTの生データをみても、4の数字の上の部分をくっつける人もいれば離して書く人も多く、特徴値の効果が低いと思われます。また、6に関して言えば、9と予測されている結果になっています。
トポロジー的特徴の検証
そこで、トポロジー的特徴が分類できているかどうか、穴が1つ（数字が0,6,9）、穴がゼロ(1,2,3,4,5,7)、穴が2つ(8)の3つに分けた場合を見てみましょう。4は穴なしとしました。
この3つの分類だけのことを言えば、86.7%の正答率が確保されていることがわかりました。そんなにわるくないのです。では、これら特徴値で不足している部分はいったいなんなのかを考えていきたいと思います。
=== Class-wise Accuracy (3-class) ===
Group A (0,6,9): accuracy = 0.8304
Group B (1,2,3,4,5,7): accuracy = 0.9012
Group C (8): accuracy = 0.7530
Overall accuracy (PH Features, 3-class): 0.8670
チューニング
パーシステントホモロジーの精度を高める
Homcloudでは、bitmap単位に計算をします。よって、穴があっても小さい場合は、寿命が短く、すぐに塞がってしまうことが考えられます。よって、この画像Dataを拡大すると結果が変わってくる可能性があります。画像を2倍、5倍にした場合の結果を見てみましょう。闇雲に細かくすればいいというわけでなく、2倍にした場合、正答率が50%を超えることができました。
- Scale = 2 (56x56) の場合、Validation accuracy SVC: 0.5045
- Scale = 5 (140x140) の場合、Validation accuracy SVC: 0.481
パーシステントホモロジーの苦手？を目印で克服する
突然ですが、QRコードはご存知かと思いますが、このコードには、四角の頂点のうち3つにマークが付いています。位置検出用パターンというのものです。今回の分類に使ったトポロジーの特徴はなんといっても連続的位置関係のみというところですから、そもそも上下左右反転や鏡像に対しても全く関係ないはずです（だったら分類に使うなと言わずに）。なので、あえて画像に位置検出用パターンみたいなものを埋め込んでみます。ただ単純に、左上に、白い空白をつけます。正方形の場合（パターン1）と長方形（パターン2）を試してみました。このマークは、0次のホモロージーのゼロ値を避ける効果（さきほどの記事参照）もありました。
| パターン1 | パターン２ |
|---|---|
| mark_height = 5, mark_width = 5 | mark_height = 15, mark_width = 3 |
 |
 |
| accuracy: 0.582 | accuracy: 0.6085 |
これにより、正答率60%達成です！マトリクスは下記のように変化しました。個別の数字の正答率を調べたところ、2,4,5以外の数字は50%以上の正答率です。
直感的な情報を付加する
正答率が低い2と5は、眺めてみると本当によく似ています。それぞれ丸い曲線と横のまっすぐな部分とか。ここは、ちょっとずるいですが、トポロジー部分以外において特徴値を加えたいと思います。特徴値を極力シンプルにするため、数字の位置情報の特徴をとりあげます。
- 特徴値 feat2 は、28x28ピクセル内で、右半分左半分の白ドットの数、上半分下半分の白ドットの数（合計4つの特徴値）
- 特徴値 feat3 は、28x28ピクセル内で、右斜め上半分、右斜め下半分の白ドットの数（合計2つの特徴値）
イメージとしては、縦横斜めのバランスは、各数字によって違いが現れるだろうという仮設のものです。本来画像全体の重心を合わせてからやるべきかと思うところですが、ここは時間切れ。
いままでの特徴値を含め、4種類（合計14個）を使って、計算すると
accuracy : 0.8135
となりました。下記に占めるマトリクスも全数字で50％を超えることができました。
まとめ
今回は、MNISTのDataの文字の形にこだわり、トポロジーの考えを利用して特徴値を取り出し、数字の分類を行いました。穴がある数字の形に関してはよく分類ができていました。さらに位置情報なや左右の数字の配置の情報を取り混ぜ、4種類（合計14個）のDataにより、最終な分類正答率は、80%を超えることができました。
参考までに、4種類の特徴値の効果を比較検証するために、画像ファイルのピクセルを2倍にし、左上にマークをつけた状態のものに対してそれぞれの正答率を計測したのが下記の表です。
| 特徴 | accuracy | |
|---|---|---|
| feat0 | 0次のホモロジー（複雑度合い） | 0.373 |
| feat1 | 1次のホモロジー（穴の大きさや数） | 0.4625 |
| feat0+1 | 0.6085 | |
| feat2 | 上下左右のバランス | 0.453 |
| feat3 | 右斜上下のバランス | 0.3625 |
| 全部入り | 0.8135 |
これらの分類の状態を可視化するために、PCAによる2要素分解をして図示しています。なんとなく、1や0はうまく分類できている様子が伺えます。
最後にピクセルの784個のDataに今回編み出した14個のDataを追加して分類結果がどうなるか？をやってみたのですが、0.05%しか貢献することができませんでした
Pixel shape: (10000, 784)
PH shape: (10000, 14)
Combined shape: (10000, 798)
Validation accuracy with pixels : 0.932
Validation accuracy with COMBINED (pixels + PH): 0.937
今回は、Homcloudの使い方を少し理解することができました。本来は物性物理学で強力なツールなのに、畑違いなところに応用してしまったので申し訳ない感じですが、普段目にしている数字がとても奥深く感じられるようになったかな。ありがとうございました。あと投稿おそくなってすみません～
参照
使用した主なコード
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.svm import SVC
from tensorflow.keras.datasets import mnist
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
# Homcloud
import homcloud.interface as hc
# ========== データ読み込み
(X_full, y_full), (X_test_full, y_test_full) = mnist.load_data()
# Data関連 Data数の定義と、2値化（白黒はっきりさせる、グレーな部分は白に寄せる）
data_num = 10000
X_full = X_full[:data_num] > 100
y_full = y_full[:data_num]
X_train, X_val, y_train, y_val = train_test_split(X_full, y_full, test_size=0.2, random_state=42)
# ========== 変換と特徴値抽出
# ========== 単純拡大変換
def upscale_mnist_with_kron(img28, scale):
ones_block = np.ones((scale, scale), dtype=img28.dtype)
return np.kron(img28, ones_block)
# ========== マーク埋込
def embed_white_mark0(img, mark_height, mark_width):
out_img = img.copy()
# 左上 corner 領域を 255 (白) に置き換え
out_img[0:mark_height, 0:mark_width] = 255
return out_img
# ========== 画像の上下左右Dataを取得
def direction_features(img28):
# img28: shape(28, 28)
# 左右分割
left = img28[:, :14] # 左半分
right = img28[:, 14:] # 右半分
# 上下分割
top = img28[:14, :]
bottom = img28[14:, :]
# 各領域で 画素の合計 or 平均 などを計算
f_left_sum = left.sum()
f_right_sum = right.sum()
f_top_sum = top.sum()
f_bottom_sum = bottom.sum()
return np.array([f_left_sum, f_right_sum,f_top_sum,f_bottom_sum])
# ========== 画像の対角線を境に(左斜め半分, 右斜め半分)の2領域に分割し画素取得
def diagonal_split_features(img28):
rows = np.arange(28).reshape(-1, 1) # shape (28,1)
cols = np.arange(28).reshape(1, -1) # shape (1,28)
row_indices = rows.repeat(28, axis=1) # shape (28,28)
col_indices = cols.repeat(28, axis=0) # shape (28,28)
# 左斜半分: r>c
left_diag_mask = (row_indices > col_indices)
# 右斜半分: r<c
right_diag_mask = (row_indices < col_indices)
left_diagonal_sum = img28[left_diag_mask].sum()
right_diagonal_sum = img28[right_diag_mask].sum()
return np.array([left_diagonal_sum,right_diagonal_sum])
# ========== ヘルパー関数群 ==========
def get_persistence_diagrams(image_array):
"""
MNIST画像からHomCloudのPython APIを用いてパーシステンスホモロジーを計算。
bitmap_levelset_persistenceを利用
"""
# 左上マーク対応
image_array = embed_white_mark0(image_array.reshape(28, 28), mark_height, mark_width)
# 拡大対応
image_array = upscale_mnist_with_kron(image_array, scale)
pdlist = hc.PDList.from_bitmap_levelset(hc.distance_transform(image_array, signed=True))
#0次ホモロジー
pd0 = np.vstack([pdlist.dth_diagram(0).births, pdlist.dth_diagram(0).deaths]).transpose()
#1次ホモロジー
pd1 = np.vstack([pdlist.dth_diagram(1).births, pdlist.dth_diagram(1).deaths]).transpose()
return pd0, pd1
def summary_features0(data):
"""
0次ホモロジー：data: (birth, death) のペアからなる numpy配列
統計量(個数、平均persistence、最小persistence、最大persistence)を返す
"""
if data.size == 0:
return np.array([100.0, 1.0, 0.0, 100.0])
pers = data[:,1] - data[:,0]
return np.array([len(pers), np.mean(pers), np.min(pers), np.max(pers)])
def summary_features1(data):
"""
1次ホモロジー：data: (birth, death) のペアからなる numpy配列
統計量(個数、平均persistence、最小persistence、最大persistence)を返す
"""
if data.size == 0:
return np.array([0.0, 0.0, 0.0,0.0])
pers = data[:,1] - data[:,0]
#return np.array([len(pers), np.mean(pers), np.max(pers)])
return np.array([len(pers), np.mean(pers), np.min(pers), np.max(pers)])
def compute_ph_features(image_array):
"""
特徴量を取得
"""
pd0, pd1 = get_persistence_diagrams(image_array)
feat0 = summary_features0(pd0)
feat1 = summary_features1(pd1)
feat2 = direction_features(image_array)
feat3 = diagonal_split_features(image_array)
return np.concatenate([feat0, feat1, feat2, feat3])
# ========== MAIN
# 初期パラメーター関連
data_num = 10000
scale = 2
mark_height = 15
mark_width = 3
# ピクセルの場合
X_train_pixels = X_train.reshape(len(X_train), -1).astype(np.float32) / 255.0 # 784次元 標準化
X_val_pixels = X_val.reshape(len(X_val), -1).astype(np.float32) / 255.0
## 標準化
scaler_pixels = StandardScaler()
X_train_pixels = scaler_pixels.fit_transform(X_train_pixels)
X_val_pixels = scaler_pixels.transform(X_val_pixels)
## 学習・評価
clf_pixels = SVC(kernel='rbf', random_state=42, C=10.0, class_weight='balanced')
clf_pixels.fit(X_train_pixels, y_train)
score_pixels = clf_pixels.score(X_val_pixels, y_val)
# パーシステントホモロジーの場合
X_train_ph = np.array([compute_ph_features(x) for x in X_train])
X_val_ph = np.array([compute_ph_features(x) for x in X_val])
## 標準化
scaler_ph = StandardScaler()
X_train_ph = scaler_ph.fit_transform(X_train_ph)
X_val_ph = scaler_ph.transform(X_val_ph)
## 学習・評価
clf_ph = SVC(kernel='rbf', random_state=42, C=10.0, class_weight='balanced')
clf_ph.fit(X_train_ph, y_train)
score_ph = clf_ph.score(X_val_pixels, y_val)
print("Validation accuracy with raw pixel features SVC:", score_pixels)
print("Validation accuracy with PH summary features SVC:", score_ph)
# ========== 特徴空間のPCAによる可視化
pca_ph = PCA(n_components=2,svd_solver='auto')
proj_ph = pca_ph.fit_transform(X_train_ph)
plt.figure(figsize=(8,6))
for label in np.unique(y_train):
mask = (y_train == label)
plt.scatter(proj_ph[mask,0], proj_ph[mask,1], label=label, alpha=0.7)
plt.legend()
plt.title("PH summary features (Train set) - PCA projection")
plt.xlabel("PC1")
plt.ylabel("PC2")
plt.show()
参考