logo

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

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

概要

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

導出

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

アルゴリズム

Filtered_Complex.png

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

20220710_234704.png

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

20220710_234904.png

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

擬似コード

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

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

dd = REMOVEPIVOTROWS(σ)(\sigma)

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

ii = maxindexdd

  • Input: チェインddを受け取る。
  • Output: テーブルTTでチェインddに含まれるすべてのsimplexの中で最大のインデックスiiを返す。例えばmaxindex(abc)の場合はab55bc66ac99の中で最大の99を返す必要がある。

kk = dimdd

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

kk = dimσ\sigma

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

kk = deg(σ)(\sigma)

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

キーワード

  • MarkMarkσ\sigmaの形で書かれ、該当するシンプレックスσ\sigmamarkedtrueに変更する。
  • StoreStorejj and dd in T[i]T[i]の形で書かれ、T[i]T[i]Jに整数jjslotにチェインddを保存する。
  • RemoveRemovexx in ddの形で書かれ、チェインddにあるxx項を削除する。

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

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

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

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

function REMOVEPIVOTROWS(σ)(\sigma)

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


実装

与えられた例でのアルゴリズムの実行結果は次のようになるべきである。 L0={(0,),(0,1),(1,1),(1,2)}L1={(2,5),(3,4)} \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

全体のコード

読めばわかるが、元の論文のノーテーションをほぼそのままにしてコードを書いた。例えば Lk=Lk{(degσi,degσj)} 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 ↩︎