• 1. 학습 목표
• 2. 비선형 연립방정식
• 3. Jacobian 행렬
• 4. 다변수 Newton–Raphson 방법
  • 4.1. 선형화
  • 4.2. 알고리듬
• 5. 손으로 계산하는 첫 번째 반복
• 6. Python 구현
  • 6.1. 두 근 구하기
  • 6.2. 수렴 기록 확인
• 7. 세 변수로 확장
• 8. 수렴과 주의 사항
• 9. 경사하강법과의 차이
• 10. 정리
• 11. 연습 문제

1. 학습 목표

이번 강의가 끝나면 다음을 할 수 있어야 한다.

• 비선형 연립방정식을 벡터 함수 $\boldsymbol F(\boldsymbol x)=\boldsymbol0$으로 나타낼 수 있다.
• Jacobian 행렬을 편미분으로 구성할 수 있다.
• 다변수 Newton–Raphson 방법의 선형 시스템을 유도할 수 있다.
• 역행렬을 직접 계산하지 않고 갱신량을 구할 수 있다.
• Python으로 다변수 비선형 연립방정식을 풀 수 있다.
• 경사하강법과 다변수 Newton–Raphson 방법의 목적을 구분할 수 있다.

2. 비선형 연립방정식

두 변수 $x,y$에 관한 두 방정식을 생각해보자.

모든 항을 왼쪽으로 옮기면

이다. 이를 하나의 벡터 함수로 나타내면

우리가 찾는 해는

을 만족한다. 각 방정식은 평면에서 하나의 곡선을 나타내며, 두 곡선의 교점이 연립방정식의 해이다. 이 예제에는 $(2,1)$과 $(-1,-2)$라는 두 해가 있다.

3. Jacobian 행렬

Jacobian 행렬은 벡터 함수의 모든 일계 편미분을 모은 행렬이다.

두 변수의 경우

앞의 예제에서는

Jacobian의 $i$번째 행은 함수 $f_i$의 기울기 벡터를 전치한 것이다.

4. 다변수 Newton–Raphson 방법

4.1. 선형화

현재 근삿값을 $\boldsymbol x^{(k)}$, 작은 변화량을 $\Delta\boldsymbol x^{(k)}$라고 하자. 벡터 함수의 1차 Taylor 근사는

이다. 새로운 점에서 함수값이 영벡터가 되도록 왼쪽을 0으로 놓으면

을 얻는다. 이 선형 연립방정식을 풀어 변화량을 구한 다음

로 갱신한다.

$J^{-1}$을 명시적으로 계산하여 곱하기보다 가우스 소거법, LU 분해 또는 np.linalg.solve로 선형 시스템을 직접 푸는 것이 좋다.

4.2. 알고리듬

1. 초기값 $\boldsymbol x^{(0)}$와 허용오차를 정한다.
2. $\boldsymbol F(\boldsymbol x^{(k)})$와 $\boldsymbol J(\boldsymbol x^{(k)})$를 계산한다.
3. $|\boldsymbol F(\boldsymbol x^{(k)})|_2$가 충분히 작으면 종료한다.
4. $J\Delta\boldsymbol x=-F$를 풀어 변화량을 구한다.
5. $\boldsymbol x^{(k+1)}=\boldsymbol x^{(k)}+\Delta\boldsymbol x$로 갱신한다.
6. 최대 반복 횟수 안에서 2–5단계를 반복한다.

5. 손으로 계산하는 첫 번째 반복

초기값을

로 정하자. 함수 벡터와 Jacobian은

이다. 따라서 변화량은

을 만족한다. 두 번째 식에서 $\Delta x=\Delta y$이고, 첫 번째 식에 대입하면

이므로

이다. 다음 근삿값은

6. Python 구현

import numpy as np

def multivariate_newton(func, jacobian, x0, tol=1e-10, max_iter=50):
    x = np.asarray(x0, dtype=float).copy()
    history = [x.copy()]
    residual_history = [np.linalg.norm(func(x))]

    for iteration in range(max_iter+1):
        F = np.asarray(func(x), dtype=float)

        if np.linalg.norm(F) <= tol:
            return x, iteration, np.array(history), np.array(residual_history)

        if iteration == max_iter:
            break

        J = np.asarray(jacobian(x), dtype=float)

        try:
            dx = np.linalg.solve(J, -F)
        except np.linalg.LinAlgError as error:
            raise np.linalg.LinAlgError(
                "Jacobian이 특이하여 변화량을 계산할 수 없습니다."
            ) from error

        x = x+dx
        history.append(x.copy())
        residual_history.append(np.linalg.norm(func(x)))

    raise RuntimeError("최대 반복 횟수 안에 수렴하지 않았습니다.")

예제의 함수와 Jacobian을 정의한다.

def func(x):
    x1, x2 = x
    return np.array([
        x1**2+x2**2-5,
        x1-x2-1,
    ])

def jacobian(x):
    x1, x2 = x
    return np.array([
        [2*x1, 2*x2],
        [1.0,   -1.0],
    ])

6.1. 두 근 구하기

초기값을 바꾸어 두 근을 구해보자.

for x0 in ([1.5, 0.5], [-0.5, -1.5]):
    root, iterations, history, residuals = multivariate_newton(
        func, jacobian, x0
    )
    print("initial:", x0)
    print("root:", root)
    print("iterations:", iterations)
    print("residual norm:", np.linalg.norm(func(root)))

첫 번째 초기값에서는 $(2,1)$, 두 번째 초기값에서는 $(-1,-2)$에 수렴한다.

6.2. 수렴 기록 확인

import matplotlib.pyplot as plt

root, iterations, history, residuals = multivariate_newton(
    func, jacobian, [1.5, 0.5]
)

plt.semilogy(range(len(residuals)), residuals, "o-")
plt.xlabel("iteration")
plt.ylabel(r"$\|\boldsymbol F(\boldsymbol x^{(k)})\|_2$")
plt.grid()
plt.show()

7. 세 변수로 확장

변수가 세 개이면 함수 벡터와 Jacobian도 세 개의 행과 열을 갖는다. 예를 들어

이면

$(1,1,1)$은 $F(1,1,1)=0$을 만족하는 해이다. 앞서 작성한 함수는 벡터의 차원에 관계없이 같은 방법으로 사용할 수 있다.

def func3(x):
    x1, x2, x3 = x
    return np.array([
        x1**2+x2+x3-3,
        x1+x2**2+x3-3,
        x1+x2+x3**2-3,
    ])

def jacobian3(x):
    x1, x2, x3 = x
    return np.array([
        [2*x1, 1.0, 1.0],
        [1.0, 2*x2, 1.0],
        [1.0, 1.0, 2*x3],
    ])

root, iterations, history, residuals = multivariate_newton(
    func3, jacobian3, [0.8, 1.2, 0.9]
)
print(root)

8. 수렴과 주의 사항

단일변수 Newton–Raphson 방법과 마찬가지로 단순근에 충분히 가까운 초기값에서는 빠른 수렴을 기대할 수 있다. 다변수 문제의 단순근에서는 해에서의 Jacobian이 가역행렬이어야 한다.

• 초기값에 따라 서로 다른 근으로 수렴할 수 있다.
• Jacobian이 특이하거나 거의 특이하면 변화량을 구하기 어렵다.
• 초기값이 해에서 멀면 발산하거나 예상하지 않은 근으로 이동할 수 있다.
• 함수와 Jacobian의 행·열 순서를 서로 일치시켜야 한다.
• 잔차와 변화량을 확인하고 최대 반복 횟수를 설정해야 한다.
• 실제 문제에서는 변화량에 감쇠계수 $0<\lambda\leq1$을 곱하는 감쇠 Newton 방법을 사용하여 안정성을 높이기도 한다.

9. 경사하강법과의 차이

10. 정리

• 비선형 연립방정식은 $F(x)=0$인 벡터 근 찾기 문제이다.
• Jacobian의 각 원소는 $J_{ij}=\partial f_i/\partial x_j$이다.
• 매 반복에서 $J\Delta x=-F$를 풀고 $x$에 $\Delta x$를 더한다.
• 역행렬을 직접 구하지 않고 선형 연립방정식을 푼다.
• 초기값과 Jacobian의 가역성이 수렴에 중요한 영향을 준다.

11. 연습 문제

1. 다음 벡터 함수의 Jacobian을 구하라.

   \[F(x,y)= \begin{bmatrix} x^2+y\\ x-y^2 \end{bmatrix}.\]

2. 예제의 $F(x,y)=[x^2+y^2-5,\ x-y-1]^T$에서 $F(2,1)$을 계산하라.

3. 같은 예제에서 $J(2,1)$을 구하라.

4. 초기값 $(1.5,0.5)$에서 풀어야 하는 선형 시스템의 우변 $-F$를 구하라.

5. 변화량이 $\Delta x=(0.2,-0.1)^T$이고 현재점이 $x=(1,2)^T$이면 다음 근삿값을 구하라.

6. 다음 코드의 빈칸을 채워 변화량을 구하라.

   dx = np.linalg.solve(____, ____)
   

7. Jacobian이 특이행렬이면 Newton 갱신이 어려운 이유를 한 문장으로 설명하라.

8. 경사하강법과 다변수 Newton–Raphson 방법의 목적을 각각 적어라.