수치 미분

차분 몫: 할선의 기울기로 접선을 어림하기

도함수의 정의(CalculusBook 3장)에서 극한을 떼어 내면 가장 간단한 근사가 나옵니다. 오른쪽 점을 쓰면 전진, 왼쪽 점을 쓰면 후진, 양쪽을 대칭으로 쓰면 중앙 차분입니다.

아래 시뮬레이터에서 점을 곡선 위로 끌고 $h$를 줄여 보세요. 할선이 접선으로 돌아 들어가는 속도가 방법마다 다릅니다.

테일러 전개로 읽는 오차 차수

예제 9.1의 "반으로 줄이면 반", "반으로 줄이면 4분의 1"은 테일러 전개(CalculusBook 12장: 나머지 항)에서 바로 나옵니다. $f$가 충분히 매끄러우면

첫 식에서 $f(x)$를 빼고 $h$로 나누면 $D_+f = f' + \frac h2 f'' + O(h^2)$입니다. 두 식을 빼면 짝수 차 항이 모두 지워져 $f(x+h) - f(x-h) = 2hf' + \frac{h^3}{3}f''' + O(h^5)$, 곧 $D_0 f = f' + \frac{h^2}{6}f''' + O(h^4)$입니다. 대칭으로 놓은 점이 $h$의 홀수 거듭제곱 항을 지운 것이 중앙 차분의 비밀입니다.

예제 9.1로 확인해 봅시다. $h = 0.01$에서 전진 오차 예측은 $\frac{0.01}{2}e^{\xi} \approx 0.005\times 2.73 = 0.0136$이고 실제 $0.01364$입니다. 중앙 오차 예측은 $\frac{10^{-4}}{6}e \approx 4.53\times10^{-5}$, 실제와 세 자리가 맞습니다. 오차의 상수까지 예측된다는 점이 중요합니다. 다음 절의 최적 간격 계산과 리처드슨 외삽이 모두 이 상수를 이용합니다.

절단 오차 대 반올림 오차: 최적의 h

범인은 1장의 자릿수 상쇄입니다. $f(x+h)$와 $f(x)$는 각각 상대 오차 $u = \varepsilon/2 \approx 1.1\times10^{-16}$ 정도로 저장됩니다. $h$가 작으면 두 값이 거의 같아서 뺄셈의 결과에는 믿을 수 있는 자리가 몇 개 남지 않고, 그 작은 차이를 다시 작은 $h$로 나누면 남은 오차가 $1/h$배로 커집니다. $f$를 계산하는 프로그램 자체의 오차가 $\delta$라면 그 효과는 더 큽니다.

$f = e^x$, $x = 1$이면 $|f| = M_2 = M_3 = e$라서 $h_+^* = \sqrt{2\varepsilon} = 2.1\times10^{-8}$, $E_+^* \approx 5.7\times10^{-8}$, $h_0^* = (1.5\varepsilon)^{1/3} = 6.9\times10^{-6}$, $E_0^* \approx 6.5\times10^{-11}$입니다. 경험칙으로 기억하면 1차 공식은 $h \approx \sqrt{\varepsilon} \approx 10^{-8}$에서 유효 숫자 약 8자리, 2차 공식은 $h \approx \sqrt[3]{\varepsilon} \approx 6\times10^{-6}$에서 약 10~11자리를 얻습니다($x$와 $f$의 크기에 맞춰 $h$를 $\sqrt{\varepsilon}\,\max(1, |x|)$처럼 비례시킵니다). 배정밀도 16자리를 다 얻지는 못한다는 것이 차분 몫의 본질적 한계입니다. 2계 도함수 $\bigl(f(x+h) - 2f(x) + f(x-h)\bigr)/h^2$은 반올림 오차가 $1/h^2$로 커지므로 $h^* \sim \varepsilon^{1/4} \approx 10^{-4}$이고 8자리 정도가 한계입니다.

아래 그래프의 가로축은 $h$, 세로축은 오차, 둘 다 로그 눈금입니다. 오른쪽에서는 절단 오차가 기울기 1(또는 2)로 줄고, 왼쪽에서는 반올림 오차가 기울기 −1로 커집니다. 둘이 만나는 바닥이 최적 $h$입니다.

고계 도함수와 차분 스텐실

점의 배치를 $x + s_jh$($s_j$는 정수 또는 실수 오프셋)라 하고 $f^{(d)}(x) \approx h^{-d}\sum_j c_j f(x + s_jh)$를 찾습니다. 각 $f(x + s_jh)$를 테일러 전개해 $h^k f^{(k)}$의 계수를 모으면 $\sum_j c_j s_j^k/k!$입니다. $k = d$일 때 1, 나머지 $k$일 때 0이 되게 하면 됩니다.

이 연립방정식은 5장의 보간과 같은 이야기입니다. 점들을 지나는 보간 다항식을 만들어 그것을 $d$번 미분한 것이 바로 이 스텐실입니다. 그래서 점을 $m$개 쓰면 $m - 1$차 다항식까지 정확히 미분합니다.

function stencil(s[1..m], d)            // 오프셋 s, 도함수 차수 d < m
  for k = 0 .. m-1
    for j = 1 .. m:  V[k][j] = s[j]^k
    rhs[k] = (k == d) ? k! : 0
  c = solve(V, rhs)                     // 유리수로 풀면 정확한 분수 계수
  q = m
  while sum_j c[j] * s[j]^q == 0:  q = q + 1
  return c, order = q - d, lead = sum_j c[j] * s[j]^q / q!

아래 계산기에서 오프셋을 골라 보세요. 분수 계수, 차수, 오차 상수를 구하고, $f(x) = e^{x}$에 적용한 실제 오차를 로그–로그로 그려 기울기를 확인합니다.

리처드슨 외삽: 서로 다른 h를 섞어 차수 올리기

근사 $A(h)$가 참값 $A$에 대해 $A(h) = A + c_1h^p + c_2h^{p'} + \cdots$($p' \gt p$)로 전개된다고 합시다. $A(h/2) = A + c_1h^p/2^p + \cdots$이므로 $2^pA(h/2) - A(h) = (2^p - 1)A + O(h^{p'})$. 양변을 $2^p - 1$로 나누면 $c_1$을 몰라도 주항이 지워집니다.

이 생각은 영국의 기상학자 리처드슨(Lewis Fry Richardson)이 20세기 초 미분방정식의 차분 해를 개선하려고 쓴 것으로 알려져 있고, 다음 장의 롬베르크 적분, 11장의 외삽 적분법이 모두 같은 구조입니다. 외삽값과 원래 값의 차이 $|A^{(1)}(h) - A(h/2)|$는 $A(h/2)$의 오차 추정으로도 쓰입니다.

잡음 데이터의 미분: 차분은 잡음을 키운다

측정값이 $y_i = f(t_i) + e_i$, $e_i$가 서로 독립이고 표준편차 $\sigma$인 잡음이라 합시다. 중앙 차분 $(y_{i+1} - y_{i-1})/(2\Delta t)$에서 잡음 부분의 표준편차는 $\sqrt{2}\sigma/(2\Delta t) = 0.707\,\sigma/\Delta t$입니다. 위 문제의 값을 넣으면 $0.707$ rad/s, 각속도가 1 rad/s 정도라면 70%의 잡음입니다. 측정을 더 촘촘히 하면 오히려 $1/\Delta t$로 나빠집니다. 정리 9.3의 반올림 오차가 $\varepsilon$에서 $\sigma$로 커진 판입니다. 미분은 나쁜 조건의 문제입니다(2장). 진폭 $\sigma$인 고주파 잡음 $\sigma\sin(\omega t)$는 미분하면 진폭 $\sigma\omega$가 되기 때문입니다.

처방은 두 가지입니다. 첫째, 차분 간격을 측정 간격보다 넓게 잡습니다. 간격 $H = m\Delta t$의 중앙 차분은 잡음 $\sigma/(\sqrt2 H)$, 절단 $H^2|f'''|/6$이라 정리 9.3과 같은 계산으로 $H^* = (3\sigma/(\sqrt2 M_3))^{1/3}$ 근처가 최적입니다. 둘째, 두 점만 쓰지 않고 창 안의 모든 점에 다항식을 최소제곱으로 맞춘 뒤(7장) 그 다항식을 미분합니다. 이를 사비츠키–골레이 필터(Savitzky–Golay filter)라 부릅니다.

$m = 10$(21점)이면 $\sum k^2 = 770$이라 잡음이 $\sigma/(27.7\Delta t) = 0.036\,\sigma/\Delta t$, 2점 중앙 차분의 약 20분의 1입니다. 대신 $\sum k^4/(6\sum k^2) = 50666/4620 = 11.0$이라 $\Delta t = 0.01$에서 절단 오차가 $1.1\times10^{-3}|f'''|$로 늘어납니다. 잡음과 해상도(빠른 변화를 따라가는 능력)를 바꾸는 거래입니다.

복소 스텝 미분: 뺄셈 없는 차분

$f$가 실수에서 실숫값을 갖는 해석 함수라면 $f(x + ih) = f(x) + ihf'(x) - \frac{h^2}{2}f''(x) - i\frac{h^3}{6}f'''(x) + \cdots$입니다. 실수부와 허수부가 갈라지므로 허수부만 취하면 $\operatorname{Im} f(x + ih) = hf'(x) - \frac{h^3}{6}f'''(x) + \cdots$.

$f = e^x$, $x = 1$이면 $\operatorname{Im} e^{1+ih}/h = e\sin h/h$이고, $h = 0.1$에서 $2.7137536$(오차 $-4.5\times10^{-3}$, 중앙 차분과 크기가 같고 부호가 반대), $h = 10^{-8}$ 이하에서는 $2.718281828459045$, 배정밀도의 $e$와 마지막 자리까지 같습니다. 스콰이어와 트랩의 시험 함수 $f(x) = e^x/\sqrt{\sin^3x + \cos^3x}$의 $x = 1.5$에서도 $h = 10^{-20}$의 복소 스텝은 $4.05342789389862$로 참값과 일치하고, 중앙 차분은 최적 $h$에서도 상대 오차 $10^{-10}$ 근처에 머뭅니다.

자동 미분: 연쇄 법칙을 기계적으로

자동 미분(automatic differentiation, AD)은 기호 미분도 수치 미분도 아닙니다. 프로그램의 실행 기록(계산 그래프)에 연쇄 법칙(CalculusBook 4장)을 한 연산씩 적용합니다. 방향이 두 가지입니다.

전진 모드와 이중수

모든 중간값 $v$에 그 도함수 $\dot v = dv/dx$를 짝지어 $(v, \dot v)$로 들고 다닙니다. 덧셈은 $(a, \dot a) + (b, \dot b) = (a + b,\ \dot a + \dot b)$, 곱셈은 $(ab,\ \dot ab + a\dot b)$, 함수는 $g(a, \dot a) = (g(a),\ g'(a)\dot a)$입니다. 이 규칙은 $\epsilon^2 = 0$인 새 수 $\epsilon$을 도입해 $a + \dot a\epsilon$으로 쓰면 저절로 나옵니다. $(a + \dot a\epsilon)(b + \dot b\epsilon) = ab + (\dot ab + a\dot b)\epsilon$. 이런 수를 이중수(dual number)라 합니다. 복소 스텝에서 $i^2 = -1$ 때문에 남던 $h^2$ 항이, $\epsilon^2 = 0$에서는 완전히 사라진다고 보면 됩니다.

후진 모드와 역전파

입력이 $n$개, 출력이 하나인 함수 $f: \mathbb{R}^n \to \mathbb{R}$(신경망의 손실 함수가 그렇습니다)의 기울기를 전진 모드로 구하면 입력마다 한 번, $n$번 돌려야 합니다. 후진 모드는 반대로 갑니다. 먼저 값을 앞으로 계산해 그래프에 저장하고(전진 패스), 출력에서 시작해 각 중간값 $v$에 대한 수반 변수(adjoint) $\bar v = \partial f/\partial v$를 연쇄 법칙으로 거꾸로 전파합니다(후진 패스). $v$가 여러 곳에 쓰였으면 각 경로의 기여를 더합니다: $\bar v = \sum_{w\,:\,v \to w}\bar w\,\partial w/\partial v$.

이 도구가 쓰이는 곳

세 가지 대표 상황에서 고르는 도구가 다릅니다. 속도 추정(엔코더, GPS, 영상 추적)은 잡음 데이터의 미분이라 평활화와 필터가 핵심입니다. 민감도 해석(설계 변수에 대한 출력의 도함수, 불확실성 전파)은 함수가 비싼 시뮬레이션이라 계산 횟수가 중요하고, 변수가 적으면 중앙 차분이나 복소 스텝, 많으면 수반(adjoint) 방법이나 후진 모드 AD를 씁니다. 신경망 학습은 입력(가중치)이 수십억 개, 출력(손실)이 하나인 극단적인 경우라 후진 모드 AD, 곧 역전파 외에는 선택지가 없습니다.