조모로디안의 알고리즘 구현
개요
조모로디안zomorodian과 칼슨Carlsson의 논문 ‘Computing Persistent Homology’에서 소개된 알고리즘의 의사코드를 설명하고 구현한다1. 추상 심플리셜 컴플렉스로 만들어진 필터드 컴플렉스를 받아 $\mathcal{P}$-인터벌을 리턴하며, 컴퓨터로 다루기 어려운 퍼시스턴트 모듈 구축을 생략하고 행렬의 리덕션으로 지속성 호몰로지를 계산한다. 더 나아가, 실제 구현에서는 행렬 연산조차 사용하지 않는다.
유도
- 조모로디안의 알고리즘 유도: 장담컨대 알고리즘의 이론적인 내용을 완전히 무시한다면 의사코드를 아무리 봐도 알고리즘을 이해할 수 없을 것이다. 아주 완벽하게 이해하는 정도는 아니더라도 왜 갑자기 행렬이 없어졌는지, 왜 마킹 같은 게 필요한지 정도는 이해할 수 있을만큼은 공부를 하고 구현에 나서도록 하자.
알고리즘
알고리즘이 작동되기 이전에 필터드 컴플렉스를 데이터로 받았다고 하자.
데이터를 저장하고 알고리즘에서 나온 정보를 기록할 수 있는 딕셔너리 혹은 테이블 $T$ 를 위와 같이 작성하자. 예를 들어 줄리아 데이터프레임으로는 데이터에 있는 숫자를 epsilon
, 알파뱃이 적힌 부분을 simplex
으로 옮겨 적고 마킹 여부를 저장할 불리언 칼럼 marked
, 체인컴플렉스의 체인을 저장할 slot
과 계산 과정에서 나오는 정수를 저장하는 칼럼 J
를 추가하는 식이다.
- 줄리아는 배열이 $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(symmetric difference)으로 대체
$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*} $$
줄리아로 구현한 결과는 다음과 같다. 인덱스가 정확하게 틀린 부분을 제외하곤 제대로 구현되었음을 확인할 수 있다.
전체 코드
읽어보면 알겠지만 원래 논문의 노테이션을 거의 그대로 따라서 코드를 작성했다. 이를테면
$$
L_{k} = L_{k} \cup \left\{ \left( \deg \sigma^{i}, \deg \sigma^{j} \right) \right\}
$$
의 줄리아다운 코드는 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 $d$
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 ↩︎