ジョモロジアンのアルゴリズムの実装
概要
ZomorodianとCarlssonの論文「Computing Persistent Homology」で紹介されたアルゴリズムの擬似コードを説明し、実装する。1 抽象的な単体複体で作られたフィルター付き複体を受け取り、-インターバルをリターンし、コンピュータで扱いにくい永続的なモジュールの構築を省略し、行列のリダクションによって永続的なホモロジーを計算する。さらに、実際の実装では行列演算さえも使用しない。
導出
- Zomorodianのアルゴリズムの導出: アルゴリズムの理論的な内容を完全に無視すると、擬似コードをいくら見てもアルゴリズムを理解することはできないだろう。完全に理解する程度でなくても、なぜ突然行列がなくなったのか、なぜマーキングのようなものが必要なのか、少なくともその程度は理解できるように勉強してから実装に取り組むべきである。
アルゴリズム
アルゴリズムが実行される前にフィルター付き複体をデータとして受け取ったとする。
データを保存し、アルゴリズムから得られた情報を記録するためのディクショナリやテーブルを上記のように作成する。例えばJuliaのデータフレームでは、データ内の数字をepsilon
、アルファベットが記された部分をsimplex
に転写し、マーキングの有無を保存するブーリアンカラムmarked
、チェインコンプレックスのチェインを保存するslot
、計算過程で出てくる整数を保存するカラムJ
を追加する。
- Juliaは配列がではなくから始まる言語であり、実装の便宜上インデックスを気にしない部分があるため、上のスクリーンショットのように数字がずつすべて異なる場合があるが、これは全く重要ではないので気にしないでほしい。
- 注意すべきは、
slot
が集合でありながらもチェインとして両方の表現を自由に行き来することである。例えば、上での計算で のような計算が行われるが、これは次のような集合演算 とも同じことである。代数的な演算での要素は集合では「ないもの」として扱われ、それをそのまま受け入れなければならない。厳密で慎重な表記を好む人には不快かもしれないが、そこまで無理なことではないので、そのまま受け入れよう。 - また、元の論文では一般的なフィールド上でアルゴリズムを導出しており、すなわちすべてのの逆元が存在して のような計算を行うが、この実装ではバイナリフィールドで十分なので、を別途計算せずにと定義されたに対して次のように代替した。
- という表現が(プログラムの)関数としても出てくるし、インデックスとしても出てくるし、多項関数の次数としても出てくるし、非常に頻繁に出てくるので、では
deg
ではなくepsilon
と表記した。実際、トポロジカルデータアナリシスの単純なレベルでは、通常このカラムの値である半径が大きくなるにつれてフィルター付き複体を構成することになるためである。 - はシンプレックスの次元に従って完全にソートされていると仮定し、それに従って
epsilon
も部分順序を持つことを期待する。
擬似コード
= COMPUTEINTERVALS
- Input: フィルター付き複体を受け取る。フィルター付き複体には、少なくともどのタイミングにどのシンプレックスが追加されたかについての情報が必要である。
- Output: に対する-インターバルの集合の集合を得る。
- Side Effect: データが記録されたテーブルの
marked
を変更する。
= REMOVEPIVOTROWS
= maxindex
- Input: チェインを受け取る。
- Output: テーブルでチェインに含まれるすべての
simplex
の中で最大のインデックスを返す。例えばmaxindex(abc)
の場合はab
の、bc
の、ac
のの中で最大のを返す必要がある。
= dim
- Input: 次元のチェインを受け取る。
- Output: 整数を返す。
= dim
- Input: 次元のシンプレックスを受け取る。
- Output: 整数を返す。
= deg
- Input: シンプレックスを受け取る。
- Output: テーブルでシンプレックスに対応する整数
epsilon
を返す。例えばdeg(cd)
の場合はcd
のepsilon
がなのでを返す必要がある。
キーワード
Mark
はMark
の形で書かれ、該当するシンプレックスのmarked
をtrue
に変更する。Store
はStore
and in の形で書かれ、のJ
に整数、slot
にチェインを保存する。Remove
はRemove
in の形で書かれ、チェインにある項を削除する。
はテーブルで番目にあるsimplex
であり、はの長さである。
function COMPUTEINTERVALS
# 初期化
for
end for
for
= REMOVEPIVOTROWS
if
# が空であることは、(ピボットではない)ゼロ列の候補である
mark
else
# の次元を計算する必要があるため、すべての項のmax
でなければならない
# はより一次元低いチェインであり、しかあり得ない
= maxindex
= dim
=
end if
end for
for
# まだマークされていない場合は、明らかにゼロ列である
if ismarked
and isempty
= dim
# からに該当、無限大処理
=
end if
end for
return
end function
function REMOVEPIVOTROWS
=
# とすると∂("abc") = ["ab", "bc", "ca"]
=
Remove
not marked
-dimensional simplex
in
while
= maxindex
if isempty
break
end if
# なので、symdiff(対称差)で代用
=
end while
return
end function
実装
与えられた例でのアルゴリズムの実行結果は次のようになるべきである。
Juliaで実装した結果は次のようになる。インデックスが正確に異なる部分を除けば、正しく実装されていることが確認できる。
全体のコード
読めばわかるが、元の論文のノーテーションをほぼそのままにしてコードを書いた。例えば
の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
Zomorodian. (2005). Computing Persistent Homology: ch4 ↩︎