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