최소제곱 근사와 직교 다항식

잡음 데이터와 보간의 실패

보간과 근사의 차이는 미지수와 조건의 개수에 있습니다. $m$개의 점 $(x_i, y_i)$를 $n+1$개의 계수를 가진 모형 $p(x) = c_0\phi_0(x) + \dots + c_n\phi_n(x)$로 설명하려 할 때, $m = n + 1$이면 보간(정사각 선형계), $m > n + 1$이면 식이 미지수보다 많은 과결정계(overdetermined system)입니다.

잡음이 있으면 모든 식을 만족하는 $\mathbf c$는 보통 없습니다. 대신 각 식의 어긋남, 곧 잔차(residual) $r_i = y_i - p(x_i)$를 모두 작게 만드는 $\mathbf c$를 고릅니다. "작게"를 어떻게 재느냐가 선택인데, 제곱합을 고르면 계산이 선형대수로 깔끔하게 떨어집니다.

왜 하필 제곱일까요? 세 가지 이유가 있습니다. 첫째, 제곱합은 $\mathbf c$의 매끄러운 이차 함수라 미분해서 0으로 놓으면 선형 방정식이 나옵니다. 둘째, 기하적으로 유클리드 거리이므로 답이 직교 사영입니다. 셋째, 잡음이 서로 독립인 정규 분포를 따르면 최소제곱 해는 최대 가능도 추정량이 되고, 가우스–마르코프 정리에 따르면 잡음의 평균이 0이고 분산이 같으며 서로 무상관이기만 해도 선형 불편 추정량 가운데 분산이 가장 작습니다. 대신 제곱은 큰 잔차를 매우 무겁게 벌하므로 이상치 하나에 결과가 끌려갑니다. 이 약점은 5절에서 다룹니다.

정규 방정식과 사영

목적 함수 $F(\mathbf c) = \|A\mathbf c - \mathbf y\|^2 = \mathbf c^\top A^\top A\mathbf c - 2\mathbf c^\top A^\top\mathbf y + \mathbf y^\top\mathbf y$의 기울기를 0으로 놓으면 $\nabla F = 2A^\top A\mathbf c - 2A^\top\mathbf y = \mathbf 0$입니다.

이 정리는 그림 한 장으로 요약됩니다. $A\mathbf c$가 만들 수 있는 벡터 전체는 $A$의 열공간 $\mathcal C(A)$이고, $\mathbb R^m$ 안의 $n$차원 평면입니다. $\mathbf y$는 보통 이 평면 밖에 있습니다. 평면 위에서 $\mathbf y$에 가장 가까운 점은 $\mathbf y$에서 평면에 내린 수선의 발 $\mathbf p = A\hat{\mathbf c}$이고, 잔차는 그 수선입니다(LinearBook 9장: 직교성과 사영).

함수 정규방정식_최소제곱(A, y)        // A: m×n, 열이 일차독립
  G ← AᵀA                               // n×n 대칭 양정치, 비용 약 m n²
  g ← Aᵀy                               // 비용 m n
  L ← 촐레스키(G)                        // G = L Lᵀ, 비용 약 n³/3
  z ← 전진대입(L, g)                     // L z = g
  c ← 후진대입(Lᵀ, z)                    // Lᵀ c = z
  반환 c

아래 시뮬레이터의 왼쪽은 데이터 평면, 오른쪽은 $m = 3$일 때의 $\mathbb R^3$입니다. 데이터 세 점의 높이 $y_1, y_2, y_3$이 $\mathbb R^3$의 한 점 $\mathbf y$가 되고, 직선 두 계수가 만들 수 있는 예측값 $(c_0 + c_1x_1, c_0 + c_1x_2, c_0 + c_1x_3)$은 평면을 이룹니다. 왼쪽 그림에서 잔차 제곱의 넓이 합을 줄이는 일과 오른쪽 그림에서 수선을 내리는 일이 같은 일임을 확인해 보세요.

정규 방정식의 조건수 제곱 문제와 QR

2장에서 본 대로 선형계 $G\mathbf x = \mathbf g$를 배정밀도로 풀면 상대 오차가 대략 $\kappa(G)\,\varepsilon_{\text{mach}}$까지 커질 수 있습니다($\varepsilon_{\text{mach}} \approx 2.2\times10^{-16}$, 2장, 12장). 직사각 행렬의 조건수는 특이값으로 정의합니다.

최소제곱 문제 자체의 민감도는 대략 $\kappa(A)$(잔차가 크면 $\kappa(A)^2$항이 붙지만 그 계수는 잔차 크기에 비례)인데, 정규 방정식이라는 알고리즘은 잔차와 상관없이 $\kappa(A)^2$의 선형계를 풉니다. $\kappa(A) = 10^8$이면 $\kappa(A^\top A) = 10^{16}$이고, 배정밀도의 16자리가 모두 사라집니다. 2장의 용어로 정규 방정식은 후방 안정하지 않은 알고리즘입니다.

처방은 $A^\top A$를 만들지 않는 것입니다. $A = QR$($Q$는 열이 정규직교인 $m\times n$, $R$은 $n\times n$ 상삼각)로 분해하면 직교 행렬은 길이를 바꾸지 않으므로

함수 QR_최소제곱(A, y)                 // 하우스홀더, A는 m×n
  for k = 1 … n
    v ← A[k:m, k]                        // k번째 열의 아래쪽
    v[1] ← v[1] + sign(v[1])·‖v‖         // 상쇄를 피하는 부호 선택
    A[k:m, k:n] ← A[k:m, k:n] − 2 v (vᵀ A[k:m, k:n]) / (vᵀv)
    y[k:m]      ← y[k:m]      − 2 v (vᵀ y[k:m]) / (vᵀv)   // Qᵀy를 같이 만든다
  end
  c ← 후진대입(A[1:n, 1:n], y[1:n])      // R c = (Qᵀy)의 위쪽 n개
  반환 c                                // 비용 약 2mn² − 2n³/3

비용은 $m \gg n$에서 정규 방정식(약 $mn^2$)의 두 배 정도입니다. 그 값으로 자릿수를 $\log_{10}\kappa(A)$개 더 지킵니다. 더 나쁜 경우(열이 거의 일차종속)에는 특이값 분해로 작은 특이값을 잘라 내는 방법이 가장 튼튼하지만 비용이 몇 배 더 듭니다(LinearBook: 그람–슈미트와 QR). 아래 표는 $[0, 1]$의 등간격 50개 점에서 계수가 모두 1인 $n$차 다항식을 잡음 없이 다시 맞춘 결과입니다.

계수 오차는 최대 절댓값 오차입니다. 정규 방정식의 오차가 대략 $\kappa(A)^2\varepsilon_{\text{mach}}$를, QR의 오차가 대략 $\kappa(A)\varepsilon_{\text{mach}}$를 따라가는 것이 보입니다. 단항식 기저 $1, x, x^2, \dots$ 자체도 범인입니다. $[0, 1]$에서 $x^9$과 $x^{10}$은 거의 같은 모양이라 열이 거의 평행합니다. 기저를 $[-1, 1]$로 옮기거나 체비쇼프 다항식으로 바꾸면 $\kappa(A)$가 극적으로 줄어듭니다. 6절에서 이 생각을 끝까지 밀고 갑니다.

다항식 차수와 과적합, 검증 오차

문제는 우리가 원하는 것이 측정한 점에서의 오차가 아니라 측정하지 않은 점에서의 오차라는 데 있습니다. 측정한 점(학습 데이터)에서 잰 오차를 학습 오차(training error), 맞춤에 쓰지 않은 점에서 잰 오차를 검증 오차(validation error)라 합니다. 차수를 올리면 학습 오차는 계속 줄지만, 어느 순간부터 모형이 추세 대신 잡음을 외우기 시작해 검증 오차가 커집니다. 이것이 과적합(overfitting)입니다.

잡음의 표준편차가 $\sigma$이면 아무리 좋은 모형도 새 측정값과의 RMS 오차를 $\sigma$보다 작게 만들 수 없습니다. 학습 오차가 $\sigma$보다 훨씬 작아졌다는 것은 잡음을 외웠다는 신호입니다. 반대로 차수가 너무 낮으면 학습·검증 오차가 모두 큽니다(과소적합). 오차를 "모형이 참 함수를 표현하지 못해서 생기는 부분(편향)"과 "잡음에 따라 맞춘 곡선이 흔들리는 부분(분산)"으로 나누면, 차수를 올릴수록 편향은 줄고 분산은 커집니다.

데이터가 적어 검증용으로 떼어 둘 여유가 없으면 교차 검증(cross-validation)을 씁니다. 데이터를 $k$묶음으로 나눠 하나씩 돌아가며 검증에 쓰고 평균을 냅니다. 선형 최소제곱에서는 한 점씩 뺀 검증(LOOCV) 오차를 다시 맞추지 않고도 $r_i/(1 - h_{ii})$로 바로 계산할 수 있습니다. 여기서 $h_{ii}$는 사영 행렬 $H = A(A^\top A)^{-1}A^\top$의 대각 원소입니다. 또 하나의 처방은 계수의 크기에 벌점을 주는 정칙화(regularization)로, $\|A\mathbf c - \mathbf y\|^2 + \lambda\|\mathbf c\|^2$을 최소화하는 능형 회귀(티호노프 정칙화)가 대표적입니다.

가중 최소제곱과 강건 회귀

먼저 측정마다 신뢰도가 다를 때입니다. $i$번째 측정의 잡음 표준편차가 $\sigma_i$로 알려져 있으면, 정밀한 측정의 잔차를 더 무겁게 쳐야 합니다.

이상치는 어느 점의 $\sigma_i$가 큰지 모르는 경우입니다. 그래서 잔차를 보고 가중치를 스스로 정합니다. 제곱 대신 큰 잔차에서 덜 가파르게 자라는 손실 $\rho(r)$을 쓰는데, 대표가 후버 손실(Huber loss)입니다.

$\sum\rho(r_i)$을 최소화하는 조건 $\sum\rho'(r_i)\,\partial r_i/\partial\mathbf c = 0$은 $\sum w(r_i)\,r_i\,\partial r_i/\partial\mathbf c = 0$으로 쓸 수 있고, 이는 가중치 $w(r_i)$를 고정했을 때의 가중 최소제곱 조건과 같습니다. 그래서 가중치를 고정하고 풀고, 새 잔차로 가중치를 다시 계산하는 반복이 자연스럽습니다.

함수 IRLS(A, y, δ)                     // 반복 재가중 최소제곱
  c ← QR_최소제곱(A, y)                  // 보통 최소제곱으로 출발
  repeat
    r ← y − A c
    w_i ← 1        (|r_i| ≤ δ)
    w_i ← δ/|r_i|  (|r_i| > δ)            // 큰 잔차의 영향을 깎는다
    c ← QR_최소제곱(√W A, √W y)
  until c가 거의 변하지 않음
  반환 c

후버 손실은 볼록이라 IRLS가 하나의 답으로 수렴하지만, 이상치가 아주 크고 많으면 여전히 영향을 받습니다. 영향을 완전히 끊는 손실(투키의 이중가중 등)은 볼록이 아니어서 출발점에 따라 답이 달라집니다. 이상치가 데이터의 절반 가까이면 RANSAC처럼 작은 부분집합을 무작위로 골라 맞추고 동의하는 점을 세는 방법을 씁니다. 컴퓨터 비전에서 두 사진의 대응점으로 변환을 추정할 때가 대표적입니다.

연속 최소제곱과 직교 다항식

합을 적분으로 바꾸면 됩니다. 구간 $[a, b]$와 양의 가중 함수 $w(x)$에 대해 내적 $\langle f, g\rangle = \int_a^b f(x)g(x)w(x)\,dx$와 노름 $\|f\| = \sqrt{\langle f, f\rangle}$을 정의하면, 함수도 벡터처럼 다룰 수 있습니다. $n$차 이하 다항식 전체 $\mathcal P_n$이 "열공간" 역할을 합니다.

단항식 $\phi_k = x^k$를 $[0, 1]$, $w = 1$로 쓰면 $\langle x^j, x^k\rangle = 1/(j + k + 1)$이고, 정규 방정식의 행렬은 유명한 힐베르트 행렬(Hilbert matrix)입니다. 그 조건수는 $5\times5$에서 $4.8\times10^5$, $10\times10$에서 $1.6\times10^{13}$으로, 10차 근사만 해도 배정밀도의 자릿수 대부분이 사라집니다. 3절과 같은 병입니다. 처방도 같습니다. 기저를 직교화합니다.

그람–슈미트 과정(LinearBook: 그람–슈미트와 QR)을 $1, x, x^2, \dots$에 $[-1, 1]$, $w = 1$의 내적으로 적용하면 차례로

이 다항식들을 $P_k(1) = 1$이 되도록 크기를 맞춘 것이 르장드르 다항식(Legendre polynomials) $P_0 = 1$, $P_1 = x$, $P_2 = \frac12(3x^2 - 1)$, $P_3 = \frac12(5x^3 - 3x)$입니다. 매번 앞의 모든 다항식에 대한 사영을 뺄 필요는 없습니다. 직교 다항식은 항상 세 항 점화식을 만족합니다.

그래서 $n$차 직교 다항식의 값은 $O(n)$번의 연산으로 안정하게 계산되고, 근사 계수는 적분 하나씩으로 따로따로 구해집니다. 차수를 하나 올려도 앞의 계수는 바뀌지 않습니다. 단항식 기저에서는 차수를 올리면 모든 계수를 다시 풀어야 했습니다.

체비쇼프 전개와 계수 감소

가중 함수를 $w(x) = 1/\sqrt{1 - x^2}$로 바꾸면 또 하나의 직교 다항식열이 나옵니다. $x = \cos\theta$로 치환하면 이 가중치는 $dx/\sqrt{1 - x^2} = -d\theta$가 되어 사라지고, 직교 다항식은 놀랍도록 단순해집니다.

점화식은 삼각함수의 합 공식 $\cos(k + 1)\theta + \cos(k - 1)\theta = 2\cos\theta\cos k\theta$ 그대로이고, 직교성은 $\int_0^\pi\cos j\theta\cos k\theta\,d\theta = 0$ 그대로입니다. 곧 체비쇼프 전개 $f(x) = \sum a_kT_k(x)$는 $f(\cos\theta)$의 코사인 푸리에 급수입니다. 그래서 계수는 체비쇼프 점 $x_j = \cos\frac{(j + 1/2)\pi}{N}$에서의 함숫값에 이산 코사인 변환을 하면 $O(N\log N)$에 얻어집니다(8장의 FFT). $T_k$는 $[-1, 1]$에서 $|T_k| \le 1$이므로 꼬리를 잘라 낸 오차는 버린 계수의 합으로 바로 묶입니다.

계수를 로그 눈금으로 그리면 함수의 성격이 그대로 보입니다. $e^x$는 어디서나 해석적이라 계수가 $1.27, 1.13, 0.27, 0.044, 0.0055, 5.4\times10^{-4}, 4.5\times10^{-5}, 3.2\times10^{-6}$으로 매 차수 열 배 가까이 줄고, 15차 근처에서 기계 엡실론에 닿습니다. $1/(1 + 25x^2)$은 $x = \pm i/5$에 극이 있어 타원이 좁고($\rho \approx 1.22$) 감소가 느립니다. 5장의 룽게 함수가 체비쇼프 노드에서도 꽤 많은 점을 요구한 이유입니다. $|x|$는 꺾인 점 때문에 $k^{-2}$, 계단 함수는 $k^{-1}$로 줄어듭니다. 이 계수의 꼴은 다음 장에서 푸리에 계수의 감소로 다시 만납니다.

잘린 체비쇼프 급수의 오차는 대체로 첫 번째 버린 항 $a_{n+1}T_{n+1}(x)$이 지배하고, $T_{n+1}$은 $[-1, 1]$에서 $+1$과 $-1$을 $n + 2$번 번갈아 찍는 함수입니다. 그래서 잘린 체비쇼프 급수의 오차 곡선은 거의 고른 진동을 합니다. $e^x$의 3차 체비쇼프 절단은 최대 오차 $6.07\times10^{-3}$으로, 첫 버린 계수 $a_4 = 5.47\times10^{-3}$과 거의 같습니다. 다음 절에서 이 "거의"를 "정확히"로 바꿉니다. 소프트웨어 Chebfun은 이 원리로 함수를 체비쇼프 계수가 기계 엡실론 아래로 떨어질 때까지 자동으로 전개해 숫자처럼 다룹니다.

최소최대 근사와 등진동: 레메즈 알고리즘

$\|f - p\|_\infty$을 최소화하는 $p^* \in \mathcal P_n$을 최소최대 근사(minimax approximation) 또는 최량 균등 근사라 합니다. 연속 함수에 대해 존재하고 유일하며, 놀랍게도 오차 곡선의 모양만 보고 판정할 수 있습니다.

등진동 정리에서 체비쇼프 다항식의 또 다른 얼굴이 나옵니다. 최고차 계수가 1인 $n$차 다항식 가운데 $[-1, 1]$에서 최댓값이 가장 작은 것은 $2^{1-n}T_n$이고 그 최댓값은 $2^{1-n}$입니다($x^n$을 $n - 1$차 이하로 최소최대 근사한 오차가 $2^{1-n}T_n$이고, 이것이 $n + 1$개 점에서 등진동하기 때문). 5장의 체비쇼프 노드가 보간 오차의 $\prod(x - x_i)$를 최소로 만든 이유가 바로 이것입니다(5장).

일반적인 경우에는 등진동점을 미리 알 수 없으므로 추측하고 고치기를 되풀이합니다. 이것이 레메즈 알고리즘(Remez exchange algorithm)입니다.

함수 레메즈(f, n, [a, b])
  x₀ < … < x_{n+1} ← 처음 기준점      // 보통 체비쇼프 극점 cos(iπ/(n+1))
  repeat
    // n+2개 미지수 (c₀…cₙ, E)의 선형계: 기준점에서 정확히 등진동하게
    c₀ + c₁xᵢ + … + cₙxᵢⁿ + (−1)ⁱ E = f(xᵢ),   i = 0 … n+1
    e(x) ← f(x) − p(x)                     // 촘촘한 격자에서 계산
    새 기준점 ← e의 극값 중 부호가 번갈아 바뀌는 n+2개 (전역 최대 포함)
  until max|e| ≈ |E|                        // 기준점에서의 높이 = 전체 최댓값
  반환 p, |E|

매 단계에서 $|E|$는 최소최대 오차의 아래 한계(드 라 발레 푸생 정리), $\max|e|$는 위 한계입니다. 두 값이 만나면 끝입니다. 매끄러운 함수에서는 이 간격이 2차로 줄어 몇 번 만에 수렴합니다. $e^x$의 $[-1, 1]$ 3차 근사는 체비쇼프 극점에서 출발하면 첫 단계에 $|E| = 5.474\times10^{-3}$, $\max|e| = 5.581\times10^{-3}$, 둘째 단계에 두 값이 $5.528\times10^{-3}$에서 넷째 자리까지 일치하고, 셋째 단계에 배정밀도 수준으로 만납니다.

표가 실무의 교훈을 줍니다. 체비쇼프 절단이나 체비쇼프 보간은 계산이 훨씬 쉬운데도 최소최대보다 10~20%만 나쁩니다. 그래서 대부분의 응용에서는 체비쇼프 근사로 충분하고, 수학 라이브러리처럼 마지막 한 비트가 중요한 곳에서만 레메즈로 다듬습니다. 실제 라이브러리는 먼저 범위 축소(예: $e^x = 2^k e^r$, $|r| \le \frac12\ln 2$)로 구간을 좁힌 뒤 그 작은 구간에서 낮은 차수의 최소최대 다항식(또는 유리 함수)을 호너 방법으로 계산합니다(4장).

이 도구가 쓰이는 곳