Vietoris-Rips フィルトレーションのチュートリアル¶

このチュートリアルでは、Vietoris-Rips フィルトレーションを用いてパーシステンス図を計算する方法を紹介します。

Vietoris-Rips フィルトレーションの利点¶

  • すべての点のペアの間の距離情報（距離行列）があれば計算可能
    • ユークリッド空間に限らず、遺伝子間の距離のような非ユークリッド空間でも適用可能です
  • アルファフィルトレーションは現在、2次元または3次元に限定されますが、Vietoris-Rips フィルトレーションはより高次元のデータにも対応できます

欠点と注意点¶

  • 計算コストが高い（組み合わせ爆発

    • 特に高次のホモロジーの計算は困難です
    • 実用的な計算のためには、計算する次数を制限する必要があります
    • また、ある距離以上の辺を無視することで、計算時間を抑える工夫が必要です

  • 幾何的な情報が近似的にしか得られない

    • 各点に球を置くような手法と比べて、幾何的な構造を正確に反映していない場合があります
    • このため、得られる結果の解釈に注意する必要があります

準備¶

必要なライブラリを読み込みます。

In [2]:
import numpy as np
import homcloud.interface as hc
from scipy.spatial import distance_matrix
import matplotlib.pyplot as plt

200次元のデータを100個、ランダムに準備します。

In [3]:
data = np.random.uniform(size=(100, 200))

距離行列 を scipy.spatial.distance_matrix を使って計算します。

In [4]:
dmatrix = distance_matrix(data, data)

distance matrix は 100x100 の配列です。 各要素は0以上の値を持ちます。

In [5]:
dmatrix.shape
Out[5]:
(100, 100)

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

homcloud.interface.PDList.from_rips_filtration を使って、距離行列から Vietoris-Rips フィルトレーションを経由してパーシステンス図を計算できます。主な引数の意味は以下の通りです：

  • 最初の引数：距離行列（2点間の距離情報を格納した対称行列）
  • maxdim：計算するホモロジーの最大次数（例：1なら 0 次と 1 次のホモロジーを計算）
  • save_to：計算結果（パーシステンス図）の保存先ファイル名
  • maxvalue（省略可能）：距離の最大値を指定（これを超える辺は無視されます）
In [6]:
%%time
hc.PDList.from_rips_filtration(dmatrix, maxdim=3, save_to="rips.pdgm")

CPU times: user 1.59 s, sys: 71.1 ms, total: 1.66 s
Wall time: 1.65 s
Out[6]:
PDList(path=rips.pdgm)

計算時間と maxdim の関係¶

maxdim の値によって、パーシステンス図の計算時間は大きく変化します：

  • maxdim = 2：ほぼ一瞬で完了
  • maxdim = 3：5〜10秒程度
  • maxdim = 4：約2分
  • maxdim > 4：未検証ですが、計算量はmaxdimに対して指数関数的に増加すると考えられます

このように、次数が高くなるほど計算コストが急激に増加するため、実用上は maxdim を適切に制限することが重要です。

これらの計算は、Ripserという高速なパーシステントホモロジー計算ライブラリによって実行されています。Ripser は、Vietoris-Rips フィルトレーションに基づくパーシステンス図の計算において、現時点で最速クラスの実装とされています。

パーシステンス図の可視化¶

以下のセルでパーシステンス図を読み込みます。

In [7]:
pdlist = hc.PDList("rips.pdgm")

1次のパーシステンス図を可視化します。

In [8]:
pd1 = pdlist.dth_diagram(1)
pd1.histogram().plot(colorbar={"type": "log"})
No description has been provided for this image

2次，3次のパーシステンス図も可視化します。

In [9]:
pd2 = pdlist.dth_diagram(2)
pd2.histogram().plot(colorbar={"type": "log"})
No description has been provided for this image
In [10]:
pd3 = pdlist.dth_diagram(3)
pd3.histogram().plot(colorbar={"type": "log"})
No description has been provided for this image

maxvalue の効果¶

ここでは、パーシステンス図を計算する際に maxvalue を指定すると何が起こるのかを見ていきましょう。

まず、距離行列の値がどのように分布しているかを確認するために、ヒストグラムを描いてみます。

重複を避けるため距離行列の上半分を取り出すには、numpy.triu_indices を使うと便利です。

In [11]:
plt.hist(dmatrix[np.triu_indices(100, 1)], 100); None
No description has been provided for this image

5〜6.5の間に分布しています。

そこで maxvalue=5.6 として計算しましょう。

In [12]:
%%time
hc.PDList.from_rips_filtration(dmatrix, maxdim=3, maxvalue=5.6, save_to="rips-5.6.pdgm")

CPU times: user 34.8 ms, sys: 942 μs, total: 35.7 ms
Wall time: 35.4 ms
Out[12]:
PDList(path=rips-5.6.pdgm)

さきほどよりかなり計算時間が早くなっています。

それでは maxvalue を指定した場合の PD をプロットしましょう。

In [13]:
pdlist56 = hc.PDList("rips-5.6.pdgm")
In [14]:
pd56_1 = pdlist56.dth_diagram(1)
pd56_1.histogram((5.1, 5.7)).plot(colorbar={"type": "log"})
pd56_2 = pdlist56.dth_diagram(2)
pd56_2.histogram((5.3, 5.9)).plot(colorbar={"type": "log"})
pd56_3 = pdlist56.dth_diagram(3)
pd56_3.histogram((5.3, 5.9)).plot(colorbar={"type": "log"})
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

PDの範囲の最大値が5.6になっています。maxvalueを指定しない場合と並べて描画してみましょう。 3x2でPDを並べていて，左がmaxvalueの指定なし，右が指定ありです。上から0次，1次，2次と並べています。

In [15]:
fig, axes = plt.subplots(3, 2, figsize=(14, 14)) 
pd1.histogram((5.1, 5.7)).plot(colorbar={"type": "log", "max": 10}, ax=axes[0][0])
pd56_1.histogram((5.1, 5.7)).plot(colorbar={"type": "log", "max": 10}, ax=axes[0][1])
pd2.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, ax=axes[1][0])
pd56_2.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, ax=axes[1][1])
pd3.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, ax=axes[2][0])
pd56_3.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, ax=axes[2][1])
No description has been provided for this image

maxvalueを指定したほうはペアが birth, death < 5.6 の範囲にあります。5.6よりも下の部分はどちらも同じ ヒストグラムになっています。

では切りとられた birth-death pairはどこに行ったのでしょうか？ 実はその一部は death time が∞であるbirth-death pairとなっています。 そのような birth-death pair の birth time は essential_births という属性に保持されています。

以下のコードで確認しましょう。

In [16]:
pd2.essential_births
Out[16]:
array([], dtype=float64)
In [17]:
pd56_2.essential_births
Out[17]:
array([5.59966087, 5.59966087, 5.59600019, 5.59598827, 5.59438658,
       5.59168386, 5.59153223, 5.59094238, 5.58698988, 5.58534145,
       5.58534145, 5.58519888, 5.58323002, 5.58296061, 5.58234787,
       5.58177948, 5.57721853, 5.57679319, 5.57349205, 5.57244968,
       5.57104015, 5.57104015, 5.56442404, 5.56233788, 5.56102133,
       5.55986166, 5.5575943 , 5.5541563 , 5.55370188, 5.55355883,
       5.5498476 , 5.53804016, 5.53659773, 5.53065109, 5.52408218,
       5.52092361, 5.49601793, 5.4945631 , 5.47305489])

maxvalue を指定しない場合、essential_births は空の配列になります。一方で、maxvalue を指定すると、いくつかの値が essential_births に現れます。

これは、maxvalue を境界として、birth が 5.6 以下かつ death が 5.6 以上のペアが、通常のパーシステンス図から essential_births に移動するためです。

また、birth と death の両方が 5.6 を超えるペアについては、maxvalue の制限により完全に除外され、情報として残りません。

以下のコードを使うと、death time が ∞（無限大）である birth-death ペアを可視化してプロットすることができます。

In [18]:
pd56_1.histogram((5.1, 5.7)).plot(colorbar={"type": "log", "max": 10}, plot_ess=True, title="1D")
pd56_2.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, plot_ess=True, title="2D")
pd56_3.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, plot_ess=True, title="3D")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

逆解析について¶

パーシステンス図上の各点は、元データに存在するホモロジー的構造（たとえば、1次元はリング、2次元は空洞など）に対応しています。これらの構造を抽出できれば、データの幾何的・トポロジカルな特徴をより深く理解することができ、解析に非常に有用です。

アルファ複体を用いる場合には、optimal volume や stable volume といった手法を使って、対応する構造を比較的効率よく取り出すことが可能です。

しかし、Vietoris-Rips 複体では単体の数が非常に多くなるため、これらの逆解析手法を適用するには計算コストが高く、実用的でない場合が多くあります。

ただし、1次のパーシステンス図に関しては、グラフの最短経路アルゴリズムを利用することで、Vietoris-Rips 複体においてもリング構造を比較的高速に抽出することが可能です。

本チュートリアルでは、この手法を活用します。これは optimal 1-cycle と呼ばれるアプローチです。

なお、このアルゴリズムで得られるリングは、あくまで近似的な構造であり、パーシステントホモロジーの情報を完全に反映しているわけではない点に注意が必要です。

optimal 1-cycle を計算するには、from_rips_filtration を呼び出す際に save_graph=True という引数を指定する必要があります。これにより、後でリング構造を抽出するためのグラフ情報が保存されます。

また、今回は 1 次のパーシステンス図のみを使用するため、maxdim=1 と設定します。

In [19]:
%%time
hc.PDList.from_rips_filtration(dmatrix, maxdim=1, save_graph=True, save_to="rips_with_graph.pdgm")

CPU times: user 12.1 ms, sys: 5 μs, total: 12.1 ms
Wall time: 11.8 ms
Out[19]:
PDList(path=rips_with_graph.pdgm)

1次のパーシステンス図を読み込んでプロットします．

In [20]:
pd1 = hc.PDList("rips_with_graph.pdgm").dth_diagram(1)
In [21]:
pd1.histogram().plot(colorbar={"type": "log"})
No description has been provided for this image

(5.3, 5.6) に一番近い点を抽出します。

In [22]:
pair = pd1.nearest_pair_to(5.3, 5.6)
pair
Out[22]:
Pair(5.312653064727783, 5.580554485321045)

Optimal 1-cycleを計算します．計算時間はそれほどかかりません．

In [23]:
optimal_1_cycle = pair.optimal_1_cycle()

boundary_points というメソッドを呼び出すとリング上の点が得られます．各点は番号で表現され，その番号は最初の点から0, 1, 2, ... と昇順で付けられます．

In [24]:
optimal_1_cycle.boundary_points()
Out[24]:
[45, 30, 68, 25, 82, 6]

注意¶

optimal 1-cycle を計算する際に、Pairs with the same birth time というエラーが発生することがあります。これは、距離行列の上半三角成分に同じ値が含まれている場合に起こります。この機能はすべての距離が異なることを前提としているため、まったく同じ距離が存在すると処理に失敗します。

このようなエラーが出る場合は、距離行列にごく小さなノイズを加えることで回避できます。その際は、行列の対称性を保ったままノイズを加えるようにしてください。

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