비선형 연립방정식과 다항식의 근

두 곡선의 교점

미지수 $n$개, 방정식 $n$개를 벡터로 씁니다. $\mathbf x = (x_1, \dots, x_n)$, $\mathbf F(\mathbf x) = (f_1(\mathbf x), \dots, f_n(\mathbf x))$이고 $\mathbf F(\mathbf x) = \mathbf 0$을 푸는 것이 목표입니다. 이 장의 기본 예제는

1차원과 결정적으로 다른 점이 하나 있습니다. 근을 가둘 수 없습니다. 1차원에서는 부호가 바뀌는 구간이 근을 보장했지만, 평면에서 $f_1$과 $f_2$의 부호 조합만으로는 두 곡선의 교점을 가두는 간단한 방법이 없습니다. 그래서 다차원의 근 찾기는 거의 모두 뉴턴 방법과 그 변형이고, 좋은 출발점과 전역 수렴 장치가 중요해집니다.

3장의 고정점 반복도 그대로 일반화됩니다. $\mathbf x = \mathbf G(\mathbf x)$ 꼴로 고쳐 쓰고 반복하면, 고정점 근처에서 $\mathbf G$의 야코비안 노름이 1보다 작을 때(더 정확히는 스펙트럼 반지름 $\rho(G'(\mathbf r)) \lt 1$일 때) 선형으로 수렴합니다. 13장의 야코비·가우스–자이델 반복은 이 틀을 선형계에 쓴 것입니다.

다변수 뉴턴 방법과 야코비안

1차원 뉴턴 방법은 $f$를 접선으로 바꿨습니다. 다변수에서는 각 $f_i$를 접평면으로 바꿉니다(CalculusBook 13장: 접평면). 모든 $f_i$의 1차 테일러 전개를 한꺼번에 쓰면 행렬 하나가 나옵니다.

선형 근사가 0이 되는 $\mathbf s$를 찾으면 됩니다. $\mathbf F(\mathbf x_k) + J(\mathbf x_k)\mathbf s = \mathbf 0$, 곧 선형 연립방정식입니다.

function newton_system(F, J, x, tol, maxit)
  for k = 1 .. maxit
    Fx = F(x)
    if norm(Fx) <= tol then return x
    A  = J(x)                     // n x n, 편도함수 n^2개
    s  = lu_solve(A, -Fx)         // 역행렬 대신 LU: 약 (2/3) n^3 번 연산
    x  = x + s
    if norm(s) <= tol * (1 + norm(x)) then return x
  fail("수렴하지 않음")

상수에 $\|J^{-1}\|$이 들어 있다는 점이 1차원의 $1/|f'(r)|$에 해당합니다. 야코비안이 특이 행렬에 가까우면(두 곡선이 거의 접하면) 수렴 반경이 좁고 근 자체의 조건수도 나쁩니다. 1차원 중근의 다차원판입니다.

비용을 세어 봅시다. 한 번의 뉴턴 반복은 $\mathbf F$ 한 번($n$개 함수), 야코비안 한 번($n^2$개 편도함수), LU 분해 약 $\tfrac23 n^3$ 연산입니다. $n = 1000$이면 LU 분해만 수억 연산이라, 큰 문제에서는 야코비안의 희소 구조를 살리고(12장, 13장), 야코비안을 덜 자주 계산하는 방법(다음 절들)을 씁니다.

감쇠 뉴턴: 직선 탐색으로 길 잃지 않기

무엇을 기준으로 "나아졌다"고 할지 정해야 합니다. 자연스러운 잣대는 잔차의 크기, 곧 가치 함수(merit function) $\phi(\mathbf x) = \tfrac12\|\mathbf F(\mathbf x)\|^2$입니다. 근에서 $\phi = 0$이고 그 밖에서는 양수입니다. 다행히 뉴턴 방향은 언제나 $\phi$를 줄이는 방향입니다.

그러니 방향은 그대로 두고 걸음 길이 $t \in (0, 1]$만 고르면 됩니다. $t = 1$부터 시도해서 잔차가 충분히 줄지 않으면 반으로 줄이는 직선 탐색(line search, backtracking)이 표준입니다. "충분히"의 기준이 아르미호 조건으로, 15장 최적화에서 다시 만납니다.

function damped_newton(F, J, x, tol)
  loop
    Fx = F(x);  if norm(Fx) <= tol then return x
    s  = lu_solve(J(x), -Fx)
    t  = 1
    while norm(F(x + t*s)) > (1 - 1e-4 * t) * norm(Fx)   // 아르미호식 충분 감소
      t = t / 2
      if t < 1e-10 then fail("‖F‖의 국소 최솟값에 갇힘?")
    x = x + t * s                // 근 근처에서는 t = 1, 곧 순수 뉴턴

브로이든 방법: 야코비안 없이

1차원 할선 기울기 $b_{k+1}$는 $b_{k+1}(x_{k+1} - x_k) = f(x_{k+1}) - f(x_k)$를 만족했습니다. 다차원판은 $\mathbf s_k = \mathbf x_{k+1} - \mathbf x_k$, $\mathbf y_k = \mathbf F(\mathbf x_{k+1}) - \mathbf F(\mathbf x_k)$에 대한 할선 조건(secant condition)

브로이든 방법은 매 반복 $\mathbf F$를 한 번만 계산합니다. 셔먼–모리슨 공식으로 $B^{-1}$을 직접 갱신하면 선형계 풀기도 $O(n^2)$으로 줄어듭니다. 대가는 수렴 속도입니다. 근 근처에서 $B_k$가 참 야코비안으로 수렴한다는 보장은 없지만, 걸음 방향에서만큼은 충분히 정확해져 초선형으로 수렴한다는 것이 알려져 있습니다(2차는 아니다). 같은 이유로 처음 $B_0$이 나쁘면(예: 단위 행렬) 느려지거나 실패할 수 있어, 보통 $B_0$만은 차분 야코비안으로 만듭니다.

다항식: 호너 방법

괄호를 안쪽부터 묶으면 됩니다. $2x^4 - 3x^3 + x - 5 = (((2x - 3)x + 0)x + 1)x - 5$. 이것이 호너 방법(Horner's method)입니다.

일반적인 다항식에 대해서는 곱셈 $n$번보다 적게 쓰는 방법이 없다는 것이 증명되어 있습니다. 반올림 오차도 작아서, 계산값은 계수를 상대적으로 약 $2n\varepsilon$만큼 흔든 다항식의 정확한 값입니다(후방 안정, 2장). 3장 예제 3.7의 $(x - 1)^7$ 잡음은 호너의 잘못이 아니라 전개된 계수의 크기 $\sum|a_k x^k| = 128$이 값 $|p(x)|$보다 엄청나게 크기 때문이었습니다.

function horner2(a, x)            // a[n], ..., a[0]: 높은 차수부터
  p = a[n]; d = 0
  for k = n-1 down to 0
    d = d * x + p                // q의 호너 → p'(x)
    p = p * x + a[k]             // p의 호너
  return (p, d)                  // 곱셈 2n번으로 p와 p'

수축: 근을 하나씩 떼어 내기와 그 오차

근 $r$을 하나 찾으면 조립제법의 몫 $q(x) = p(x)/(x - r)$은 차수가 하나 낮고 나머지 근을 모두 가집니다. 이제 $q$의 근을 찾고, 또 나누고… 이렇게 차수를 줄여 가는 것을 수축(deflation)이라 합니다. 문제는 찾은 $r$이 정확하지 않다는 데 있습니다. $r$의 작은 오차가 $q$의 계수를 흔들고, 그 흔들림이 다음 근들로 옮아갑니다.

윌킨슨(James H. Wilkinson)의 분석에 따르면 높은 차수부터 계수를 만드는 앞쪽 수축(조립제법)은 크기가 작은 근부터 떼어 낼 때 안정하고, 큰 근부터 떼면 작은 근들이 크게 망가질 수 있습니다. 실무에서는 두 가지를 더 합니다. ① 수축된 다항식에서 찾은 근을 원래 다항식에 대고 뉴턴 몇 번으로 다듬는다(polishing). ② 가능하면 수축 대신 아래의 동반 행렬처럼 모든 근을 한꺼번에 구한다.

동반 행렬: 모든 근을 고유값으로

방향을 거꾸로 생각해 봅시다. 행렬의 고유값은 특성 다항식의 근입니다. 그렇다면 주어진 다항식을 특성 다항식으로 갖는 행렬을 만들고, 그 고유값을 안정한 행렬 알고리즘으로 구하면 다항식의 근이 모두 나옵니다. 14장에서 보듯 고유값을 특성 다항식으로 구하는 것은 나쁜 생각이지만, 이 반대 방향은 좋은 생각입니다.

MATLAB의 roots와 NumPy의 numpy.roots가 바로 이렇게 동반 행렬의 고유값을 계산합니다. 고유값은 이동 QR 알고리즘으로 구합니다(14장). $H - \mu I = QR$로 분해하고 $H \leftarrow RQ + \mu I$로 바꾸는 닮음 변환을 되풀이하면 아래 대각 성분이 0으로 가고 대각에 고유값이 드러납니다. 동반 행렬은 이미 헤센베르크 꼴이라 한 단계가 $O(n^2)$이고 전체 비용은 $O(n^3)$, 대신 출발점 걱정도 수축 순서 걱정도 없습니다.

복소 평면의 뉴턴 방법과 프랙털

뉴턴 방법은 복소수에서도 똑같이 쓸 수 있습니다. $z_{k+1} = z_k - p(z_k)/p'(z_k)$. 실근이 없는 $z^2 + 1$도 복소 출발점에서 $\pm i$로 수렴합니다. 사실 이것은 2변수 뉴턴과 같은 것입니다. $z = x + iy$, $p = u + iv$로 쓰면 $p(z) = 0$은 $u(x, y) = 0$, $v(x, y) = 0$의 연립방정식이고, 코시–리만 방정식 $u_x = v_y$, $u_y = -v_x$ 때문에 그 야코비안 $\begin{bmatrix} u_x & u_y \\ v_x & v_y\end{bmatrix}$은 복소수 $p'(z)$를 곱하는 회전·확대 행렬이 됩니다(MathBook: 복소수). 2변수 뉴턴 시뮬레이터의 둘째 연립방정식이 바로 $z^3 = 1$의 실수·허수 부분입니다.

$z^3 - 1$의 세 근 $1$, $e^{\pm 2\pi i/3}$의 끌개 영역을 색칠하면 놀라운 그림이 나옵니다. 근 가까이에서는 각 영역이 넓고 매끄럽지만 경계에서는 세 색이 끝없이 뒤섞입니다. 이 질문은 케일리(Arthur Cayley)가 1879년에 2차식은 쉽고 3차식은 어렵다고 적으며 던진 것으로 알려져 있고, 컴퓨터 그래픽이 생긴 뒤에야 그 모양이 보였습니다. 이유는 복소 동역학의 정리가 말해 줍니다. 유리 함수 반복에서 각 끌개 영역의 경계는 모두 같은 집합(쥘리아 집합)이라서, 경계의 어느 점이든 그 아무리 작은 근방에도 세 영역이 모두 들어 있습니다.

• 실수축 위의 출발점은 실수축을 떠나지 않으므로 실근 1로 간다(0과 0으로 가는 점들 제외). $z_0 = -1$이면 $-0.3333$, $2.7778$, $1.8950$, $1.3562$, $1.0854$, $1.0065$, $1.00004$.
• $z_0 = 1 + i$는 $0.6667 + 0.5i$, $0.5788 - 0.1275i$, $1.2470 + 0.3136i$, $1.0090 + 0.1137i$, $\dots$로 돌아서 1로 가고, $z_0 = 0.1 + 0.01i$는 0 가까이라 $p'$이 작아 첫 걸음에 $32.4 - 6.53i$로 크게 튄 뒤 13번 만에 1로 돌아온다.
• $z^3 - 2z + 2$에서는 3장의 순환 $0 \to 1 \to 0$이 끌어당기는 주기 궤도라서, 그 둘레에 어떤 근으로도 수렴하지 않는 출발점의 열린 영역이 생긴다(시뮬레이터의 회색 영역).

연속법: 쉬운 문제에서 어려운 문제로

매개변수 $t \in [0, 1]$로 두 문제를 잇는 호모토피(homotopy) $\mathbf H(\mathbf x, t)$를 만듭니다. $t = 0$에서는 해를 알고, $t = 1$에서는 원래 문제 $\mathbf F$가 됩니다. 대표적인 두 가지는

해 $\mathbf x(t)$는 $\mathbf H(\mathbf x(t), t) = \mathbf 0$을 만족하는 곡선(해 경로)을 그립니다. 연속법(continuation method)은 이 경로를 작은 걸음으로 따라갑니다. 각 걸음은 예측(앞 해, 또는 접선 방향으로 외삽)과 수정(그 점에서 $\mathbf H(\cdot, t_{k+1}) = 0$에 대한 뉴턴 몇 번)입니다. 뉴턴 호모토피를 $t$로 미분하면 $J\dot{\mathbf x} + \mathbf F(\mathbf x_0) = 0$, 곧 $\dot{\mathbf x} = -J^{-1}\mathbf F(\mathbf x_0)$입니다. 경로를 따라 잔차가 $\mathbf F(\mathbf x(t)) = (1 - t)\mathbf F(\mathbf x_0)$로 방향을 바꾸지 않고 줄어드니, 연속법은 뉴턴 방향을 아주 작은 걸음으로 따라가는 "연속 뉴턴"인 셈입니다.

function continuation(H, Hx, x0, N)
  x = x0                           // H(x0, 0) = 0
  for k = 1 .. N
    t = k / N
    repeat 3~5 times                 // 수정자: t를 고정한 뉴턴
      x = x - lu_solve(Hx(x, t), H(x, t))
  return x                         // t = 1: F(x) = 0의 해

$t$를 단조롭게 늘리는 이 단순한 방법은 경로가 되돌아갈 때($t$가 최댓값을 찍고 줄어드는 전환점) 실패합니다. 전환점에서는 $\partial\mathbf H/\partial\mathbf x = J$가 특이 행렬이 되어 $t$를 더 늘린 문제에 가까운 해가 없습니다. 해결책은 $t$도 미지수로 보고 경로의 호의 길이를 매개변수로 삼는 의사 호 길이 연속법입니다. 아래 시뮬레이터가 이 방법을 씁니다.

응용: 회로의 동작점, 로봇 팔, GPS

회로의 동작점: 지수 함수와 걸음 제한

5 V 전원 — 저항 $R_1 = 1$ kΩ — 마디 1 — 다이오드 $D_1$ — 접지, 그리고 마디 1 — $R_2 = 1$ kΩ — 마디 2 — $D_2$ — 접지인 회로를 생각합니다. 다이오드 전류는 $I_D = I_s(e^{V/V_T} - 1)$($I_s = 10^{-14}$ A, $V_T = 25.852$ mV)이고, 마디마다 들어오고 나가는 전류의 합이 0입니다(KCL).

로봇 팔의 역기구학

두 마디 팔(길이 $L_1, L_2$, 관절 각 $\theta_1, \theta_2$)의 끝 위치는 $x = L_1\cos\theta_1 + L_2\cos(\theta_1 + \theta_2)$, $y = L_1\sin\theta_1 + L_2\sin(\theta_1 + \theta_2)$입니다(정기구학). 원하는 위치 $(p_x, p_y)$를 주고 각을 구하는 역기구학은 $\mathbf F(\theta_1, \theta_2) = (x - p_x, y - p_y) = \mathbf 0$이고, 야코비안의 행렬식은 $\det J = L_1L_2\sin\theta_2$입니다. 팔이 완전히 펴지거나 접히면($\theta_2 = 0, \pi$) 특이 자세라 뉴턴이 불안정해집니다.

$L_1 = L_2 = 1$, 목표 $(1, 1)$, 출발 $(\theta_1, \theta_2) = (0.3, 1.0)$이면 뉴턴 반복은 $(-0.0675, 1.7115)$, $(-0.00243, 1.5752)$, $(-1.96\times10^{-6}, 1.570801)$, $(-4.3\times10^{-12}, 1.5707963268)$으로 4번 만에 $(0, \pi/2)$에 닿습니다. 목표가 팔이 닿지 않는 곳이면 해가 없고 뉴턴은 진동합니다. 이때는 $\|\mathbf F\|^2$를 최소화하는 감쇠 최소제곱 걸음 $(J^TJ + \lambda I)\mathbf s = -J^T\mathbf F$(레벤버그–마쿼트, 15장)가 팔을 목표 쪽으로 쭉 뻗게 해 줍니다.

GPS: 의사거리 방정식

위성 $i$의 위치 $\mathbf s_i$를 알고 신호의 도착 시각으로 거리 $\rho_i$를 재는데, 수신기의 값싼 시계에는 오차 $b$(거리로 환산)가 있습니다. 미지수는 위치 $\mathbf x$와 $b$이고, $f_i(\mathbf x, b) = \|\mathbf x - \mathbf s_i\| + b - \rho_i = 0$입니다. 3차원이면 미지수 4개라 위성이 4개 이상 필요합니다. 야코비안의 $i$번째 행은 $(\text{위성에서 수신기로의 단위 벡터}, 1)$입니다. 위성이 더 많으면 최소제곱(가우스–뉴턴, 7장·15장)으로 풉니다.

이 도구가 쓰이는 곳