보간·근사·수치적분

보간

왜 필요한가. 불규칙한 시각에 측정한 신호를 균일한 격자로 옮겨야 FFT를 쓰거나 다른 신호와 나란히 놓을 수 있다. 시뮬레이터가 격자점에서만 준 값을 임의의 위치에서 읽어야 할 때도 있다. 둘 다 점들 사이의 값을 채우는 문제다.

정의. 서로 다른 노드(node) $x_0, \dots, x_n$과 값 $y_i = f(x_i)$가 주어졌을 때, 모든 노드를 정확히 지나는 함수 $p$, 즉 $p(x_i) = y_i$인 $p$를 찾는 일을 보간(interpolation)이라 한다. $p$를 $n$차 이하 다항식으로 제한하면 그런 다항식은 정확히 하나 있다. 노드 범위 $[\min x_i, \max x_i]$ 안의 값을 구하면 보간, 밖의 값을 구하면 외삽(extrapolation)이다.

오차 공식. $f$가 $n+1$번 미분 가능하면, 보간 다항식의 오차는 어떤 점 $\xi$(노드 범위와 $x$를 포함하는 구간 안)에 대해

이다. 오른쪽은 두 부분으로 나뉜다. 앞의 분수는 함수가 얼마나 구불구불한지를, 뒤의 곱 $\prod (x - x_i)$(노드 다항식)는 노드가 어떻게 놓였는지를 나타낸다. 등간격 노드에서 노드 다항식은 구간 가운데에서 작고 양 끝에서 크다. 그리고 $x$가 노드 범위를 벗어나면 모든 인수 $(x - x_i)$가 같은 부호로 커지므로 곱이 폭발한다. 외삽이 위험한 이유가 이 식에 그대로 있다.

Runge 현상. 매끄러워 보이는 함수 $f(x) = 1/(1 + 25x^2)$를 $[-1, 1]$의 등간격 노드로 보간하면, 노드를 늘릴수록 양 끝의 진동이 커진다. 이를 Runge 현상이라 한다. 이 함수는 복소평면의 $x = \pm i/5$에 극을 가져서 고계 도함수가 $n!$보다 빠르게 커지고, 등간격 노드 다항식은 양 끝에서 그것을 상쇄할 만큼 작지 않다.

대응은 두 가지다. 첫째, 노드를 양 끝에 몰아 배치한다. Chebyshev 노드 $x_j = \cos\big(\pi(2j+1)/(2n+2)\big)$는 단위원 위의 등간격 점을 $x$축에 내린 것으로, 노드 다항식의 최대 크기를 가능한 가장 작게($2^{-n}$) 만든다. 둘째, 다항식 하나로 전체를 덮지 않고 구간마다 낮은 차수의 다항식을 이어 붙인다. cubic spline은 이웃한 노드 사이마다 3차식을 쓰고, 노드에서 값·1계·2계 도함수가 연속이도록 이은 함수다. 차수가 고정되어 있으므로 노드를 늘려도 진동이 커지지 않고, 매끄러운 함수에 대해 오차가 노드 간격 $h$의 4제곱에 비례해 준다.

import numpy as np
from scipy.interpolate import BarycentricInterpolator, CubicSpline

def runge(x):
    return 1.0 / (1.0 + 25.0 * x**2)

x_test = np.linspace(-1, 1, 2001)                      # (2001,)
for n in [5, 10, 15, 20]:
    x_eq = np.linspace(-1, 1, n + 1)                   # (n+1,) 등간격
    x_ch = np.cos(np.pi * (2 * np.arange(n + 1) + 1) / (2 * n + 2))   # (n+1,) Chebyshev
    err_eq = np.max(np.abs(BarycentricInterpolator(x_eq, runge(x_eq))(x_test) - runge(x_test)))
    err_ch = np.max(np.abs(BarycentricInterpolator(x_ch, runge(x_ch))(x_test) - runge(x_test)))
    err_sp = np.max(np.abs(CubicSpline(x_eq, runge(x_eq))(x_test) - runge(x_test)))
    print(f"n = {n:2d}  등간격 {err_eq:9.2e}  Chebyshev {err_ch:8.2e}  spline {err_sp:8.2e}")

n =  5  등간격  4.33e-01  Chebyshev 5.56e-01  spline 4.32e-01
n = 10  등간격  1.92e+00  Chebyshev 1.09e-01  spline 2.20e-02
n = 15  등간격  2.11e+00  Chebyshev 8.31e-02  spline 3.09e-02
n = 20  등간격  5.98e+01  Chebyshev 1.53e-02  spline 3.18e-03

등간격 보간의 최대 오차는 $n = 20$에서 60에 이른다. 같은 수의 Chebyshev 노드는 오차를 꾸준히 줄이고, 등간격 노드 그대로의 cubic spline은 더 빨리 줄인다. 여기서 다항식 계산에 쓴 BarycentricInterpolator는 보간 다항식을 계수로 펼치지 않고 무게중심(barycentric) 공식으로 계산해 반올림 오차에 강하다. 즉 위의 60은 반올림 탓이 아니라 등간격 보간 다항식 자체의 오차다.

보간 안과 밖. 같은 다항식이라도 노드 범위 안과 밖에서의 성질은 전혀 다르다. $[0, 1]$에서 잡음(표준편차 0.01)이 섞인 측정값 11개를 10차 다항식과 cubic spline으로 보간하고, 구간 안팎에서 참값과 비교한다.

rng = np.random.default_rng(0)
x_in = np.linspace(0, 1, 11)                           # (11,) 측정 구간 [0, 1]
y_in = np.sin(2 * np.pi * x_in) + 0.01 * rng.normal(size=11)   # 측정 잡음 0.01
for x_out in [0.5, 0.95, 1.05, 1.1, 1.2]:
    p10 = BarycentricInterpolator(x_in, y_in)(x_out)  # 10차 다항식
    sp = CubicSpline(x_in, y_in)(x_out)
    print(f"x = {x_out:4.2f}  참값 {np.sin(2 * np.pi * x_out):7.3f}  10차 다항식 {p10:7.3f}  spline {sp:7.3f}")

x = 0.50  참값   0.000  10차 다항식   0.004  spline   0.004
x = 0.95  참값  -0.309  10차 다항식  -0.309  spline  -0.322
x = 1.05  참값   0.309  10차 다항식   0.140  spline   0.325
x = 1.10  참값   0.588  10차 다항식  -0.295  spline   0.651
x = 1.20  참값   0.951  10차 다항식  -7.771  spline   1.196

구간 안($x \le 0.95$)에서는 두 방법 모두 잡음 수준의 오차를 낸다. 구간 끝에서 0.05만 나가도 10차 다항식은 참값의 절반도 맞추지 못하고, 0.2 나가면 $-7.8$로 부호까지 틀린다. 잡음 0.01이 노드 다항식의 큰 인수와 곱해져 증폭된 결과다. spline은 덜 나쁘지만 역시 구간 밖에서는 오차가 거리와 함께 커진다. 어떤 보간법도 외삽을 보장하지 않는다. 차이는 얼마나 빨리 나빠지느냐뿐이다.

어떤 보간법을 쓸까. 실무의 기본 선택지는 다음과 같다.

농도·확률처럼 음수가 될 수 없거나 단조여야 하는 양을 보간할 때 cubic spline은 노드 사이에서 음수나 역전을 만들 수 있다. 이때는 PCHIP을 쓴다.

근사

왜 필요한가. 측정값에는 잡음이 있다. 모든 점을 정확히 지나는 보간 함수는 잡음까지 따라가므로, 점 사이와 바깥에서 참함수와 크게 어긋난다. 잡음 있는 데이터에서 원하는 것은 점을 지나는 함수가 아니라 점들에 가까운 단순한 함수다.

정의. 기저 함수 $\phi_0, \dots, \phi_m$($m < n$)의 결합 $p(x) = \sum_k c_k \phi_k(x)$ 중에서 제곱 오차의 합

을 최소로 하는 계수 $\mathbf{c}$를 고르는 것을 최소제곱 근사(least-squares approximation)라 한다. $\mathbf{V}$는 $(n+1) \times (m+1)$ 행렬이고, 기저가 단항식 $1, x, x^2, \dots$이면 Vandermonde 행렬이다. 보간은 $m = n$인 특수한 경우로, 잔차가 0이 된다.

직관. 차수 $m$은 모델의 유연성이다. 너무 낮으면 참함수의 모양을 따라가지 못하고(편향), 너무 높으면 잡음을 따라간다(분산). 데이터 잔차는 차수를 올릴수록 줄어들기만 하므로, 잔차만 보고 차수를 고르면 항상 가장 높은 차수를 고르게 된다.

기저의 선택. 같은 다항식 공간이라도 어떤 기저로 표현하느냐에 따라 $\mathbf{V}$의 조건수가 크게 달라진다. 단항식 $x^k$는 $[-1, 1]$에서 $k$가 커질수록 서로 비슷한 모양(양 끝에서만 크고 가운데에서 평평)이 되어 열들이 거의 평행해진다. Chebyshev 다항식 $T_k(x) = \cos(k \arccos x)$나 Legendre 다항식 같은 직교다항식(orthogonal polynomials)은 서로 직교하도록 만든 기저라 이 문제가 훨씬 덜하다. 9.2장에서 본 것처럼 조건수가 크면 같은 데이터에서도 계수가 반올림 오차에 민감해진다.

rng = np.random.default_rng(1)
x = np.linspace(-1, 1, 21)                             # (21,)
y = np.exp(x) + 0.05 * rng.normal(size=x.size)         # 잡음 표준편차 0.05
x_fine = np.linspace(-1, 1, 1001)                      # (1001,)
for deg in [2, 4, 8, 20]:
    p = np.polynomial.Chebyshev.fit(x, y, deg)
    rms_fit = np.sqrt(np.mean((p(x) - y) ** 2))
    err = np.max(np.abs(p(x_fine) - np.exp(x_fine)))
    print(f"차수 {deg:2d}  데이터 잔차 RMS {rms_fit:.4f}  참함수와 최대 오차 {err:.4f}")
V_mono = np.vander(x, 21, increasing=True)             # (21, 21) 단항식 기저
V_cheb = np.polynomial.chebyshev.chebvander(x, 20)     # (21, 21) Chebyshev 기저
print(f"cond: 단항식 {np.linalg.cond(V_mono):.1e}  Chebyshev {np.linalg.cond(V_cheb):.1e}")

차수  2  데이터 잔차 RMS 0.0411  참함수와 최대 오차 0.0828
차수  4  데이터 잔차 RMS 0.0261  참함수와 최대 오차 0.0267
차수  8  데이터 잔차 RMS 0.0240  참함수와 최대 오차 0.0365
차수 20  데이터 잔차 RMS 0.0000  참함수와 최대 오차 148.7384
cond: 단항식 8.3e+08  Chebyshev 8.6e+03

이 결과는 세 가지를 보여 준다.

• 데이터 잔차는 차수와 함께 단조롭게 줄어 20차(보간)에서 0이 된다. 참함수와의 오차는 4차에서 가장 작고, 20차에서는 148로 폭발한다. 데이터에 잘 맞는 것과 참함수에 가까운 것은 다르다.
• 4차 이후 잔차가 거의 줄지 않는다(0.026 → 0.024). 잔차가 잡음 수준에서 멈추는 지점이 적절한 차수의 신호다. 실제로는 교차검증(16장)이나 정보 기준으로 고른다.
• 같은 21개 점에서 단항식 기저의 조건수는 $8 \times 10^8$, Chebyshev 기저는 $9 \times 10^3$이다. np.polyfit 대신 np.polynomial.Chebyshev.fit을 쓰는 이유다. 다만 기저를 바꿔도 21개 점을 지나는 20차 다항식은 하나뿐이므로, 20차의 폭발 자체는 기저로 해결되지 않는다.

연구에서의 의미. 다항식 근사는 신경망 회귀의 가장 단순한 친척이다. 모델 유연성과 잡음 사이의 균형, 학습 데이터 잔차가 일반화 오차를 말해 주지 않는다는 사실, 기저(특징 표현)가 수치적 안정성을 좌우한다는 사실은 그대로 이어진다. 차이는 다항식에서는 이 모든 것을 정확한 식으로 계산할 수 있다는 점이다.

수치적분

왜 필요한가. 대부분의 적분은 원시함수를 식으로 쓸 수 없다. 측정 데이터의 적분(스펙트럼 피크의 넓이, 누적 강수량, 전류를 적분한 전하량)은 애초에 함수식이 없다. 8.3장에서 scipy.integrate.quad를 블랙박스로 썼다면, 이 절은 그 안에서 무슨 일이 일어나는지와 언제 믿을 수 있는지를 다룬다.

정의. 구적 규칙(quadrature rule)은 적분을 함숫값의 가중합으로 근사한다.

노드 $x_i$와 가중치 $w_i$를 고르는 방법에 따라 규칙이 달라진다. 구간을 폭 $h$인 작은 구간으로 나누어 같은 규칙을 반복 적용한 것을 합성(composite) 규칙이라 한다.

• 사다리꼴 규칙(trapezoidal rule): 작은 구간마다 양 끝을 잇는 직선 아래 넓이. 합성 규칙의 오차는 $-\frac{(b-a)h^2}{12} f''(\xi)$, 즉 $O(h^2)$다.
• Simpson 규칙: 이웃한 두 구간마다 세 점을 지나는 2차식 아래 넓이. 가중치는 $\frac{h}{3}(1, 4, 1)$이고, 합성 규칙의 오차는 $-\frac{(b-a)h^4}{180} f^{(4)}(\xi)$, 즉 $O(h^4)$다. 2차식으로 맞췄는데도 대칭성 덕분에 3차 다항식까지 정확하다.
• Gauss 구적(Gauss quadrature): 노드를 등간격으로 고정하지 않고 노드와 가중치를 모두 최적화한다. $n$점 Gauss–Legendre 규칙은 자유도 $2n$개로 $2n-1$차 다항식까지 정확히 적분한다. 노드는 Legendre 다항식 $P_n$의 근이다.

수렴 차수(order of convergence)는 이 오차식의 $h$ 지수다. 오차가 $Ch^p$이면 $h$를 절반으로 할 때 오차는 $2^p$분의 1이 된다. 그래서 두 해상도의 오차 비율로 $p$를 실험적으로 잴 수 있다.

from scipy.integrate import simpson, trapezoid

f, exact = np.exp, np.e - 1.0                          # ∫_0^1 e^x dx
print(" n   사다리꼴    Simpson    Gauss(n점)")
for n in [4, 8, 16]:
    x = np.linspace(0, 1, n + 1)                       # (n+1,) 등간격 n구간
    xg, wg = np.polynomial.legendre.leggauss(n)        # [-1, 1]의 노드·가중치
    gauss = 0.5 * np.sum(wg * f(0.5 * (xg + 1)))       # [0, 1]로 옮김
    e_g = abs(gauss - exact)
    e_g = f"{e_g:.2e}" if e_g > 1e-13 else "< 1e-13"   # 반올림 한계 아래는 판정만
    print(f"{n:2d}  {abs(trapezoid(f(x), x=x) - exact):.2e}  "
          f"{abs(simpson(f(x), x=x) - exact):.2e}  {e_g}")

 n   사다리꼴    Simpson    Gauss(n점)
 4  8.94e-03  3.70e-05  9.33e-10
 8  2.24e-03  2.33e-06  < 1e-13
16  5.59e-04  1.46e-07  < 1e-13

$n$을 두 배로 할 때 사다리꼴 오차는 4.0배, Simpson 오차는 15.9배 줄어 각각 $p = 2$, $p = 4$를 확인한다. Gauss 구적은 점 4개(7차 다항식까지 정확)로 $10^{-9}$, 점 8개로 이미 배정밀도의 반올림 한계에 닿는다(그 아래의 값은 기계마다 달라지므로 판정만 찍었다). 매끄러운(해석적인) 함수에서 Gauss 구적의 오차는 $h$의 거듭제곱이 아니라 점 수에 대해 지수적으로 준다.

차수는 가정 위에 서 있다. 위의 오차식은 $f''$, $f^{(4)}$가 유계라는 가정에서 나온다. 가정이 깨지면 차수도 깨진다. $\int_0^1 \sqrt{x}\,dx$는 $x = 0$에서 도함수가 무한대라 Simpson 규칙도 $O(h^{1.5})$로 느려지고, 불연속인 피적분함수에서는 어떤 고차 규칙도 $O(h)$를 넘지 못한다. 반대로 매끄러운 주기함수를 한 주기에 걸쳐 적분하면 사다리꼴 규칙이 지수적으로 수렴한다. 이론 차수를 믿지 말고, $h$를 바꾸어 관측 차수를 재는 습관이 필요하다.

적응 구적. scipy.integrate.quad는 Fortran 라이브러리 QUADPACK을 감싼 것으로, Gauss–Kronrod 규칙(Gauss 노드에 점을 추가해 두 근사의 차이로 오차를 추정하는 규칙)을 쓰고, 오차 추정이 큰 구간만 반복해서 나눈다. 그래서 특이점이나 급변 구간이 있어도 그곳에 점을 모은다. 반환값의 두 번째 원소는 추정 오차이고, 보장된 오차 한계가 아니다. 피적분함수에 좁은 피크가 있는데 처음 표본점이 피크를 모두 비켜 가면, quad는 틀린 값과 작은 오차 추정을 함께 돌려준다. 예를 들어 폭 0.001인 가우시안 피크 $e^{-((x-7.3)/0.001)^2}$를 $[0, 100]$에서 적분하면 quad는 값 0과 오차 추정 0을 돌려준다(정확값 $1.77 \times 10^{-3}$). 피크 위치를 알면 피크를 감싸는 좁은 구간의 경계를 points=[7.29, 7.3, 7.31]처럼 알려 주어야 정확한 값이 나온다. 피크 중심 하나만 points=[7.3]으로 주면 여전히 틀린다. Gauss 계열 노드는 구간 끝점을 쓰지 않으므로, 피크가 넓은 구간의 끝에 붙어 있으면 여전히 노드가 피크를 비켜 가기 때문이다.

고차원 적분

왜 필요한가. 확률모형의 기댓값, 베이즈 추론의 주변우도, 통계역학의 분배함수는 수십에서 수백만 차원의 적분이다. 1차원 구적 규칙을 각 축에 적용하는 방식(텐서곱 격자)은 이 영역에서 무너진다.

차원의 저주. 축마다 $m$점을 쓰면 $d$차원 격자의 점은 $N = m^d$개다. 점 예산 $N$이 정해져 있으면 축당 점 수는 $m = N^{1/d}$이고, $p$차 규칙의 오차는

이다. $d$가 커지면 지수 $p/d$가 0에 가까워져 점을 아무리 늘려도 오차가 거의 줄지 않는다. 이를 차원의 저주(curse of dimensionality)라 한다.

Monte Carlo 적분. 균등분포에서 뽑은 $N$개의 점 $\mathbf{x}_j$로

를 계산하면, 중심극한정리에 의해 오차의 표준편차는 $\sigma_f/\sqrt{N}$이다($\sigma_f$는 $f(\mathbf{x})$의 표준편차). 수렴은 느리지만($p/d$ 대신 $1/2$) 비율에 $d$가 없다. 또한 같은 표본의 표준편차로 오차를 추정할 수 있다. 같은 점 4096개로 $\int_{[0,1]^d} \cos(x_1 + \cdots + x_d)\,d\mathbf{x}$(정확값을 식으로 알 수 있다)를 계산해 비교한다.

def exact_cos_sum(d):                                  # ∫ cos(x_1+…+x_d) dx
    return (((np.exp(1j) - 1) / 1j) ** d).real

N = 4096
rng = np.random.default_rng(0)
for d in [1, 2, 3, 4, 6, 12]:
    m = round(N ** (1 / d))                            # 축마다 m점 (m^d = N)
    g = (np.arange(m) + 0.5) / m                       # 중점 규칙 노드
    grid = np.stack(np.meshgrid(*[g] * d), axis=-1).reshape(-1, d)   # (N, d)
    ex = exact_cos_sum(d)
    err_mid = abs(np.mean(np.cos(grid.sum(axis=1))) - ex)
    mc = [np.mean(np.cos(rng.random((N, d)).sum(axis=1))) for _ in range(20)]
    err_mc = np.sqrt(np.mean((np.array(mc) - ex) ** 2))   # 20회 반복의 RMS 오차
    print(f"d = {d:2d}  축당 {m:4d}점  중점격자 오차 {err_mid:.1e}  MC 오차(RMS) {err_mc:.1e}")

d =  1  축당 4096점  중점격자 오차 2.1e-09  MC 오차(RMS) 2.4e-03
d =  2  축당   64점  중점격자 오차 1.0e-05  MC 오차(RMS) 4.7e-03
d =  3  축당   16점  중점격자 오차 3.0e-05  MC 오차(RMS) 7.6e-03
d =  4  축당    8점  중점격자 오차 9.2e-04  MC 오차(RMS) 5.3e-03
d =  6  축당    4점  중점격자 오차 1.2e-02  MC 오차(RMS) 4.4e-03
d = 12  축당    2점  중점격자 오차 7.7e-02  MC 오차(RMS) 4.1e-03

격자의 오차는 1차원 $2 \times 10^{-9}$에서 12차원 $8 \times 10^{-2}$로 일곱 자릿수 이상 나빠진다. MC 오차는 모든 차원에서 $10^{-3}$ 수준에 머문다. 이 함수에서는 4차원과 6차원 사이에서 MC가 앞선다. 교차점은 함수의 매끄러움과 규칙의 차수에 따라 달라지지만, 차원이 커지면 결국 MC 계열 방법이 이긴다는 결론은 일반적이다. MC 오차를 한 번의 실행이 아니라 20번 반복의 RMS로 잰 점에 주의한다. 한 번의 MC 결과는 운 좋게 작을 수 있다.

중간 영역에는 두 가지 절충이 있다. 준몬테카를로(quasi-Monte Carlo)는 무작위 점 대신 공간을 고르게 채우는 저불일치(low-discrepancy) 수열(Sobol 수열 등, scipy.stats.qmc)을 써서, 매끄러운 함수에서 오차를 대략 $(\log N)^d/N$으로 줄인다(Caflisch, 1998). sparse grid는 텐서곱 격자에서 기여가 작은 점들을 덜어 내 점 수의 증가를 늦춘다. 수백 차원 이상의 확률 모형에서는 균등 표본 대신 중요도 표집이나 MCMC(13장)를 쓴다.

참고 문헌

• Trefethen, L. N. (2020). Approximation Theory and Approximation Practice, Extended Edition. SIAM. ISBN 978-1-61197-593-2. — Chebyshev 보간, Runge 현상, 무게중심 공식, 보간과 구적의 수렴 이론
• Trefethen, L. N. (2008). Is Gauss quadrature better than Clenshaw–Curtis? SIAM Review, 50(1), 67–87. doi:10.1137/060659831 — Gauss 구적과 Chebyshev 점 기반 구적의 실제 정확도 비교
• Caflisch, R. E. (1998). Monte Carlo and quasi-Monte Carlo methods. Acta Numerica, 7, 1–49. doi:10.1017/S0962492900002798 — Monte Carlo 적분의 수렴과 준몬테카를로 방법