logo

ジョモロジアンのアルゴリズムの実装 📂位相データ分析

ジョモロジアンのアルゴリズムの実装

概要

ZomorodianとCarlssonの論文「Computing Persistent Homology」で紹介されたアルゴリズムの擬似コードを説明し、実装する。1 抽象的な単体複体で作られたフィルター付き複体を受け取り、$\mathcal{P}$-インターバルをリターンし、コンピュータで扱いにくい永続的なモジュールの構築を省略し、行列のリダクションによって永続的なホモロジーを計算する。さらに、実際の実装では行列演算さえも使用しない。

導出

  • Zomorodianのアルゴリズムの導出: アルゴリズムの理論的な内容を完全に無視すると、擬似コードをいくら見てもアルゴリズムを理解することはできないだろう。完全に理解する程度でなくても、なぜ突然行列がなくなったのか、なぜマーキングのようなものが必要なのか、少なくともその程度は理解できるように勉強してから実装に取り組むべきである。

アルゴリズム

Filtered_Complex.png

アルゴリズムが実行される前にフィルター付き複体をデータとして受け取ったとする。

20220710_234704.png

データを保存し、アルゴリズムから得られた情報を記録するためのディクショナリやテーブル$T$を上記のように作成する。例えばJuliaのデータフレームでは、データ内の数字をepsilon、アルファベットが記された部分をsimplexに転写し、マーキングの有無を保存するブーリアンカラムmarked、チェインコンプレックスのチェインを保存するslot、計算過程で出てくる整数を保存するカラムJを追加する。

20220710_234904.png

  • Juliaは配列が$0$ではなく$1$から始まる言語であり、実装の便宜上インデックスを気にしない部分があるため、上のスクリーンショットのように数字が$1~2$ずつすべて異なる場合があるが、これは全く重要ではないので気にしないでほしい。
  • 注意すべきは、slotが集合でありながらもチェインとして両方の表現を自由に行き来することである。例えば、$\mathbb{Z}_{2}$上での計算で $$ a + (a - b) = 2 a - b = b \pmod{2} $$ のような計算が行われるが、これは次のような集合演算 $$ \left\{ a \right\} \cup \left\{ a, -b \right\} = \left\{ b \right\} $$ とも同じことである。代数的な演算での要素$0$は集合では「ないもの」として扱われ、それをそのまま受け入れなければならない。厳密で慎重な表記を好む人には不快かもしれないが、そこまで無理なことではないので、そのまま受け入れよう。
  • また、元の論文では一般的なフィールド$F$上でアルゴリズムを導出しており、すなわちすべての$q \in F$の逆元$q^{-1} \in F$が存在して $$ d = d - q^{-1} T[i] $$ のような計算を行うが、この実装ではバイナリフィールド$\mathbb{Z}_{2}$で十分なので、$q^{-1}$を別途計算せずに$A \Delta B := \left( A \cup B \right) \setminus \left( A \cap B \right)$と定義された$\Delta$に対して次のように代替した。 $$ d = d \Delta T[i] $$
  • $\deg$という表現が（プログラムの）関数としても出てくるし、インデックスとしても出てくるし、多項関数の次数としても出てくるし、非常に頻繁に出てくるので、$T$ではdegではなくepsilonと表記した。実際、トポロジカルデータアナリシスの単純なレベルでは、通常このカラムの値である半径$\varepsilon > 0$が大きくなるにつれてフィルター付き複体を構成することになるためである。
  • $T$はシンプレックスの次元に従って完全にソートされていると仮定し、それに従ってepsilonも部分順序を持つことを期待する。

擬似コード

$\left\{ L_{k} \right\}$ = COMPUTEINTERVALS$(K)$

  • Input: フィルター付き複体$K$を受け取る。フィルター付き複体には、少なくともどのタイミング$\deg \in \mathbb{Z}$にどのシンプレックス$\sigma$が追加されたかについての情報が必要である。
  • Output: $k = 0, \cdots , \dim K$に対する$\mathcal{P}$-インターバルの集合$L_{k}$の集合$\left\{ L_{k} \right\}_{k=0}^{\dim K}$を得る。
  • Side Effect: データが記録されたテーブル$T$のmarkedを変更する。

$d$ = REMOVEPIVOTROWS$(\sigma)$

  • Input: $k$次元のシンプレックス$\sigma$を受け取る。
  • Output: $(k-1)$次元のチェイン、つまり$k$次元のシンプレックス同士を演算したある$\mathsf{C}_{k-1}$の要素を得る。

$i$ = maxindex$d$

  • Input: チェイン$d$を受け取る。
  • Output: テーブル$T$でチェイン$d$に含まれるすべてのsimplexの中で最大のインデックス$i$を返す。例えばmaxindex(abc)の場合はabの$5$、bcの$6$、acの$9$の中で最大の$9$を返す必要がある。

$k$ = dim$d$

  • Input: $k$次元のチェイン$d$を受け取る。
  • Output: 整数$k$を返す。

$k$ = dim$\sigma$

  • Input: $k$次元のシンプレックス$\sigma$を受け取る。
  • Output: 整数$k$を返す。

$k$ = deg$(\sigma)$

  • Input: シンプレックス$\sigma$を受け取る。
  • Output: テーブル$T$でシンプレックス$\sigma$に対応する整数epsilonを返す。例えばdeg(cd)の場合はcdのepsilonが$2$なので$2$を返す必要がある。

キーワード

  • MarkはMark$\sigma$の形で書かれ、該当するシンプレックス$\sigma$のmarkedをtrueに変更する。
  • StoreはStore$j$ and $d$ in $T[i]$の形で書かれ、$T[i]$のJに整数$j$、slotにチェイン$d$を保存する。
  • RemoveはRemove$x$ in $d$の形で書かれ、チェイン$d$にある$x$項を削除する。

$\sigma^{i}$はテーブル$T$で$i$番目にあるsimplexであり、$m$は$T$の長さである。

function COMPUTEINTERVALS$(K)$
　　# 初期化
　　for $k \in 0:\dim K$
　　　　$L_{k} := \emptyset$
　　end for

　　for $j \in 0:(m-1)$
　　　　$d$ = REMOVEPIVOTROWS$\left( \sigma^{j} \right)$
　　　　if $d = \emptyset$
　　　　　　# $d$が空であることは、（ピボットではない）ゼロ列の候補である
　　　　　　mark $\sigma^{j}$
　　　　else
　　　　　　# $d$の次元を計算する必要があるため、すべての項のmaxでなければならない
　　　　　　# $d$は$\sigma^{j}$より一次元低いチェインであり、$i < j$しかあり得ない
　　　　　　$i$ = maxindex$d$
　　　　　　$k$ = dim$d$
　　　　　　$L_{k}$ = $L_{k} \cup \left\{ \left( \deg \sigma^{i}, \deg \sigma^{j} \right) \right\}$
　　　　end if
　　end for

　　for $j \in 0:(m-1)$
　　　　# まだマークされていない場合は、明らかにゼロ列である
　　　　if $\sigma^{j}$ ismarked and $T[j]$ isempty
　　　　　　$k$ = dim$d$
　　　　　　# $H_{k-1}$から$\sum^{\hat{e}_{i}} F[t]$に該当、無限大処理
　　　　　　$L_{k}$ = $L_{k} \cup \left\{ \left( \deg \sigma^{i}, \infty \right) \right\}$
　　　　end if
　　end for
　　return $\left\{ L_{k} \right\}$
end function

function REMOVEPIVOTROWS$(\sigma)$

　　$k$ = $\dim \sigma$
　　# $\partial abc = ab - bc + ca$とすると∂("abc") = ["ab", "bc", "ca"]
　　$d$ = $\partial_{k} \sigma$
　　Remove not marked $(k-1)$-dimensional simplex in $d$
　　while $d \ne \emptyset$
　　　　$i$ = maxindex$d$
　　　　if $T[i]$ isempty
　　　　　　break
　　　　end if
　　　　# $\mathbb{Z}_{2}$なので、symdiff（対称差）で代用
　　　　$d$ = $d \Delta T[i]$
　　end while
　　return $d$
end function

実装

与えられた例でのアルゴリズムの実行結果は次のようになるべきである。 $$ \begin{align*} L_{0} =& \left\{ \left( 0, \infty \right) , \left( 0,1 \right) , \left( 1,1 \right) , \left( 1,2 \right) \right\} \\ L_{1} =& \left\{ \left( 2,5 \right) , \left( 3,4 \right) \right\} \end{align*} $$

Juliaで実装した結果は次のようになる。インデックスが正確に異なる部分を除けば、正しく実装されていることが確認できる。

20220711_005505.png

全体のコード

読めばわかるが、元の論文のノーテーションをほぼそのままにしてコードを書いた。例えば $$ L_{k} = L_{k} \cup \left\{ \left( \deg \sigma^{i}, \deg \sigma^{j} \right) \right\} $$ のJulia風のコードはpush!(L_[k], (deg(σⁱ), deg(σʲ)))であるが、論文とほぼ同じように見えるようにL_[k] = L_[k] ∪ [(deg(σⁱ), deg(σʲ))]として実装した。

using DataFrames

data = DataFrame([
0 "a"
0 "b"
1 "c"
1 "d"
1 "ab"
1 "bc"
2 "cd"
2 "ad"
3 "ac"
4 "abc"
5 "acd"
], ["epsilon", "simplex"])
T = copy(data)
T[!, :"marked"] .= false
T[!, :"slot"] .= [[]]
T[!, :"J"] .= 0

dimK = 2
m = nrow(T)

test_σ = "abc"
dim(σ) = length(σ)

function deg(σ)
    return T.epsilon[findfirst(T.simplex .== σ)]
end
deg(test_σ)

function ∂(σ)
    k = dim(σ)
    return [σ[(1:k)[Not(t)]] for t = 1:k]
end
∂(test_σ)

function maxindex(chain)
    return (T.simplex .∈ Ref(chain)) |> findall |> maximum
end
maxindex(∂(test_σ))

function REMOVEPIVOTROWS(σ)
    k = dim(σ); d = ∂(σ)
    d = d[d .∈ Ref(T[T.marked,:simplex])] # Remove unmarked terms in ▷eq029◁
    while !(d |> isempty)
        i = maxindex(d)
        if T[i,:slot] |> isempty break end
        d = symdiff(d, T[i,:slot])
        # print("d in empty")
    end
    return d
end
REMOVEPIVOTROWS(test_σ)

L_ = [[] for k = 0:dimK]
for j0 = 0:(m-1)
    j = j0+1
    σʲ = T[j,:simplex]
    d = REMOVEPIVOTROWS(σʲ)
    if d |> isempty
        T[j,:marked] = true
    else
        i = maxindex(d); k = dim(σʲ)
        σⁱ = T[i,:simplex]
        T[i,[:J,:slot]] = j0,d
        L_[k] = L_[k] ∪ [(deg(σⁱ), deg(σʲ))]
    end
end
for j0 = 0:(m-1)
    j = j0+1
    σʲ = T[j,:simplex]
    if (T[j,:marked]) && (T[j,:slot] |> isempty) && (T[j,:J] |> iszero)
        k = dim(σʲ); L_[k] = L_[k] ∪ [(deg(σʲ), Inf)]
        print("j: $j")
    end
end

  1. Zomorodian. (2005). Computing Persistent Homology: ch4 ↩︎