Vietoris-Rips filtrationのチュートリアル¶

このチュートリアルでは Vietoris-Rips filtrationを使ってパーシステント図を計算する例です。

Vietoris-Rips filtrationは以下のような利点があります。

  • 2点間の距離の情報(distance matrix)だけあれば計算できる
    • 遺伝子間の距離情報のようにユークリッド空間以外の空間でも計算できる
    • alpha filtrationは2次元、3次元でしか現在は計算できないが、より高次元の点でも計算できる

一方で以下のような欠点もあります。

  • 組み合わせ爆発により、計算コストが高い
    • 特に高次のホモロジーの計算が困難
    • そのため、計算する次数を制限することは必須
    • また、実用的な計算時間のためある距離より大きな辺は無視する必要がある場合がある
  • データの幾何的な情報を近似的にしか反映していない
    • 各点に球を置くやりかたとはけっこう大きなずれが生じる

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

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

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

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

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

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

distance matrix は 100x100 の行列です

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

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

homcloud.interface.PDList.from_rips_filtration でパーシステント図を計算します。引数の意味は以下の通りです

  • 最初の引数は距離行列
  • maxdim で最大次数を指定
  • save_to で保存先を指定
  • ここでは指定していませんが、maxvalue で取り扱う birth/death time の最大値を指定します
In [5]:
%%time
hc.PDList.from_rips_filtration(dmatrix, maxdim=3, save_to="rips.pdgm")

CPU times: user 1 s, sys: 83.8 ms, total: 1.09 s
Wall time: 1.08 s
Out[5]:
PDList(path=rips.pdgm)

計算時間は maxdim をどうするかでかなり変わります。

  • 2 - 一瞬
  • 3 - 5-10 sec
  • 4 - 2min
  • それ以上は試してない。maxdimに対して指数的に計算量が増大すると思われる

これらの計算は ripser というプログラムで計算されています。この ripser は Vietoris-Rips filtration による パーシステント図を計算するには(おそらく)最速のプログラムです。

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

以下のようにしてパーシステント図を読み込みます。

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

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

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

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

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

maxvalue の効果¶

さて、ここで PD を計算するときに maxvalue を指定すると何がおきるかを見ていきましょう。 まず、distance matrixの値がどのように分布しているか、histogramを見てみましょう。

numpy.triu_indicesを使うと距離行列の上半分だけ取り出すことができます。

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

5〜6.5の間に分布しているようです。

ここでは 5.6 を上限として計算しましょう。

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

CPU times: user 18.5 ms, sys: 0 ns, total: 18.5 ms
Wall time: 17.9 ms
Out[11]:
PDList(path=rips-5.6.pdgm)

さきほどよりかなり計算時間が早くなったことがわかると思います。

それでは PD をプロットしましょう。

In [12]:
pdlist56 = hc.PDList("rips-5.6.pdgm")
In [13]:
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 [14]:
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 というattributeに保持されています。では見てみましょう。

In [15]:
pd2.essential_births
Out[15]:
array([], dtype=float64)
In [16]:
pd56_2.essential_births
Out[16]:
array([5.59986973, 5.59731436, 5.5969243 , 5.59630442, 5.59630442,
       5.59550047, 5.59471798, 5.59387875, 5.59387875, 5.59291601,
       5.59288454, 5.59288454, 5.58777428, 5.58478022, 5.58478022,
       5.58323336, 5.58320761, 5.58302402, 5.57969189, 5.57728481,
       5.57637548, 5.57592154, 5.57303667, 5.57085276, 5.57085276,
       5.56939459, 5.56820822, 5.56689215, 5.56289721, 5.56282139,
       5.56219244, 5.56159353, 5.5589695 , 5.5589695 , 5.55742836,
       5.55691719, 5.55553007, 5.55448961, 5.55394936, 5.55394936,
       5.54576492, 5.54538584, 5.54364109, 5.53661442, 5.53520107,
       5.52763557, 5.50185394, 5.50104332, 5.49944115, 5.48618507,
       5.45888186])

このように、maxvalueを指定しなかった場合にはessential_birthsは空ですが、指定した場合にはいくつかの値が得られます。 birth が 5.6 以下で death が 5.6 以上の birth-death pair がここに移動します。birth, death がともに 5.6以上の場合には それらの birth-death pair は痕跡もなく消されます。

これらの death time が ∞ な birth-death pair を並べてプロットすることもできます。

In [17]:
pd56_1.histogram((5.1, 5.7)).plot(colorbar={"type": "log", "max": 10}, plot_ess=True)
pd56_2.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, plot_ess=True)
pd56_3.histogram((5.3, 5.9)).plot(colorbar={"type": "log", "max": 10}, plot_ess=True)
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 の利用でこれが可能になります．しかし Vietoris-Rips 複体は単体の個数が多いためそのような構造を抽出するのは難しい問題になります．しかし1次元のPDに関してのみはある種の近似により比較的簡単にそういった構造の抽出が可能です．ここではその構造(optimal 1-cycleと呼びます)を抽出する方法について説明しましょう．

optimal 1-cycleを計算するためには，from_rips_filtrationを呼ぶ際にsave_graph=Trueという引数を付ける必要があります．1次元PDしか使わないのでmaxdimも1で計算しましょう．

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

CPU times: user 8.49 ms, sys: 441 µs, total: 8.93 ms
Wall time: 8.15 ms
Out[18]:
PDList(path=rips_with_graph.pdgm)

1次元のPDを読み込んでプロットします．

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

(5.3, 5.6) に一番近い点を取り出します．

In [21]:
pair = pd1.nearest_pair_to(5.3, 5.6)
pair
Out[21]:
Pair(5.288362503051758, 5.525618553161621)

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

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

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

In [23]:
optimal_1_cycle.boundary_points()
Out[23]:
[93, 50, 48, 7, 3]

注意¶

optimal 1-cycle を計算するときに， Pairs with the same birth time というエラーが発生する場合があります．実は距離行列の対称でない要素が同じ値を持っているとこの機能はうまく動きません．このようなエラーが出場合には距離行列に小さなノイズを加えておくとうまくいくでしょう．対称性を保ったままノイズを加えてください．

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