# 다변수 시스템에서의 비선형 근사(뉴턴 방법 확장)

#### 다변수 뉴턴 방법의 기본 아이디어

다변수 뉴턴 방법은 다차원 공간에서 정의된 비선형 함수의 해를 근사하기 위한 강력한 알고리즘이다. 일반적으로 하나의 독립변수와 하나의 종속변수가 있는 스칼라 방정식 $f(x)=0$을 풀기 위해 사용하는 뉴턴 방법을 여러 변수로 확장한 형태이다. 비선형 시스템이 다음과 같은 형태라고 하자.

$$
\begin{align} \mathbf{f}(\mathbf{x}) =  \begin{pmatrix} f\_1(\mathbf{x}) \ f\_2(\mathbf{x}) \ \vdots \ f\_n(\mathbf{x}) \end{pmatrix} = \mathbf{0} \end{align}
$$

여기에서 $\mathbf{x} \in \mathbb{R}^n$이고 $\mathbf{f}(\mathbf{x}) \in \mathbb{R}^n$이다. 스칼라 뉴턴 방법은 단순히 도함수를 역수로 곱하는 방식으로 근사를 진행한다. 다변수 뉴턴 방법은 그 대신 야코비(Jacobian) 행렬을 사용하여 방향 벡터를 구한다.

#### 야코비 행렬과 잔차 벡터

야코비 행렬 $J(\mathbf{x})$는 $\mathbf{x}=(x\_1,x\_2,\dots,x\_n)$에서 각 성분함수 $f\_i(\mathbf{x})$에 대한 편미분으로 구성되는 $n\times n$ 행렬이다. 즉

$$
\begin{align} J(\mathbf{x}) = \frac{\partial \mathbf{f}}{\partial \mathbf{x}} = \begin{pmatrix} \frac{\partial f\_1}{\partial x\_1}(\mathbf{x}) & \frac{\partial f\_1}{\partial x\_2}(\mathbf{x}) & \cdots & \frac{\partial f\_1}{\partial x\_n}(\mathbf{x})  \
$$

6pt] \frac{\partial f\_2}{\partial x\_1}(\mathbf{x}) & \frac{\partial f\_2}{\partial x\_2}(\mathbf{x}) & \cdots & \frac{\partial f\_2}{\partial x\_n}(\mathbf{x}) \\

$$
6pt] \vdots & \vdots & \ddots & \vdots  \
$$

6pt] \frac{\partial f\_n}{\partial x\_1}(\mathbf{x}) & \frac{\partial f\_n}{\partial x\_2}(\mathbf{x}) & \cdots & \frac{\partial f\_n}{\partial x\_n}(\mathbf{x}) \end{pmatrix}. \end{align}

$$
$\mathbf{x}\_k$에서의 잔차 벡터(residual vector)는 $\mathbf{f}(\mathbf{x}\_k)$로 나타낼 수 있으며, 이는 차후 수렴 판단에도 중요한 지표로 쓰인다.

### 뉴턴-라프슨(Newton-Raphson) 공식의 다변수 확장

스칼라 뉴턴 방법에서 한 번의 반복은 $x\_{k+1} = x\_k - \frac{f(x\_k)}{f'(x\_k)}$로 정의된다. 이 개념을 다변수로 확장하려면, 편도함수의 역수를 야코비 행렬의 역행렬로 바꾸면 된다. 즉 반복 공식은 다음과 같다.
$$

\begin{align} \mathbf{x}\_{k+1} = \mathbf{x}\_k - \[J(\mathbf{x}\_k)]^{-1} \mathbf{f}(\mathbf{x}\_k). \end{align}

$$
야코비 행렬이 정칙(가역)이어야만 다음 반복 단계를 진행할 수 있다. 수치 구현에서는 매 반복마다 $J(\mathbf{x}\_k)$를 계산하고 이를 통해 $\Delta \mathbf{x}\_k = \[J(\mathbf{x}\_k)]^{-1} \mathbf{f}(\mathbf{x}*k)$를 구한 뒤, $\mathbf{x}*{k+1} = \mathbf{x}\_k - \Delta \mathbf{x}\_k$로 갱신한다.

### 선형화와 수렴 특징

다변수 뉴턴 방법에서, $\mathbf{f}(\mathbf{x})$를 $\mathbf{x}\_k$ 부근에서 선형화하면
$$

\begin{align} \mathbf{f}(\mathbf{x}\_k + \Delta \mathbf{x}) \approx \mathbf{f}(\mathbf{x}\_k) + J(\mathbf{x}\_k),\Delta \mathbf{x}. \end{align}

$$
이 근사를 $\mathbf{0}$과 동일시하여 $\mathbf{f}(\mathbf{x}\_k + \Delta \mathbf{x}) \approx \mathbf{0}$이라고 두면
$$

\begin{align} J(\mathbf{x}\_k),\Delta \mathbf{x} \approx -\mathbf{f}(\mathbf{x}\_k). \end{align}

$$
따라서 $\Delta \mathbf{x} = - \[J(\mathbf{x}\_k)]^{-1} \mathbf{f}(\mathbf{x}\_k)$가 되고, 이로부터 뉴턴 방정식이 바로 도출된다. 적절한 초기값 $\mathbf{x}\_0$ 근방에서 야코비 행렬이 적당히 잘 정의되어 있고 역행렬이 안정적으로 계산된다면, 다변수 뉴턴 방법은 국소적으로 2차 수렴(quadratic convergence) 특성을 가진다.

### 뉴턴 방법의 알고리즘적 흐름

mermaid를 통해 반복 흐름을 간단히 표현하면 아래와 같다.

```mermaid
flowchart TB
    A((Start)) --> B[Initialize<br>x₀]
    B --> C["Compute<br>J(x₀) and f(x₀)"]
    C --> D["Solve<br>J(x₀)*Δx = -f(x₀)"]
    D --> E[Update x₁ = x₀ + Δx]
    E --> F{Check<br>convergence?}
    F -- Yes --> G((End))
    F -- No --> H[Set x₀ = x₁ and repeat]
```

이 과정에서는 매 반복마다 야코비 계산 및 역행렬(또는 이에 준하는 연산)을 수행해야 하므로 연산 비용이 상당할 수 있다. 따라서 매우 높은 차원의 경우, $\[J(\mathbf{x}\_k)]^{-1}$를 직접 구하는 대신 선형대수 알고리즘(가우스 소거법, LU 분해, QR 분해, 등)을 이용하여 $\Delta \mathbf{x}$를 직접 구하는 것이 일반적이다.

### 간단한 예시(Octave)

다변수 뉴턴 방법의 구현이 어떻게 동작하는지를 살펴보기 위해 간단한 예제를 Octave로 작성하면 아래와 같다.

```octave
% 예: F(x) = 0를 만족하는 x를 찾는 다변수 뉴턴 방법
% x = (x1, x2)^T
% f1(x1,x2) = x1^2 + x2^2 - 4
% f2(x1,x2) = exp(x1) + x1*x2 - 1
function y = f(x)
  x1 = x(1);
  x2 = x(2);
  y = [ x1^2 + x2^2 - 4;
        exp(x1) + x1*x2 - 1 ];
endfunction
function J = Jf(x)
  x1 = x(1);
  x2 = x(2);
  J = [ 2*x1, 2*x2;
        exp(x1) + x2, x1 ];
endfunction
x = [1; 1];  % 초기 추정값
tol = 1.0e-8;
maxiter = 50;
for k = 1:maxiter
  fx = f(x);
  if norm(fx) < tol
    break
  endif
  Jx = Jf(x);
  dx = - Jx \ fx;
  x = x + dx;
endfor
disp("근사해:")
disp(x)
disp("잔차:")
disp(norm(f(x)))
```

이 코드는 $\mathbf{x}\_0 = (1, 1)^T$를 초기값으로 두고, 다변수 뉴턴 방법을 통해 매 반복마다 야코비 행렬 $J(\mathbf{x}\_k)$를 계산하고 $\[J(\mathbf{x}\_k)]^{-1}\mathbf{f}(\mathbf{x}*k)$를 구한 뒤 $\mathbf{x}*{k+1} = \mathbf{x}\_k - \Delta \mathbf{x}\_k$ 방식으로 해를 갱신한다. 일정 기준(위 예제에서는 잔차 벡터의 놈이 $10^{-8}$ 이하)이 만족되면 반복을 중단한다.

### 야코비 역행렬 대신 직접 선형시스템 해 풀기

$\[J(\mathbf{x}\_k)]^{-1}$를 직접 구하는 것은 수치적으로 비효율적일 수 있다. 일반적으로는
$$

\begin{align} J(\mathbf{x}\_k),\Delta \mathbf{x} = -\mathbf{f}(\mathbf{x}\_k) \end{align}

$$
를 직접 푸는 과정을 통해 $\Delta \mathbf{x}$를 구한 뒤
$$

\begin{align} \mathbf{x}\_{k+1} = \mathbf{x}\_k + \Delta \mathbf{x} \end{align}

$$
형태로 갱신하는 것이 더 안정적이다. LU 분해나 QR 분해와 같은 고전적인 선형대수 기법 또는 희소 구조를 활용한 특수 알고리즘(고차원 문제에서 중요)을 적절히 적용하여 반복 효율을 높일 수 있다.

### 댐핑(damping) 기법과 수정 뉴턴 방법

다변수 뉴턴 방법은 국소적으로 빠른 수렴 특성을 가지지만, 초기값이 해 근방에 놓이지 않는 경우 발산하거나 수렴 속도가 저하될 수 있다. 이를 보완하기 위해 댐핑(damping) 혹은 선형 탐색(line search) 기법을 접목하여, 업데이트 벡터의 스텝 크기를 적절히 조절하는 방식을 적용할 수 있다.
댐핑 뉴턴(damped Newton) 알고리즘은
$$

\begin{align} \mathbf{x}\_{k+1} = \mathbf{x}\_k - \alpha\_k \[J(\mathbf{x}\_k)]^{-1} \mathbf{f}(\mathbf{x}\_k) \end{align}

$$
형태로 $\alpha\_k \in (0,1]$ 범위의 스칼라 계수를 적절히 선택한다. 일반적으로 $\alpha\_k$를 너무 크게 잡으면 오버슈팅(overshoot)으로 인한 발산이 발생할 수 있고, 너무 작게 잡으면 수렴 속도가 심각하게 저하된다. 따라서 적절한 $\alpha\_k$를 결정하기 위해 하향선형탐색(backtracking line search)나 해 발산을 방지하기 위한 체계적인 전략(예: Wolfe 조건, Armijo 조건 등)을 적용한다.

#### 선형 탐색(backtracking line search)

선형 탐색 방식은 $\alpha\_k$를 한 번에 결정하기 위해 특정한 목적 함수를 정의하고, 이 함숫값이 충분히 감소하는 $\alpha\_k$를 찾는 기법을 말한다. 예를 들어, 다음과 같은 표준적인 감소 조건(Armijo 조건)을 사용할 수 있다.
$$

\begin{align} |\mathbf{f}(\mathbf{x}\_k - \alpha\[J(\mathbf{x}\_k)]^{-1}\mathbf{f}(\mathbf{x}\_k))| \le (1 - c,\alpha) |\mathbf{f}(\mathbf{x}\_k)|\end{align}

$$
여기서 $c$는 $0 < c < 1$ 범위 내의 상수이고, $\alpha$는 적당히 큰 값에서 시작하여(예: $\alpha=1$) 점진적으로 줄여가며 조건을 만족할 때까지 찾는다. 이 과정을 통해 전역적 안정성과 수렴 가능성을 어느 정도 확보할 수 있다.

#### 수정 뉴턴 방법(modified Newton method)

수정 뉴턴 방법은 반복마다 야코비 행렬을 재계산하지 않고 일정 횟수마다만 갱신하거나, 야코비의 근사 행렬을 유지하면서 단계별로 조금씩 보정하여 계산 비용을 절약하는 방법이다. 특히 차원이 큰 문제에서, 매 반복마다 야코비을 정확하게 구하고 역행렬(혹은 이에 준하는 연산)을 진행하는 것은 연산 부담이 매우 커질 수 있다. 대표적으로는 다음 두 가지 접근이 자주 사용된다.

1. $J(\mathbf{x}\_0)$를 초기 스텝에서만 계산하고, 반복 과정에서는 동일한 야코비(또는 그 근사치)을 그대로 사용한다. 이 경우, 실제 뉴턴보다 수렴 속도가 느려지지만 반복당 계산 비용을 줄일 수 있어 전체 수행 시간 측면에서 이점이 있을 수 있다.
2. 야코비을 주기적으로 재계산한다(예: 5\~10번의 반복에 한 번씩). 또한 반복 과정 중에서 기울어진 정도가 커지거나 수렴 세기가 약해지면 그 시점에서 야코비을 다시 계산하여 정확도를 높인다.

#### 야코비의 조건수와 민감도

다변수 뉴턴에서 $\[J(\mathbf{x})]^{-1}$를 구하는 과정은 수치해석적으로 야코비 행렬의 조건수(condition number)와 밀접한 관련이 있다. 조건수가 큰 야코비 행렬은 작은 오차도 큰 업데이트 오차로 증폭시킬 수 있으므로, 수치 안정성이 악화된다. 이 경우에는
$$

\kappa(J(\mathbf{x})) = | J(\mathbf{x}) | ,| \[J(\mathbf{x})]^{-1} |

$$
가 크게 나타날 것이고, 실제로 $\Delta \mathbf{x}$가 불안정하게 계산될 수 있다. 이러한 민감도 문제는 선형 대수학적 기법(LU 분해 시의 부분 피벗팅, 정규 방정식 대신 QR 분해 도입 등)이나, 특정 규제(Reguralization) 기법을 적용하여 완화할 수 있다.

#### 대체 접근: 준-뉴턴(Quasi-Newton) 방법

실제 응용에서 야코비이 단순히 구할 수 없는 형태라면, 야코비 없이도 뉴턴 비슷한 수렴 특성을 확보하기 위한 준-뉴턴(Quasi-Newton) 계열 알고리즘을 사용할 수 있다. 가장 널리 쓰이는 방법으로는 Broyden류의 방법들이 있으며, 그중 BFGS(Broyden-Fletcher-Goldfarb-Shanno) 방법이 대표적이다. BFGS는 2차 도함수 정보를 직접 구하지 않아도, 반복 과정에서 얻는 $\mathbf{f}(\mathbf{x}\_{k+1}) - \mathbf{f}(\mathbf{x}*k)$와 $\mathbf{x}*{k+1} - \mathbf{x}\_k$의 변화량을 토대로 야코비 역행렬의 근사치(또는 헤시안 역행렬의 근사치)를 갱신한다.
만약 $B\_k \approx \[J(\mathbf{x}*k)]^{-1}$라고 할 때, $B*{k+1}$는 다음 조건을 만족하도록 업데이트된다.
$$

\begin{align} B\_{k+1} \bigl\[\mathbf{f}(\mathbf{x}\_{k+1}) - \mathbf{f}(\mathbf{x}*k)\bigr] = \mathbf{x}*{k+1} - \mathbf{x}\_k. \end{align}

$$
이와 같은 '계산된 변화량을 보존한다'는 제약을 최대로 만족하면서 $B\_{k+1}$를 구하는 방식이 Broyden 계열의 핵심 아이디어이다. BFGS 및 제한메모리 BFGS(L-BFGS) 방식들은 대규모 문제에도 적용할 수 있는 효율적 준-뉴턴 방법으로 잘 알려져 있다.

### 불완전 뉴턴(Inexact Newton) 방법

다변수 뉴턴 방법을 실제로 구현할 때, 반복마다 선형계 $J(\mathbf{x}\_k),\Delta \mathbf{x} = -\mathbf{f}(\mathbf{x}\_k)$를 아주 정밀하게(즉, 기계 오차 수준까지) 풀 필요가 없을 때가 있다. 과도하게 정밀한 해를 찾는 것은 반복마다의 비용이 증가하여 결과적으로 전체 알고리즘 수행 시간 측면에서 비효율적일 수 있다. 이를 보완하는 접근이 불완전 뉴턴(Inexact Newton) 방법이다.
불완전 뉴턴 방법은
$$

\begin{align} J(\mathbf{x}\_k),\Delta \mathbf{x}\_k \approx -\mathbf{f}(\mathbf{x}\_k) \end{align}

$$
를 근사적으로 풀되, 특정 허용오차(톨러런스) $\eta\_k$를 만족하도록 정한다. 예를 들어
$$

\begin{align} | J(\mathbf{x}\_k),\Delta \mathbf{x}\_k + \mathbf{f}(\mathbf{x}\_k)| \le \eta\_k,|\mathbf{f}(\mathbf{x}\_k)| \end{align}

$$
같은 조건이 만족되면 충분하다고 보고, $\Delta \mathbf{x}\_k$ 계산 정확도를 제한한다. 이렇게 하면 매 반복에서의 연산량을 대폭 절약할 수 있으며, 특히 $\eta\_k$를 점차 줄여가는 방식을 적용하면 국소 근방에서는 더욱 정밀한 계산을 수행하여 최종 수렴을 이룬다.

### 야코비 근사: 수치 미분과 자동 미분

다변수 뉴턴 방법에서 요구되는 야코비 행렬은 함수 $f\_i(\mathbf{x})$들의 편미분으로 구성된다. 어떤 경우에는 명시적인 편미분 식을 쉽게 구할 수 있으나, 복잡한 함수 구조나 외부 라이브러리를 연계한 경우라면 기호 미분(symbolic differentiation)은 물론이고 직접적인 편미분 코드 작성도 쉽지 않을 수 있다. 이때는 다음과 같은 기법을 활용한다.

#### 수치 미분(Finite Difference Approximation)

수치 미분은 편미분을 근사하기 위해 유한 차분 기법을 사용하는 방법이다. 예컨대, 일변수 함수 $g(x)$에 대한 전진 차분(forward difference) 근사는
$$

\begin{align} \frac{\partial g}{\partial x}(x) \approx \frac{g(x + h) - g(x)}{h}, \end{align}

$$
이고, 중심 차분(central difference) 근사는
$$

\begin{align} \frac{\partial g}{\partial x}(x) \approx \frac{g(x + h) - g(x - h)}{2h}. \end{align}

$$
다변수의 경우에도 편미분을 위와 유사한 방식으로 계산할 수 있다. 즉, $J(\mathbf{x})$의 $(i,j)$ 성분은
$$

\begin{align} \frac{\partial f\_i}{\partial x\_j}(\mathbf{x}) \approx \frac{f\_i(\mathbf{x} + h,\mathbf{e}\_j) - f\_i(\mathbf{x})}{h} \end{align}

$$
(전진 차분 기준)으로 근사할 수 있다. 여기서 $\mathbf{e}\_j$는 표준 기저벡터이고, $h$는 작은 스텝 크기이다. 중심 차분을 쓰면 오차차수가 더 높아지는 반면, 함수 평가가 그만큼 더 많이 필요해진다. $h$ 값 선택은 수치 오차(반올림 오차와 근사 오차의 절충) 관점에서 주의 깊게 결정해야 한다.

#### 자동 미분(Automatic Differentiation)

자동 미분(AD)은 컴퓨터 프로그램에 나타난 계산 그래프를 활용하여 미분 연산을 기계적으로 처리해주는 기법이다. 심볼릭 미분에 비해 표현식이 지나치게 커지거나 여러 특수함수를 다뤄야 하는 문제를 회피할 수 있으며, 수치 미분에 비해 오차가 극히 작고 계산 효율이 일반적으로 우수하다. 또, AD는 순방향(forward-mode)과 역방향(reverse-mode)으로 나누어지는데, 문제 차원이나 요구되는 미분의 개수에 따라 모드를 선택하면 이점이 있다.
예컨대 역방향 자동 미분은 스칼라 스칼라 함수의 그레디언트를 구할 때 효과적이고, 순방향 자동 미분은 스칼라 다차원 함수(다수의 종속변수)를 다룰 때 비교적 간편하게 적용 가능하다. 실제 응용에서도 C++, Python, MATLAB 등에서 AD 라이브러리를 이용해 손쉽게 야코비을 구할 수 있다.

### 초기값 설정과 반복 종료 기준

다변수 뉴턴 방법에서는 초기값 $\mathbf{x}\_0$에 따라 수렴성과 수렴 속도가 크게 달라진다. 비선형 시스템이 여러 해를 갖거나 복잡한 지형(topology)을 갖는 경우, 초기값을 어떻게 선택하느냐에 따라 전혀 다른 해로 수렴하거나 발산할 수도 있다.
반복 종료 기준은 다음과 같은 기준들을 조합하여 적용한다.
$$

|\mathbf{f}(\mathbf{x}*k)| < \varepsilon\_f,\quad |\mathbf{x}*{k+1} - \mathbf{x}*k| < \varepsilon\_x,\quad k > k*{\mathrm{max}},

$$
등을 고려해서 임의의 임계값(예: $10^{-8}$) 이하로 작아지면 중단한다. $k\_{\mathrm{max}}$는 반복 한계 횟수를 의미하고, 계산 자원 소모가 지나치게 커지는 것을 방지하기 위함이다.

### 전역적 수렴 전략(Global Convergence)

다변수 뉴턴 방법은 기본적으로 국소 수렴(local convergence)이 보장된다. 즉, 초깃값이 해 근방에 있을 때 이차 수렴(quadratic convergence)이 발휘되지만, 멀리 떨어진 지점에서는 발산하거나 해와 무관한 지점으로 치우칠 우려가 있다. 이때 전역적 수렴(최소한의 어느 정도 넓은 범위에서 해로 유도)을 어느 정도 확보하기 위해 댐핑 기법이나 신뢰영역(trust-region) 방법 등을 접목할 수 있다.

#### 신뢰영역(Trust-Region) 접근

신뢰영역 방법은 업데이트 방향과 크기를 결정할 때, 2차 근사(또는 1차 근사) 모델이 '믿을 만하다'고 간주되는 영역을 한정하는 방식이다. 예컨대, $\Delta \mathbf{x}$의 크기를 $\Delta$ 이하로 제한하고, $\Delta$를 조절하는 반복 알고리즘을 통해 수렴을 안정화한다. 소위 Levenberg-Marquardt 방법도 신뢰영역과 유사한 원리를 활용한다.
신뢰영역 접근의 기본 아이디어는 다음과 같은 보조 문제를 매 스텝 풀어 $\Delta \mathbf{x}\_k$를 구하는 것이다.
$$

\begin{align} \min\_{\Delta \mathbf{x}} \quad m\_k(\Delta \mathbf{x}) = |\mathbf{f}(\mathbf{x}\_k) + J(\mathbf{x}\_k),\Delta \mathbf{x}|^2,\end{align}

$$
단,
$$

\begin{align} |\Delta \mathbf{x}| \le \Delta\_k.\end{align}

$$
여기서 $\Delta\_k$는 신뢰영역의 반경이다. $\Delta \mathbf{x}$를 구하고 나면, 실제 함수값의 개선 정도와 모델의 예측 개선 정도를 비교하여 $\Delta\_k$를 늘리거나 줄이는 전략을 취함으로써 전역적 안정성을 높인다.

### 복합 방정식 구조와 실전 적용상의 유의점

실제 대규모 문제에서 다변수 뉴턴 방법은 선형 시스템 풀기에 드는 비용이 매우 커질 수 있다. 행렬-벡터 연산이 지배적인 상황에서, 다음 사항들이 성능과 안정성에 중요한 역할을 한다.
직교 분해(QR), LU, 혹은 Cholesky(양의 정부호 행렬인 경우) 등으로 분해할 때, 부분 피벗팅(partial pivoting)이나 완전 피벗팅(full pivoting)을 적용한다면 수치 안정성을 크게 개선할 수 있다.
희소(sparse) 행렬 구조를 최대한 활용하기 위해선, 희소 LU나 Krylov 하부공간 기법(GMRES, CG 등)을 이용한 반복적 선형계 풀기를 병행할 수 있다.
야코비을 직접 구할 수 없는 경우, 수치 미분 또는 자동 미분을 적극 검토해야 한다.
비선형 시스템이 상호결합된 복합방정식(예: 유체역학, 전자기학 연립방정식 등)으로 나타날 때는 문제에 특화된 전처리(Preconditioning)를 통해 해 검색 효율을 높일 수 있다.
$$