6장 — 해밀토니안과 위상공간

leapfrog: 넓이를 지키는 적분법

지금까지는 해밀턴 흐름을 정확한 해로 따라갔다. 하지만 HMC처럼 공을 굴리는 계산을 컴퓨터로 하려면 시간을 Δt 간격으로 잘라야 한다. 잘라서 한 걸음씩 옮기는 계산도 해밀턴 흐름처럼 넓이와 에너지를 지킬까?

오일러 방법과 leapfrog

가장 먼저 떠오르는 방법은 지금의 속도로 위치와 운동량을 한꺼번에 옮기는 오일러 방법(Euler method)이다. 다른 방법은 운동량을 반 걸음 옮기고, 새 운동량으로 위치를 한 걸음 옮긴 뒤, 새 위치에서 운동량을 다시 반 걸음 옮기는 leapfrog 방법 (개구리가 번갈아 뛰어넘듯 위치와 운동량을 엇갈려 갱신하는 적분법)이다. 두 방법으로 그네(E = p²/2 + q²/2)를 Δt = 0.1로 적분하고, 에너지와 한 걸음이 넓이를 몇 배로 바꾸는지 비교해 보자.

import numpy as np

def grad_U(q):            # U(q) = q²/2 (그네, m = ω = 1)
    return q

def euler(q, p, dt):      # 위치와 운동량을 동시에 옛 값으로 갱신
    return q + dt * p, p - dt * grad_U(q)

def leapfrog(q, p, dt):   # 반 걸음 차기 → 한 걸음 이동 → 반 걸음 차기
    p = p - 0.5 * dt * grad_U(q)
    q = q + dt * p
    p = p - 0.5 * dt * grad_U(q)
    return q, p

H = lambda q, p: 0.5 * p**2 + 0.5 * q**2
dt = 0.1
for step in (euler, leapfrog):
    q, p = 1.0, 0.0
    for n in range(1, 1001):
        q, p = step(q, p, dt)
        if n in (100, 1000):
            print(step.__name__, n, round(H(q, p), 4))
    # 한 걸음 변환의 야코비 행렬식 (넓이 배율), 수치 미분으로
    e = 1e-6
    J = np.array([np.subtract(step(1 + e, 0.3, dt), step(1, 0.3, dt)),
                  np.subtract(step(1, 0.3 + e, dt), step(1, 0.3, dt))]).T / e
    print(step.__name__, "det J =", round(np.linalg.det(J), 6))
# euler 100 1.3524
# euler 1000 10479.5778
# euler det J = 1.01
# leapfrog 100 0.4996
# leapfrog 1000 0.4997
# leapfrog det J = 1.0

처음 에너지는 0.5인데, 오일러 방법으로는 100걸음 뒤에 1.3524, 1000걸음 뒤에는 10479.58이 되어 그네가 하늘 끝까지 날아간다. 이유는 행렬식에 있다. 이 그네에서 오일러 방법의 한 걸음은 (q, p)에 행렬 [[1, Δt], [−Δt, 1]]을 곱하는 것이고, 그 행렬식은 1 + Δt² = 1.01이다. 한 걸음마다 넓이가 1% 늘어나니 상태는 원을 따라 돌지 못하고 바깥으로 나선을 그리며, 에너지도 정확히 매 걸음 1.01배가 되어 0.5 × 1.01^100 = 1.3524가 된다.

leapfrog의 행렬식이 1인 까닭

반면 leapfrog의 행렬식은 정확히 1이다. 같은 기울기를 쓰고 갱신 순서만 바꿨을 뿐인데 왜 넓이가 늘지 않을까? 세 단계가 각각 운동량만 바꾸거나 위치만 바꾸는 전단 변환(한 변수를 다른 변수의 함수만큼 밀어 주는 변환, shear map)이고, 그런 변환의 야코비 행렬은 대각선이 모두 1인 삼각행렬이라 행렬식이 1이기 때문이다. 이렇게 넓이를 정확히 지키는 적분법을 심플렉틱 적분법 (위상공간의 넓이 구조를 보존하는 적분법, symplectic integrator)이라 한다. 에너지는 0.5에서 조금씩 흔들리지만 10만 걸음을 가도 0.49875와 0.5 사이를 벗어나지 않는다. 사실 이 그네에서 leapfrog는 조금 변형된 에너지 p²/(2(1 − Δt²/4)) + q²/2를 정확히 보존하므로 진짜 에너지도 멀리 벗어날 수 없다.

역사: 여러 번 다시 발견된 적분법

leapfrog는 한 사람이 만든 방법이 아니다. 오랜 시간에 걸친 궤도를 손으로, 뒤에는 컴퓨터로 따라가야 했던 사람들이 저마다 같은 계산법에 이르렀다. 1790년대 초 프랑스의 천문학자 들랑브르는 로그표와 천문표를 계산하며 이미 이 방법을 썼고, 1907년 노르웨이의 수학자 스퇴르메르는 오로라를 설명하려고 태양에서 날아온 전기 띤 알갱이가 지구 자기장 속에서 그리는 궤도를 이 방법을 다듬은 셈법으로 계산했다. 1909년 영국의 코웰과 크로멜린은 이듬해 돌아올 핼리 혜성의 궤도를 이 방법으로 계산했고, 1960년대에는 프랑스의 물리학자 베를레가 분자들이 서로 밀고 당기는 액체를 컴퓨터로 흉내 내며(분자 동역학) 다시 찾아냈다. 그래서 이 방법은 오늘날 스퇴르메르–베를레 방법이라고도 불린다. 궤도를 오래 따라가야 한다는 같은 고민이 같은 답을 여러 번 불러낸 셈이다.

ML에서: 넓이를 지키는 결합층

이산 normalizing flow의 초기 모델인 NICE(Dinh 외, 2014)의 덧셈 결합층(additive coupling layer)은 입력을 두 부분으로 나눠 한쪽을 그대로 두고 다른 쪽에 신경망 출력을 더하므로, 야코비 행렬식이 1인 부피 보존 변환이다. 이것은 leapfrog의 반 걸음, 곧 위치는 그대로 두고 운동량에 −(Δt/2)∇U(q)를 더하는 전단 변환과 같은 모양이다. leapfrog 한 걸음은 신경망 대신 정해진 함수(위치에너지의 기울기와 운동량)를 쓴 결합층 세 개를 쌓은 것이라고 볼 수 있다.

문제 12. 걸음을 줄이면 해결될까

그네(ℋ = p²/2 + q²/2)를 q = 1, p = 0에서 출발시켜 오일러 방법(q ← q + Δt·p, p ← p − Δt·q를 동시에)으로 적분하면 Δt = 0.1에서 에너지가 매 걸음 1.01배가 된다. (가) 걸음을 Δt = 0.01로 줄여 시간 10까지 가면 에너지는 얼마인가? 시간 1000까지 가면? (나) 걸음 크기는 0.1 그대로 두고 두 줄의 순서만 바꿔, 운동량을 먼저 갱신한 뒤 새 운동량으로 위치를 옮기면 에너지는 어떻게 되는가? 이때 한 걸음은 넓이를 몇 배로 바꾸는가? (풀어 본 뒤 위젯 3의 「문제 12 불러오기」로 확인해 보자.)

걸음이 너무 커서 그래요. Δt = 0.01로 줄여서 같은 시간 10까지 1000걸음 돌렸더니 0.5526이에요. 훨씬 낫죠. 더 줄이면 0.5에 붙을 거예요.

한 걸음마다 에너지가 정확히 1 + Δt²배가 되니까, 시간 t까지 가면 (1 + Δt²)^(t/Δt), 대략 e^(Δt·t)배야. Δt = 0.01이면 시간 10에서는 e^0.1 ≈ 1.105배라 괜찮아 보이는데, 시간 1000까지 가면 e^10배야.

돌려 볼게요… 10만 걸음 뒤에 11007.7이요. 어, 그네가 또 날아갔어요.

걸음을 줄이면 날아가는 시점을 늦출 뿐 막지는 못해요. HMC나 분자 시뮬레이션처럼 오래 굴려야 하는 곳에서는 치명적이죠. 민준 학생, 코드에서 한 줄만 바꿔 볼까요? 운동량을 먼저 갱신하고, 위치는 새 운동량으로 옮기는 거예요.

p ← p − Δt·q를 먼저 하고 q ← q + Δt·p를 새 p로… Δt = 0.1로 10만 걸음 돌렸는데 에너지가 0.476과 0.526 사이에서만 왔다 갔다 해요. 한 줄 순서만 바꿨는데요?

민준아, 행렬식을 봐. 운동량만 바꾸는 단계는 [[1, 0], [−Δt, 1]], 위치만 바꾸는 단계는 [[1, Δt], [0, 1]]이라 둘 다 행렬식이 1이야. 동시에 갱신한 원래 방법은 [[1, Δt], [−Δt, 1]]이라 1 + Δt²이었고.

그래요. 넓이가 매 걸음 1%씩 늘면 상태는 나선을 그리며 바깥으로 나갈 수밖에 없어요. 넓이를 정확히 지키는 적분법이 심플렉틱 적분법이고, 방금 민준 학생이 만든 것과 leapfrog가 대표적이에요.

그런데 넓이를 지킨다고 에너지까지 지키는 건 아니잖아요. 0.476에서 0.526까지 흔들리는 걸 보면요.

맞아요. 넓이 보존과 에너지 보존은 다른 성질이에요. 다만 심플렉틱 적분법은 원래 에너지와 조금 다른 「변형된 에너지」를 거의 정확히 보존해서, 진짜 에너지가 한쪽으로 흘러가지 않고 그 근처에서 흔들리기만 해요. 그리고 leapfrog는 앞뒤가 대칭이라 흔들리는 폭도 훨씬 작아요. 같은 Δt에서 0.49875에서 0.5 사이였죠.

과제에서 시뮬레이션이 터지면 늘 학습률 줄이듯이 Δt부터 줄였는데, 적분법 자체를 바꿔야 하는 경우가 있네요.