Interactive Differential Equations

정의. 상미분방정식(ODE)은 미지 함수 \(x(t)\)와 그 도함수들 사이의 관계식이다. 가장 흔한 1계 형태는 \(\frac{dx}{dt}=f(t,x)\)이며, 최고차 도함수의 차수를 그 방정식의 계수(order)라 한다. "상(ordinary)"은 독립변수가 하나뿐이라는 뜻으로, 변수가 여럿이면 편미분방정식이 된다.

방정식이지만 답이 함수다. \(x^2-4=0\)의 답은 수 \(\pm 2\)지만, \(\dot{x}=2x\)의 답은 함수 \(x(t)=Ce^{2t}\)다. 게다가 상수 \(C\)마다 하나씩, 해가 하나가 아니라 함수의 무리다. 오른쪽 그림에서 여러 곡선이 동시에 그려지는 것이 그 때문이며, 그중 하나를 고르려면 다음 절의 초기조건이 필요하다.

국소 규칙이 전역 궤적을 만든다. \(f(t,x)\)는 각 점에서의 기울기만 알려준다. 그 기울기를 짧은 선분으로 평면에 촘촘히 그린 것을 기울기장(direction field)이라 하고, 해곡선은 이 선분들에 매 순간 접하며 지나가는 곡선이다. 방정식을 푸는 법을 몰라도 해가 어떤 모양일지는 이 그림만으로 읽어낼 수 있다는 점이 중요하다.

어디서 나오는가. 뉴턴의 제2법칙 \(m\ddot{x}=F\)가 2계 ODE이고, 개체군 모형 \(\dot{P}=kP\)는 "증가 속도가 현재 개체수에 비례한다"를 그대로 옮긴 1계 ODE다. 물리·화학·생물·경제에서 "변화율"이라는 말이 나오는 순간 미분방정식이 세워진다.

머신러닝에서. 학습률을 0으로 보낸 극한에서 경사하강법은 \(\dot\theta=-\nabla L(\theta)\)라는 ODE가 된다(경사흐름). Neural ODE는 신경망 자체를 \(f\)로 두어 \(\dot{h}=f_\theta(h,t)\)를 풀고, 그 해를 층 대신 사용한다. 그리고 이 글의 종착지인 확산 모델의 생성 과정도 결국 학습된 ODE 하나를 푸는 일이다.

\[\begin{aligned} \frac{dx}{dt} &= f(t, x) \\ \dot{x} = 2x &\;\Rightarrow\; x(t) = C e^{2t} \\ m\ddot{x} &= F, \qquad \dot{P} = kP \end{aligned}\]

프리셋으로 방정식을 바꿔보세요. 배경의 짧은 선분들이 기울기장이고, 색색의 곡선이 서로 다른 상수 C에 대응하는 해입니다. 모든 해곡선이 자기 위치의 선분에 정확히 접한다는 점을 확인해보세요 — 방정식이 말하는 것은 그것뿐입니다.

정의. 초기값 문제(IVP)는 미분방정식과 초기조건을 함께 주는 문제다: \(\dot{x}=f(t,x),\; x(t_0)=x_0\). \(n\)계 방정식이라면 \(x(t_0)\)부터 \(x^{(n-1)}(t_0)\)까지 \(n\)개의 조건이 필요하다.

기하적으로: 곡선 하나 고르기. 기울기장 위에는 서로 겹치지 않는 해곡선들이 빼곡히 깔려 있다. 초기조건 \((t_0,x_0)\)는 평면 위의 점 하나를 찍는 일이고, 그 점을 지나는 곡선이 유일한 해가 된다. 오른쪽에서 점을 드래그하면 해곡선이 그 점을 따라 미끄러진다.

존재성과 유일성. 아무 \(f\)에서나 되는 것은 아니다. 피카르–린델뢰프 정리는 \(f\)가 연속이고 \(x\)에 대해 립시츠 조건을 만족하면 국소적으로 해가 존재하고 유일하다고 말한다. 립시츠 조건이 깨지면 유일성이 무너질 수 있다 — \(\dot{x}=\sqrt{|x|},\,x(0)=0\)은 \(x\equiv 0\)도 해이고 \(x=t^2/4\)도 해다. 오른쪽 프리셋에서 이 경우를 직접 볼 수 있다.

유한 시간 폭발. 해가 존재해도 영원히 존재하지는 않을 수 있다. \(\dot{x}=x^2,\,x(0)=1\)의 해는 \(x=1/(1-t)\)로 \(t\to 1\)에서 무한대로 발산한다. 방정식은 멀쩡한데 해가 유한 시간에 사라지는 것이며, 수치 시뮬레이션에서 값이 NaN이 되는 이유 중 하나다.

머신러닝에서. Neural ODE에서 초기값은 입력이고 최종 시각의 값이 출력이다. 해의 유일성이 보장되기에 이 대응이 함수로서 잘 정의되며, 게다가 궤적이 서로 교차할 수 없다는 성질 덕분에 이 사상은 자동으로 가역이다 — 정규화 흐름의 연속시간 버전이 성립하는 근거가 정확히 이것이다.

\[\begin{aligned} \dot{x} &= f(t,x), \qquad x(t_0) = x_0 \\ \text{립시츠: } &|f(t,x)-f(t,y)| \le L|x-y| \\ &\Rightarrow \text{국소적으로 해가 유일} \end{aligned}\]

초기점을 드래그하면 그 점을 지나는 해곡선 하나가 선택됩니다. "유일성 붕괴" 프리셋에서는 같은 초기점에서 두 개의 서로 다른 해가 뻗어 나오고, "유한시간 폭발"에서는 해가 특정 시각에 수직으로 치솟습니다.

정의. \(\dot{x}=f(x)\)처럼 우변에 \(t\)가 명시적으로 나타나지 않으면 자율계(autonomous), \(\dot{x}=f(t,x)\)처럼 나타나면 비자율계다. 자율계는 "물리 법칙이 어제나 오늘이나 같다"는 상황이고, 비자율계는 외부에서 시간에 따라 변하는 힘이 들어오는 상황이다.

자율계는 시간 평행이동에 불변이다. \(x(t)\)가 해라면 \(x(t+c)\)도 해다. 그러니 "언제 출발했는가"는 중요하지 않고 "어디서 출발했는가"만 중요하다. 덕분에 시간축을 지워버리고 상태공간만 그려도 정보가 손실되지 않는데, 그 그림이 다음 절의 위상 초상이다.

궤적이 교차하지 않는다. 자율계에서 각 점의 화살표는 하나뿐이므로 두 궤적이 한 점에서 만나면 그 이후로 완전히 같아져야 한다. 따라서 서로 다른 궤적은 결코 교차하지 않는다. 반면 비자율계는 같은 위치라도 시각에 따라 화살표가 달라지므로, \(x\)-\(t\) 평면이 아닌 \(x\)만의 공간에 그리면 궤적이 교차하는 것처럼 보인다.

비자율계를 자율계로 바꾸기. 시간을 상태의 일부로 승격시키면 된다. \(s=t\)라는 변수를 추가해 \(\dot{s}=1,\;\dot{x}=f(s,x)\)로 쓰면 우변에 \(t\)가 사라져 자율계가 된다. 차원이 하나 늘어나는 대신 이론이 단순해지는 거래다.

머신러닝에서: 왜 시간을 입력으로 넣는가. 확산 모델의 노이즈 예측 신경망은 \(\epsilon_\theta(x, t)\)처럼 시각 \(t\)를 반드시 입력으로 받는다. 순방향 노이즈 스케줄이 시간에 따라 변하는 비자율계이기 때문이며, 같은 \(x\)라도 초기 단계인지 후기 단계인지에 따라 해야 할 일이 전혀 다르다. Neural ODE에서 \(f_\theta(h,t)\)에 \(t\)를 넣을지 말지가 표현력과 직결되는 것도 같은 이유다.

\[\begin{aligned} \dot{x} &= f(x) &&\text{자율계} \\ \dot{x} &= f(t,x) &&\text{비자율계} \\ \dot{s} = 1,\; \dot{x} &= f(s,x) &&\text{자율화} \end{aligned}\]

두 계를 번갈아 보세요. 자율계에서는 시각이 달라도 같은 높이에서 기울기가 똑같아 기울기장의 각 가로줄이 동일합니다. 비자율계에서는 같은 높이라도 시각마다 기울기가 달라 해곡선의 모양이 출발 시각에 따라 바뀝니다.

정의. 2차원 자율계 \(\dot{x}=f(x,y),\;\dot{y}=g(x,y)\)에 대해, 상태공간 \((x,y)\) 위에 벡터장 \((f,g)\)와 대표 궤적들을 그린 그림이 위상 초상이다. 앞 글의 벡터장이 여기서 그대로 재등장한다 — 다만 이제 그 흐름선에 "시간에 따른 운동"이라는 해석이 붙는다.

고정점. \(f=g=0\)인 점에서는 화살표가 사라진다. 거기 놓인 상태는 영원히 움직이지 않으므로 고정점(평형점)이라 한다. 계의 장기 거동은 대체로 "어느 고정점으로 끌려가는가"로 결정된다.

선형화가 분류를 준다. 고정점 근처에서 계를 야코비안으로 근사하면 \(\dot{\delta}\approx J\delta\)가 되고, \(J\)의 고유값이 종류를 결정한다. 실수 고유값이 둘 다 음수면 안정 마디(빨려 들어감), 둘 다 양수면 불안정, 부호가 엇갈리면 안장(한 방향으로는 끌리고 다른 방향으로는 밀려남)이다. 고유값이 복소수면 회전이 섞여 나선이 되고, 실수부의 부호가 감기는지 풀리는지를 정한다. 선형대수의 고유값 하나가 동역학의 운명을 읽어준다.

왜 고유값이 운명을 결정하는가. \(\dot\delta=J\delta\)의 해는 고유벡터 방향으로 \(\delta(t)=e^{\lambda t}v\) 꼴이다. 즉 고유벡터 방향에서는 운동이 단순한 지수함수 하나로 줄어든다. 그러니 \(\lambda<0\)이면 그 방향 성분이 지수적으로 사라지고(끌려 들어감), \(\lambda>0\)이면 지수적으로 커진다(밀려남). 일반적인 출발점은 고유벡터들의 합이므로 각 방향이 자기 \(\lambda\)에 따라 따로따로 줄거나 늘고, 결국 \(\lambda\)가 가장 큰(가장 덜 음수인) 방향이 장기적으로 그림을 지배한다. 안장점에서 궤적이 한 방향으로 접근하다 결국 다른 방향으로 튀어 나가는 것이 이 경쟁의 결과다.

복소 고유값이 왜 회전인가. \(\lambda=a\pm bi\)일 때 \(e^{\lambda t}=e^{at}(\cos bt \pm i\sin bt)\)이므로, 실수부 \(a\)는 크기를 조절하고 허수부 \(b\)는 \(\cos\)·\(\sin\)을 통해 각도를 돌린다. 두 역할이 곱해져 "돌면서 커지거나 작아지는" 나선이 나온다. 회전이 나타나는 조건은 기하적으로도 읽을 수 있다 — 실수 고유벡터가 있다는 것은 변환이 방향을 보존하는 직선이 존재한다는 뜻인데, 순수한 회전에는 그런 직선이 아예 없다. 그래서 실수 고유벡터가 없는 경우(고유값이 복소수인 경우)가 곧 회전이 섞인 경우다. \(a=0\)이면 크기가 변하지 않아 닫힌 타원 궤도가 되고, GAN 학습이 수렴하지 못하고 진동하는 것이 정확히 이 상황이다.

극한 순환. 고정점만이 종착지는 아니다. 어떤 계는 닫힌 궤도로 수렴하는데, 이를 극한 순환이라 한다. 심장 박동이나 포식자–피식자 진동처럼 스스로 유지되는 리듬이 이런 구조다.

머신러닝에서. 경사흐름 \(\dot\theta=-\nabla L\)의 고정점이 곧 임계점이고, 그 야코비안이 헤시안이다. 헤시안이 양정부호면 안정 마디(국소 최소), 부호가 섞이면 안장점이 되는데 — 앞 글의 임계점 판정이 동역학의 언어로 다시 쓰인 것이다. 한편 GAN처럼 두 목적함수가 경쟁하는 계는 야코비안이 대칭이 아니라 복소 고유값을 갖고, 그래서 수렴 대신 나선형 진동이 나타난다.

\[\begin{aligned} \dot{x} &= f(x,y), \quad \dot{y} = g(x,y) \\ \dot{\delta} &\approx J\delta \quad \text{(고정점 근처)} \\ \lambda_{1,2} &\;\Rightarrow\; \text{마디 · 안장 · 나선} \end{aligned}\]

그림을 클릭하면 그 점에서 출발한 궤적이 그려집니다. 프리셋마다 고정점의 성격이 달라지며, 리드아웃에 야코비안의 고유값과 그로부터 판정한 종류가 표시됩니다. 포식자–피식자 프리셋에서는 궤적이 닫힌 고리를 그리며 영원히 순환합니다.

정의. \(x_{n+1}=x_n+\Delta t\, f(t_n,x_n)\). 현재 위치에서의 기울기를 읽어 그 방향으로 \(\Delta t\)만큼 직진하고, 도착한 곳에서 다시 기울기를 읽는다. 이것을 반복한다.

테일러 1차 절단이다. \(x(t+\Delta t)=x(t)+\Delta t\,\dot{x}(t)+\frac{\Delta t^2}{2}\ddot{x}+\cdots\)에서 1차까지만 남긴 것이 오일러법이다. 앞 글의 테일러 전개가 여기서 알고리즘이 된다 — 그리고 버린 항이 곧 오차다.

국소 오차와 전역 오차. 한 걸음의 오차는 \(O(\Delta t^2)\)이지만, 구간 \([0,T]\)를 가려면 \(T/\Delta t\)걸음이 필요하므로 누적 오차는 \(O(\Delta t)\)가 된다. 그래서 오일러법을 1차 정확도 방법이라 한다. \(\Delta t\)를 절반으로 줄이면 오차도 대략 절반이 되는데, 이 느린 수렴이 더 좋은 방법을 찾게 만든 동기다.

기하적으로 왜 어긋나는가. 오일러법은 매 걸음 출발점의 기울기만 쓴다. 그런데 구간 안에서 기울기가 변하므로, 곡선이 휘는 쪽에서 항상 바깥으로 밀려난다. 오른쪽에서 \(\Delta t\)를 키우면 근사 궤적이 참 해의 바깥쪽으로 체계적으로 벗어나는 것을 볼 수 있다 — 무작위 오차가 아니라 편향이다.

머신러닝에서. DDIM 샘플러를 비롯한 확산 모델의 결정론적 샘플러들이 본질적으로 확률흐름 ODE에 대한 오일러법(또는 그 변형)이다. 샘플링 단계 수를 줄인다는 것은 \(\Delta t\)를 키운다는 뜻이고, 그러면 이 절단 오차가 그대로 이미지 품질 저하로 나타난다. DPM-Solver 같은 고차 샘플러가 등장한 이유가 바로 다음 절에 있다.

\[\begin{aligned} x_{n+1} &= x_n + \Delta t\, f(t_n, x_n) \\ \text{국소 오차} &= O(\Delta t^2) \\ \text{전역 오차} &= O(\Delta t) \quad \text{(1차 정확도)} \end{aligned}\]

Δt를 크게 잡으면 계단식 근사(주황)가 참 해(회색)에서 눈에 띄게 벗어납니다. 걸음마다 그려지는 짧은 접선이 "그 지점의 기울기를 읽어 직진한다"는 규칙을 그대로 보여줍니다. Δt를 절반으로 줄일 때 오차도 대략 절반이 되는지 리드아웃에서 확인해보세요.

정의 (RK4). 한 걸음마다 네 곳에서 기울기를 재고 가중평균한다: 시작점 \(k_1\), 중점을 \(k_1\)으로 추정해 잰 \(k_2\), 중점을 \(k_2\)로 다시 추정해 잰 \(k_3\), 그리고 끝점의 \(k_4\). 최종 갱신은 \(x_{n+1}=x_n+\frac{\Delta t}{6}(k_1+2k_2+2k_3+k_4)\)다.

왜 중점을 두 번 쓰는가. 오일러법의 결함은 구간 시작점의 기울기로 구간 전체를 대표한 데 있다. 구간 중간의 기울기를 쓰면 훨씬 나은 대표값이 되고(중점법, 2차), 여러 지점을 가중평균하면 테일러 전개의 더 높은 항까지 상쇄된다. 가중치 \((1,2,2,1)/6\)은 심슨 적분 규칙과 같은 계수인데, 우연이 아니라 같은 원리다.

차수의 위력. RK4는 전역 오차가 \(O(\Delta t^4)\)다. \(\Delta t\)를 절반으로 줄이면 오차가 16분의 1이 된다. 한 걸음에 함수 평가를 4번 하므로 걸음당 비용은 4배지만, 같은 정확도에 필요한 걸음 수가 비교할 수 없이 적어 실질적으로는 훨씬 싸다. 오른쪽에서 세 방법을 같은 \(\Delta t\)로 비교하면 이 차이가 몇 자릿수로 벌어진다.

공짜는 아니다. 고차 방법은 해가 매끄러울 때만 이득을 준다. \(f\)가 급격히 변하거나 불연속이면 고차 항의 가정이 깨져 차수가 무너진다. 또 다음 절에서 볼 강성(stiff) 문제에서는 차수가 높아도 \(\Delta t\)를 키울 수 없다 — 정확도와 안정성은 별개의 문제다.

머신러닝에서. Neural ODE는 내부적으로 적응형 RK 계열 솔버(Dormand–Prince 등)를 쓰며, 그래서 "층 수"가 고정되지 않고 해의 난이도에 따라 정해진다. 확산 모델의 DPM-Solver·DEIS 같은 고속 샘플러도 확률흐름 ODE의 특수한 구조(선형 항을 정확히 적분하고 비선형 항만 근사)를 이용한 고차 방법이며, 그 덕에 1000단계를 10~20단계로 줄인다.

\[\begin{aligned} k_1 &= f(t_n, x_n) \\ k_2 &= f(t_n + \tfrac{\Delta t}{2}, x_n + \tfrac{\Delta t}{2}k_1) \\ k_3 &= f(t_n + \tfrac{\Delta t}{2}, x_n + \tfrac{\Delta t}{2}k_2) \\ k_4 &= f(t_n + \Delta t, x_n + \Delta t\,k_3) \\ x_{n+1} &= x_n + \tfrac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}\]

같은 Δt로 세 방법을 비교해보세요. 아래 로그 스케일 오차 그래프에서 기울기가 곧 정확도 차수입니다 — 오일러는 1, 중점법은 2, RK4는 4의 기울기로 내려갑니다. Δt를 절반으로 줄일 때 RK4의 오차가 16분의 1이 되는지 확인해보세요.

시험 방정식. \(\dot{x}=\lambda x\;(\lambda<0)\)의 참 해는 \(x(t)=x_0e^{\lambda t}\)로 0에 수렴한다. 여기에 오일러법을 적용하면 \(x_{n+1}=(1+\lambda\Delta t)x_n\)이므로 \(x_n=(1+\lambda\Delta t)^n x_0\)이다. 이 값이 0으로 가려면 \(|1+\lambda\Delta t|<1\), 즉 \(\Delta t < 2/|\lambda|\)여야 한다.

경계를 넘으면 무슨 일이 일어나는가. \(\Delta t\)가 \(2/|\lambda|\)를 넘는 순간 배율의 절댓값이 1보다 커지고, 매 걸음 부호가 뒤집히며 크기가 커진다 — 참 해는 단조 감소하는데 수치해는 부호를 번갈아 튀며 폭발한다. 오른쪽에서 \(\Delta t\)를 임계값 너머로 밀면 이 현상이 즉시 나타난다. 정확도를 아무리 높여도 이 문제는 해결되지 않는다.

강성(stiffness). 한 계 안에 아주 빠른 성분과 아주 느린 성분이 함께 있으면(고유값의 크기 차이가 크면) 곤란해진다. 관심 있는 것은 느린 거동인데, step size는 가장 빠른 성분이 강제하는 안정성 한계에 묶인다. 이런 문제를 강성 문제라 하고, 명시적 방법으로는 사실상 풀 수 없다.

암시적 방법. 후진 오일러 \(x_{n+1}=x_n+\Delta t\,f(t_{n+1},x_{n+1})\)는 도착점의 기울기를 쓴다. 시험 방정식에서는 배율이 \(1/(1-\lambda\Delta t)\)가 되어 \(\lambda<0\)이면 \(\Delta t\)가 아무리 커도 절댓값이 1 미만이다(무조건 안정). 대신 매 걸음 방정식을 풀어야 해서 비용이 든다 — 안정성을 계산량으로 사는 셈이다.

머신러닝에서. 경사하강법의 학습률 상한이 정확히 같은 계산에서 나온다. 이차 손실에서 \(\theta_{n+1}=(I-\eta H)\theta_n\)이므로 안정 조건은 \(\eta<2/\lambda_{\max}(H)\)다. "학습률이 너무 크면 손실이 진동하며 발산한다"는 경험칙이 이 부등식 그 자체다. 그리고 앞 글의 조건수가 여기서 재등장한다 — \(\lambda_{\max}\)가 학습률을 묶고 \(\lambda_{\min}\)이 수렴 속도를 정하니, 조건수가 크다는 것은 곧 강성 문제라는 뜻이다.

\[\begin{aligned} \text{전진 오일러: } &|1+\lambda\Delta t| < 1 \;\Rightarrow\; \Delta t < 2/|\lambda| \\ \text{후진 오일러: } &\left|\tfrac{1}{1-\lambda\Delta t}\right| < 1 \;\; \forall \Delta t > 0 \\ \text{경사하강: } &\eta < 2/\lambda_{\max}(H) \end{aligned}\]

Δt를 임계값(빨간 선) 너머로 밀어보세요. 참 해(회색)는 얌전히 0으로 가는데 전진 오일러(주황)는 부호를 번갈아 튀며 폭발합니다. 후진 오일러(보라)는 같은 Δt에서도 멀쩡합니다. λ를 더 음수로 만들면 임계값이 좁아지는 것이 강성 문제입니다.

정의. 자율계 \(\dot{x}=f(x)\)의 흐름(flow)은 초기 상태를 시간 \(t\) 뒤의 상태로 보내는 사상 \(\Phi_t: x_0 \mapsto x(t)\)다. 궤적 하나가 아니라 공간 전체를 한꺼번에 옮기는 변환으로 본다는 점이 시점의 전환이다.

군 구조. 흐름은 \(\Phi_0=\mathrm{id}\)이고 \(\Phi_s\circ\Phi_t=\Phi_{s+t}\)를 만족한다 — \(t\)만큼 흘린 뒤 \(s\)만큼 더 흘리는 것은 \(s+t\)만큼 흘리는 것과 같다. 따라서 \(\Phi_t^{-1}=\Phi_{-t}\)로, 모든 흐름은 시간을 거꾸로 돌려 되돌릴 수 있다. 이 가역성이 이 글 전체에서 반복해 쓰인다.

부피는 어떻게 변하는가. 흐름을 따라 옮겨진 영역의 부피는 발산이 지배한다: \(\frac{d}{dt}\log V=\nabla\cdot f\). 발산이 0인 벡터장(비압축성)은 부피를 정확히 보존하고, 양이면 팽창시키며, 음이면 수축시킨다. 오른쪽에서 사각형 하나가 흐름을 따라 변형되는 모습으로 이것을 볼 수 있다.

머신러닝에서: Neural ODE. \(f\)를 신경망으로 두면 \(\Phi_1\)이 하나의 층처럼 작동한다. ResNet의 \(h_{n+1}=h_n+f_\theta(h_n)\)이 정확히 \(\Delta t=1\)인 오일러 한 걸음이라는 관찰이 출발점이었다 — 잔차 연결은 이산화된 흐름이었던 셈이다. 층 수 대신 적분 구간과 솔버 정밀도가 계산량을 정하고, 메모리는 층 수와 무관해진다(수반법).

연속 정규화 흐름. 흐름이 가역이므로 확률변수를 흘려보내 분포를 변환할 수 있다. 이때 로그밀도는 \(\frac{d\log p}{dt}=-\nabla\cdot f\)를 따르는데, 야코비안 행렬식 대신 대각합만 계산하면 된다는 점이 결정적이다. 마지막 절의 Probability Flow ODE가 정확히 이 구조를 갖는다.

\[\begin{aligned} \Phi_t &: x_0 \mapsto x(t) \\ \Phi_s\circ\Phi_t &= \Phi_{s+t}, \qquad \Phi_t^{-1} = \Phi_{-t} \\ \frac{d}{dt}\log V &= \nabla\cdot f, \qquad \frac{d\log p}{dt} = -\nabla\cdot f \end{aligned}\]

시간 t를 움직이면 격자 전체와 그 안의 사각형이 흐름을 따라 변형됩니다. 발산이 0인 장에서는 사각형이 찌그러져도 넓이가 유지되고, 수축하는 장에서는 넓이가 줄어듭니다. t를 음수로 돌리면 정확히 원래 자리로 돌아온다는 점 — 흐름의 가역성입니다.

정의. 편미분방정식(PDE)은 여러 변수의 미지 함수와 그 편미분들 사이의 관계식이다. 대표적으로 \(u(t,x)\)에 대해 \(\frac{\partial u}{\partial t}=\alpha\frac{\partial^2 u}{\partial x^2}\)(열방정식)처럼 쓴다.

ODE와 무엇이 다른가. ODE의 해는 곡선 하나지만 PDE의 해는 함수 전체가 시간에 따라 변형되는 것이다. ODE가 "점 하나의 운동"이라면 PDE는 "모양 전체의 진화"다. 그래서 초기조건도 값 하나가 아니라 초기 함수 \(u(0,x)\) 전체이고, 공간이 유한하면 경계조건까지 필요하다.

세 가지 원형. 2계 선형 PDE는 크게 셋으로 분류된다. 포물형(열방정식 \(u_t=\alpha u_{xx}\))은 매끄럽게 만들고 정보를 잃으며 시간을 거꾸로 돌리기 어렵다. 쌍곡형(파동방정식 \(u_{tt}=c^2u_{xx}\), 수송방정식)은 모양을 유지한 채 유한한 속도로 전파하며 가역적이다. 타원형(라플라스 방정식 \(\Delta u=0\))은 시간이 없는 평형 상태를 기술한다. 이 글에서 만날 것은 앞의 둘과 그 조합이다.

대부분 손으로 풀 수 없다. ODE보다 훨씬 심하다. 그래서 실무에서는 공간을 격자로 잘라 편미분을 차분으로 바꾸고, 결과로 나온 거대한 ODE 계를 시간에 대해 적분한다(선의 방법). 오른쪽 데모들이 정확히 이 방식으로 계산된다.

머신러닝에서. 이 글의 주인공인 Fokker–Planck 방정식이 PDE이며, 확률밀도가 시간에 따라 어떻게 변형되는지를 기술한다. 확산 모델을 이해한다는 것은 결국 이 PDE와 그것을 만들어내는 SDE의 관계를 이해하는 일이다. 한편 물리 방정식 자체를 신경망으로 푸는 PINN·Neural Operator 계열도 활발한 분야다.

\[\begin{aligned} \frac{\partial u}{\partial t} &= \alpha\frac{\partial^2 u}{\partial x^2} &&\text{포물형 (확산)} \\ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} &= 0 &&\text{쌍곡형 (수송)} \\ \Delta u &= 0 &&\text{타원형 (평형)} \end{aligned}\]

같은 초기 모양에 세 방정식을 각각 적용해 비교해보세요. 수송은 모양을 그대로 옮기고, 확산은 뭉개며, 파동은 좌우로 갈라집니다. 시간 슬라이더를 되돌려보면 확산만 정보를 잃어 원래 모양을 복원하지 못한다는 점이 드러납니다.

정의. \(\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0\). 상수 속도 \(c\)의 경우 해는 \(u(t,x)=u_0(x-ct)\)로, 초기 모양을 오른쪽으로 \(ct\)만큼 민 것이다. 모양은 조금도 변하지 않는다.

특성곡선. 왜 그런지는 관찰자를 흐름에 태워보면 알 수 있다. \(x(t)=x_0+ct\)를 따라 움직이며 \(u\)를 보면 \(\frac{d}{dt}u(t,x(t))=u_t+c\,u_x=0\)이므로 그 관찰자에게 \(u\)는 상수다. 이런 곡선을 특성곡선이라 하고, PDE가 특성곡선을 따라 ODE로 환원되는 이 구조가 쌍곡형 방정식의 핵심이다.

속도가 위치에 따라 다르면. \(c=c(x)\)면 특성곡선이 휘고, 모양이 늘어나거나 압축된다. 특성곡선끼리 충돌하면(빠른 뒤쪽이 느린 앞쪽을 따라잡으면) 해가 다치기(shock) 시작하는데, 이것이 비선형 보존법칙의 어려움이다.

정보를 잃지 않는다. 수송은 완전히 가역적이다. \(c\)의 부호를 뒤집어 같은 시간만큼 흘리면 정확히 원래 모양으로 돌아온다. 다음 절의 확산과 대비되는 결정적 성질이며, 이 차이가 마지막 절에서 "결정론적 생성이 가능한가"라는 질문으로 되돌아온다.

머신러닝에서. 수송방정식은 다음 절 연속방정식의 특수한 경우(밀도가 보존되고 속도장이 일정한 경우)다. 그리고 Flow Matching이 하는 일이 정확히 이것의 일반화다 — 데이터 분포를 목표 분포로 실어 나르는 속도장을 학습한다. 최적수송(optimal transport) 이론과 생성 모델이 만나는 지점이기도 하다.

\[\begin{aligned} \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} &= 0 \\ u(t,x) &= u_0(x - ct) \\ \frac{dx}{dt} = c &\;\Rightarrow\; u \text{ 는 그 곡선 위에서 상수} \end{aligned}\]

속도 c를 바꾸며 재생해보세요. 아래 시공간 그림에서 기울어진 줄무늬가 특성곡선이고, 그 기울기가 곧 속도입니다. c를 음수로 하면 왼쪽으로 흐르고, 시간을 되감으면 모양이 완벽하게 복원됩니다 — 수송은 정보를 잃지 않습니다.

정의. 밀도 \(\rho(t,x)\)와 속도장 \(v(t,x)\)에 대해 \(\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho v)=0\). 여기서 \(j=\rho v\)를 플럭스(흐름량)라 하며, 식은 "밀도의 시간 변화 + 플럭스의 발산 = 0"으로 읽는다.

왜 발산인가. 작은 상자를 하나 생각하자. 상자 안의 총량이 줄어드는 유일한 방법은 경계로 흘러나가는 것뿐이다. 경계를 통한 알짜 유출량이 바로 발산 \(\nabla\cdot j\)였으므로, 보존은 곧 \(\partial_t\rho=-\nabla\cdot j\)가 된다. 앞 글의 "샘과 싱크" 그림이 여기서 물리적 의미를 얻는다.

총량은 항상 보존된다. 전체 공간에 대해 적분하면 발산 항이 (경계에서 사라진다는 가정 아래) 0이 되어 \(\frac{d}{dt}\int\rho\,dx=0\)이다. 확률밀도에 적용하면 "전체 확률은 언제나 1"이라는 뜻이 된다 — 확률 글의 정규화 조건이 동역학 아래에서도 자동으로 유지되는 이유다.

수송방정식과의 관계. \(v\)가 상수이고 \(\nabla\cdot v=0\)이면 \(\nabla\cdot(\rho v)=v\cdot\nabla\rho\)가 되어 앞 절의 수송방정식으로 환원된다. 반대로 \(v\)가 위치에 따라 변하면 밀도가 압축되거나 희박해지는 항이 추가로 생긴다.

머신러닝에서: 여기가 다리다. 입자들이 ODE \(\dot{x}=v(x)\)를 따라 움직일 때 그 입자들의 밀도는 정확히 이 연속방정식을 따른다. 즉 "궤적의 언어"와 "밀도의 언어"를 잇는 사전이 연속방정식이다. 세 절 뒤의 "입자에서 밀도로"에서 이것을 직접 확인하고, 마지막 절의 Probability Flow ODE는 이 관계를 거꾸로 이용한다 — 원하는 밀도 변화를 주고 그것을 만들어내는 속도장을 찾는 것이다.

\[\begin{aligned} \frac{\partial \rho}{\partial t} + \nabla\cdot(\rho v) &= 0 \\ j &= \rho v \quad \text{(플럭스)} \\ \frac{d}{dt}\int \rho\,dx &= 0 \quad \text{(총량 보존)} \end{aligned}\]

속도장을 바꿔가며 재생해보세요. 밀도(파랑)가 속도장을 따라 실려가며 수렴하는 곳에서는 쌓이고 발산하는 곳에서는 옅어집니다. 아래 막대가 총 질량인데, 모양이 아무리 변해도 언제나 1로 유지되는 것을 확인해보세요.

정의. \(\frac{\partial u}{\partial t}=\alpha\,\Delta u\). \(u\)가 온도면 열방정식, 농도나 확률밀도면 확산방정식이라 부르지만 수학적으로 완전히 같은 식이다. \(\alpha>0\)는 확산계수다.

라플라시안이 하는 일. 앞 글에서 라플라시안은 "이 점의 값이 주변 평균보다 얼마나 낮은가"였다. 그러니 이 방정식은 "주변보다 낮으면 올라가고 높으면 내려간다"를 그대로 옮긴 것이다. 봉우리는 깎이고 골은 메워지며, 시간이 충분히 지나면 전부 평평해진다.

기본해는 가우시안이다. 초기 조건이 한 점에 몰린 델타 함수라면 해는 \(u(t,x)=\frac{1}{\sqrt{4\pi\alpha t}}e^{-x^2/4\alpha t}\)로, 분산이 \(2\alpha t\)로 선형 증가하는 가우시안이다. 폭이 \(\sqrt{t}\)에 비례해 퍼진다는 이 사실이 뒤에서 브라운 운동을 만날 때 그대로 재등장한다 — 우연이 아니라 같은 현상의 두 얼굴이다.

왜 잔주름이 먼저 사라지는가. 라플라시안이 "주변 평균과의 차이"였다는 것을 떠올리면 이유가 바로 보인다. 같은 높이의 물결이라도 촘촘한 물결은 봉우리에서 골까지 가는 거리가 짧아 훨씬 급하게 휘어 있고, 곧 주변 평균과의 차이가 크다. 차이가 크면 그만큼 빨리 깎이므로, 잔주름은 순식간에 뭉개지고 완만한 큰 굴곡만 오래 남는다. 이 직관을 수식으로 확인하면 파수 \(k\)인 성분 \(\sin kx\)의 2차 미분이 \(-k^2\sin kx\)이므로 감쇠율이 \(k^2\)에 비례하고, 따라서 \(e^{-\alpha k^2 t}\)로 죽는다 — 주파수가 2배면 감쇠는 4배 빠르다. 사진에 가우시안 블러를 걸면 미세한 질감이 먼저 날아가고 큰 덩어리 형태는 남는 것이 같은 현상이다.

비가역적이다. 확산은 정보를 파괴한다. 고주파 성분일수록 빠르게(파수 \(k\)에 대해 \(e^{-\alpha k^2 t}\)로) 사라지므로, 시간이 조금만 지나도 세부 구조는 복원 불가능해진다. 시간을 거꾸로 돌린 방정식 \(u_t=-\alpha\Delta u\)는 고주파가 폭발해 수치적으로 다룰 수 없다. 그런데 확산 모델은 바로 이 "역방향 확산"을 해낸다 — 어떻게 가능한지가 이 글 후반부의 이야기다(답을 미리 말하면, 데이터 분포에 대한 추가 정보인 스코어를 학습하기 때문이다).

머신러닝에서. 확산 모델의 순방향 과정은 데이터 분포를 이 방정식에 따라 뭉개 가우시안으로 만드는 것이다. 이미지 처리의 가우시안 블러가 열방정식을 짧은 시간 푼 것이고, 그래프 위에서는 그래프 라플라시안에 의한 확산이 GNN의 메시지 전달과 over-smoothing 현상으로 나타난다.

\[\begin{aligned} \frac{\partial u}{\partial t} &= \alpha\,\Delta u \\ u(t,x) &= \tfrac{1}{\sqrt{4\pi\alpha t}}e^{-x^2/(4\alpha t)} \quad \text{(점 초기조건)} \\ \hat{u}(t,k) &= \hat{u}_0(k)\,e^{-\alpha k^2 t} \quad \text{(고주파가 먼저 죽는다)} \end{aligned}\]

재생하면 뾰족한 봉우리들이 먼저 무너지고 이어서 넓은 구조가 사라집니다. 아래 주파수 그래프에서 높은 주파수가 지수적으로 빠르게 죽는 것이 보입니다 — 이것이 확산이 비가역인 이유이고, 되감기로 원래 모양이 복원되지 않는 이유입니다.

두 서술. 라그랑주 관점은 입자를 따라다닌다: \(\dot{x}_i=v(x_i,t)\). 오일러 관점은 고정된 자리에서 지나가는 것을 센다: \(\partial_t\rho+\nabla\cdot(\rho v)=0\). 같은 물리를 기술하지만 미지수가 각각 궤적과 밀도로 다르다.

번역 규칙. \(N\)개의 입자를 속도장 \(v\)를 따라 움직이면, 그 경험분포는 \(N\to\infty\)에서 연속방정식의 해로 수렴한다. 반대로 밀도가 주어지면 입자를 그 밀도에서 뽑아 흘려보내면 된다. 오른쪽 데모에서 위쪽 점 구름(입자)과 아래쪽 곡선(밀도)이 같은 속도장 아래에서 함께 움직이며 서로 일치하는 것을 볼 수 있다.

왜 두 언어가 다 필요한가. 밀도는 계 전체의 완전한 정보를 담지만 고차원에서 저장·계산이 불가능하다(결합분포의 차원 저주와 같은 문제다). 입자는 고차원에서도 다룰 수 있지만 유한 표본이라 잡음이 있다. 이론은 밀도로 세우고 계산은 입자로 한다는 것이 실무의 전략이며, 확산 모델이 정확히 그렇게 작동한다 — Fokker–Planck로 정당화하고 샘플 경로로 실행한다.

결정론이어도 밀도는 퍼진다. 개별 입자가 무작위성 없이 움직여도, 속도장이 발산을 가지면 밀도는 퍼지거나 모인다. 즉 밀도의 확산이 반드시 무작위성에서 오는 것은 아니다. 이 관찰이 마지막 절 Probability Flow ODE의 씨앗이다 — 무작위한 SDE와 똑같은 밀도 변화를 결정론적 ODE로 재현할 수 있는 이유가 여기에 있다.

머신러닝에서. 생성 모델은 이 번역을 양방향으로 쓴다. 학습할 때는 밀도(가능도, 스코어)를 다루고, 샘플링할 때는 입자(하나의 경로)를 흘려보낸다. Flow Matching은 원하는 밀도 경로를 먼저 정하고 그것을 실현하는 속도장을 회귀로 학습하는데, 이것이 가능한 근거가 바로 이 절의 대응이다.

\[\begin{aligned} \dot{x}_i &= v(x_i, t) &&\text{(입자, 라그랑주)} \\ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho v) &= 0 &&\text{(밀도, 오일러)} \\ \hat\rho_N = \tfrac{1}{N}\textstyle\sum_i \delta_{x_i} &\;\xrightarrow{N\to\infty}\; \rho \end{aligned}\]

입자 수 N을 늘리며 재생해보세요. 위쪽 입자들의 히스토그램이 아래쪽 연속방정식 해(주황)에 점점 겹쳐집니다. N이 작으면 울퉁불퉁하지만 커질수록 매끄러워지는 것 — 확률 글의 샘플링에서 본 수렴이 여기서 동역학과 만납니다.