조모로디안의 알고리즘 구현
개요
조모로디안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 ↩︎