시뮬레이션

폭발적 퍼콜레이션 — 수직으로 튀는 곡선도 사실은 연속이었다

2 min read

폭발적 퍼콜레이션은 무작위로 간선을 추가하되, 매번 후보 두 개 중 잇는 두 클러스터 크기의 곱이 작은 쪽만 골라 성장을 억누르는 그래프 생성 규칙이다. 2009년 Science는 이 규칙에서 거대 클러스터가 순식간에 튀어나온다며 불연속 상전이로 보고했다. 2011년 Riordan과 Warnke는 이 전이가 연속임을 증명했다.

화면 속에서 곡선이 수직으로 튀는 순간

왼쪽 패널에서 연구자가 화면 속 바닥을 기다 수직으로 치솟는 곡선을 가리키고 붉은 배지가 '수직으로 보인다'를 단다. 오른쪽 패널의 돋보기 안에서 같은 구간이 이어진 계단형 경사로 드러나고 배지가 '실제로는 연속'을 단다

그래프 시뮬레이션을 걸어 두고 최대 클러스터가 자라는 곡선을 지켜보는 상황을 떠올려 보자. 간선을 하나씩 넣는 동안 곡선은 한참 바닥에 붙어 있다. 그러다 어느 지점에서 화면을 거의 수직으로 가른다. 이때 대부분은 같은 판정을 내린다 — 여기가 불연속 전이다. 같은 판정이 2009년 Science에 실렸고, 2년 만에 뒤집혔다.

거대 클러스터(giant component)는 노드 수 $N$에 비례하는 크기로 자라는 연결 덩어리다. 노드가 100배 늘면 이 덩어리도 대략 100배 커진다. 아래에서 $C$는 그 노드 수이고, 실제로 보는 값은 점유율 $C/N$이다.

느리게 자라는 그래프, 억지로 늦춘 그래프

기준점부터 잡는다. 노드 $N$개에 간선을 무작위로 하나씩 잇는 것이 에르되시-레니(ER) 무작위 그래프다. 간선을 $m$개 넣었을 때 진행도를 $r = m/N$으로 잰다. 거대 클러스터는 $r = 1/2$에서 처음 나타난다. 나타나는 방식이 완만하다는 점이 핵심이다 — 임계점에서 점유율은 0에서 출발해 기울기 4로 연속적으로 올라간다. 진행도가 얼마일 때 곡선이 어느 높이에 있는지는 식 하나를 풀면 나온다.

(1) $$ \rho = 1 – e^{-2r\rho} $$

여기서 $\rho$는 $C/N$의 극한값이고, 진행도 $r$이 그대로 식에 들어간다. $r=1$을 넣고 식 1을 풀면 0.7968이 나온다.

Achlioptas product rule은 이 과정에 한 줄을 더한다. 매 스텝 후보 간선을 하나가 아니라 두 개 뽑고, 각 후보가 이을 두 클러스터의 크기를 곱한 다음, 곱이 작은 쪽만 실제로 잇는다. 큰 덩어리끼리 붙는 일을 계속 미루는 규칙이다. 후보를 하나 더 볼 뿐인데 성장 곡선의 모양이 통째로 달라진다.

union-find로 직접 돌려본 지연-폭발 곡선

규칙을 코드로 옮기는 절차는 짧다. 필요한 것은 클러스터의 현재 크기를 조회하고 두 클러스터를 합치는 union-find 하나뿐이다.

  1. 후보 간선 2개를 무작위로 뽑는다.
  2. 각 후보가 이을 두 클러스터의 현재 크기를 union-find로 조회해 곱을 구한다.
  3. 곱이 작은 후보를 채택하고 나머지는 버린다.
  4. 채택한 간선으로 union 연산을 수행하고 크기를 갱신한다.

곱을 비교하는 그 한 순간이 이 규칙의 전부이니, 스텝 하나만 떼어 놓고 본다. 크기 10·6·1인 세 덩어리에서 후보 두 개를 견줘 보면 선택이 어떻게 갈리는지 드러난다.

def find(p, x):
    while p[x] != x:
        p[x] = p[p[x]]
        x = p[x]
    return x

def union(p, size, a, b):
    ra, rb = find(p, a), find(p, b)
    if ra == rb:
        return
    if size[ra] < size[rb]:
        ra, rb = rb, ra
    p[rb] = ra
    size[ra] += size[rb]

# 노드 17개를 세 덩어리로: A(0~9, 10개) · B(10~15, 6개) · C(16, 1개)
parent = list(range(17))
size = [1] * 17
for i in range(1, 10):
    union(parent, size, 0, i)
for i in range(11, 16):
    union(parent, size, 10, i)

def sz(x):
    return size[find(parent, x)]

candidates = {"A-C": (0, 16), "A-B": (0, 10)}
for name, (u, v) in candidates.items():
    print(f"후보 {name}: {sz(u)} x {sz(v)} = {sz(u) * sz(v)}")

name = min(candidates, key=lambda k: sz(candidates[k][0]) * sz(candidates[k][1]))
u, v = candidates[name]
union(parent, size, u, v)
print(f"채택: {name}")
print(f"최대 클러스터: {max(sz(i) for i in range(17))}")
후보 A-C: 10 x 1 = 10
후보 A-B: 10 x 6 = 60
채택: A-C
최대 클러스터: 11

A-B를 골랐다면 최대 클러스터가 단번에 16이 됐을 자리에서, 규칙은 11만 허용했다. 이 억제가 매 스텝 쌓인 결과가 그림 1다.

ER 무작위 그래프와 Achlioptas product rule의 최대 클러스터 점유율 곡선을 겹쳐 그린 선 그래프
그림 1. N=100,000·실현 20개 평균. ER은 r=0.5부터 완만히 자라고, product rule은 tc≈0.888까지 억제된 뒤 급등한다.

우리 재현에서 ER의 $r=1$ 점유율은 0.7965였다. 식 1의 이론값과 셋째 자리까지 맞으니 구현은 건전하다. 같은 코드에 product rule만 얹으면 거대 클러스터의 등장이 $r \approx 0.888$까지 미뤄지고, 그 지점을 지나면 곡선이 거의 세로로 선다. $r=1$ 점유율은 오히려 ER보다 높은 0.9172다.

폭발이 어떻게 준비되는가

수직으로 보이는 구간에서 실제로 무슨 일이 벌어지는지는 최대 클러스터 하나만 봐서는 알 수 없다. 클러스터 크기 분포 전체를 시간축과 함께 펼쳐야 보인다.

클러스터 크기 구간을 가로축, 진행도 r을 세로축으로 두고 노드 점유율을 밝기로 표시한 히트맵
그림 2. N=65,536 단일 실현. 임계점 아래에서 중간 크기 클러스터에 질량이 쌓이고, 임계점에서 그 덩어리들이 한꺼번에 병합된다.

그림 2에서는 임계점 한참 전부터 중간 크기 클러스터들이 차곡차곡 자란다. product rule이 막는 것은 큰 덩어리끼리의 병합이지, 중간 덩어리가 생기는 일이 아니기 때문이다. 그래서 임계점에 이르면 비슷한 크기의 덩어리 여럿이 한꺼번에 합쳐진다. 폭발은 그 자리에서 갑자기 생긴 것이 아니다. 그 전 구간 내내 준비된 것이다.

2009년의 판정 도구: Δ 기준

Achlioptas·D’Souza·Spencer는 이 곡선의 성격을 눈이 아니라 숫자로 판정하려 했다. 그 잣대가 Δ 기준이다. 최대 클러스터가 $\sqrt{N}$을 마지막으로 밑돈 스텝을 $t_0$라 하자. 처음으로 $0.5N$을 넘은 스텝은 $t_1$이다. 그 사이에 들어간 간선 수를 잰다.

(2) $$ \Delta = t_1 – t_0 $$

논리는 단순하다. 전이가 연속이면 $\sqrt{N}$에서 $N/2$까지 가는 데 $N$에 비례하는 간선이 필요하다. 그래서 $\Delta/N$이 상수로 남는다. 반대로 이 값이 $N$을 키울수록 0으로 가면 극한에서는 폭이 없는 점프가 된다. Science 2009는 product rule에서 이 점프가 $2N^{2/3}$ 스텝 안에 끝난다고 보고했다.

우리 실측도 같은 그림을 낸다. 식 2를 그대로 재면 $N=4{,}096$에서 $\Delta/N$이 0.0891이다. 노드를 256배 늘려 $N=1{,}048{,}576$까지 가면 이 값이 0.0116으로 떨어진다. 9점 로그-로그 적합의 기울기는 −0.368, 즉 $\Delta \sim N^{0.632}$이니 선형 증가에서 뚜렷이 벗어난다. 문헌이 주장한 $N^{2/3}$ 스케일과 부합하는 결과다(지수 차 0.03). 같은 격자의 ER 대조군은 0.259에서 0.224로 거의 그대로다.

반전 1 — 그 잣대는 고전 퍼콜레이션도 오판한다

문제는 이 잣대가 통과시키는 대상에 있다. 리뷰 문헌은 같은 기준을 2차원 격자의 고전 본드 퍼콜레이션에 적용해 본 결과를 보고한다. 연속 전이라는 것이 오래전에 밝혀진 계다. 그런데 여기서도 $\Delta/N$이 시스템 크기를 키울수록 멱법칙으로 줄어든다.

기준대로라면 고전 퍼콜레이션도 불연속이 된다. 그럴 리가 없으니 틀린 것은 계가 아니라 기준이다.

$\Delta/N \to 0$은 불연속의 증거가 아니다. 명백히 연속인 계도 이 조건을 만족한다. 곡선이 급해 보인다는 관찰과 곡선이 실제로 끊겨 있다는 주장은 다른 층위이고, 둘을 잇는 다리는 유한한 시뮬레이션 안에 없다.

반전 2 — 정리가 뒤집었다, 그리고 착시의 정체

유한한 시뮬레이션 안에서 판정이 나지 않으면 밖에서 찾는 수밖에 없다. 2011년 Riordan과 Warnke가 이 문제를 증명 쪽에서 닫았다.

Achlioptas 과정의 상전이는 연속이다.

Riordan & Warnke, “Achlioptas process phase transitions are continuous”, Ann. Appl. Probab. 22(4):1450–1464 (2012) — 논문 제목

조건은 느슨하다. 매 스텝 뽑는 정점의 수가 $N$과 무관하게 고정돼 있기만 하면 된다. product rule은 후보 간선 2개, 즉 정점 4개를 뽑으므로 이 안에 들어간다. 그러니 $N$을 키울수록, 최대 클러스터가 무시할 만한 크기에서 $N$에 비례하는 크기로 몇 스텝 만에 뛸 확률은 0으로 간다.

그러면 왜 시뮬레이션에서는 그렇게 안 보였나. 임계점 근처의 성장 속도를 정하는 지수 $\beta$가 이 계열에서 극단적으로 작기 때문이다. 관련 모델의 수치해가 $\beta = 0.0555(1)$, 대략 $1/18$이다.

지수가 이 정도면 임계점 직후 곡선은 사실상 수직으로 솟았다가 완만해진다. 연속이지만, 눈으로도 유한 표본으로도 불연속과 구별되지 않는다.

구별이 얼마나 안 되는지는 점프 크기의 감쇠 속도로 잰다. 한 스텝에 최대 클러스터가 늘어난 양의 최댓값을 $N$으로 나눈 값이 그림 3다.

최대 단일스텝 점프 크기를 N에 대해 로그-로그로 그린 산점도
그림 3. N=2^12~2^20, 로그-로그 축. 실현 수는 N에 따라 5~20개다.

$N$을 256배 키우는 동안 이 비율은 2할가량 줄어드는 데 그쳤다. 진짜 불연속이라면 상수로 버텨야 하고, 평범한 연속 전이라면 눈에 띄게 무너져야 한다. 실제로는 어느 쪽으로도 기울지 않는다. 개별 점은 실현 수 차이로 요동이 있으니, 읽을 것은 스윕 전체의 감소 총량이다. 문헌이 20억 노드를 돌려도 시뮬레이션만으로는 판별할 수 없다고 보고한 것이 이 상황이다.

판정법은 남았다 — 유한크기 스케일링

이 논쟁이 남긴 것은 "수직으로 보인다"를 대체할 절차다. 한 크기에서 곡선 모양을 보는 대신, $N$을 바꿔 가며 점프가 어떻게 변하는지를 본다. 점프 비율이 멱법칙으로 죽으면 연속이고, 크기를 바꿔도 상수로 남으면 불연속 쪽이다. 판정의 대상이 곡선의 인상에서 지수 하나로 옮겨 간다.

자기 시뮬레이션에서 이 판정을 해 보려면 $N$을 두세 배씩 늘려 가며 같은 양을 다시 재고 로그-로그로 찍으면 된다. 기울기가 0에서 떨어지지 않는다면 아직 판정할 때가 아니다. 위의 union-find 코드에 스텝별 최대 증가량만 기록해 두면 그대로 돌아간다.

경계 조건은 남는다. 정리가 요구하는 것은 매 스텝 표본 수가 고정이라는 조건이고, 이 조건을 깨는 규칙은 실제로 불연속 전이를 만든다.

그러니 남는 교훈은 폭발이 없었다는 것이 아니다. 계가 연속인지 불연속인지는 $N$을 키우는 동안 점프 비율이 무너지느냐 버티느냐에서 갈린다. 화면 속 수직 구간은 그 감쇠를 확인하기 전까지 판정을 보류한 채로 남는다.

더 읽기