이번 글에서는 Newton’s method (뉴턴의 방법)을 구현한 코드를 공개하도록 하겠다. Newton’s method 는 워낙 유명한 방법이니까 알고리즘 설명은 생략한다. 진짜 궁금하다면 위키에 있는 Newton’s method 에 대한 설명을 보도록 하자. 그러면 이제 본격적으로 Newton’s method (뉴턴의 방법)에 대해 잘 알아보도록 하겠습니다요.
Newton’s method (뉴턴의 방법) 구현
inverse matrix를 구해서 업데이트 하는 부분이 필요한데 나는 그것보다 gaussian eliminiation 을 이용해서 해를 구한 후 업데이트 하는 방식을 택했다. 이를 위해 gaussian elimination 방법을 택했다.
Gaussian elimination 구현
이제 이 Gaussian elimination 을 바탕으로 Newton’s method (뉴턴의 방법)을 아래와 같이 구현하였다.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
|
import numpy as np
import matplotlib.pyplot as plt
# partially pivoting
def gaussian_elimination(A, b):
#Ax=b 의 solution 구하기
#입력변수 A,b
n = len(b)
aug_mat = np.hstack((A, b.reshape(–1, 1)))
for i in range(n):
pivot_row = np.argmax(np.abs(aug_mat[i:, i])) + i
aug_mat[[i, pivot_row]] = aug_mat[[pivot_row, i]]
pivot = aug_mat[i, i]
aug_mat[i] /= pivot
for j in range(i + 1, n):
divided_factor = aug_mat[j, i] / aug_mat[i, i]
aug_mat[j] –= divided_factor * aug_mat[i]
# back substitution
x = np.zeros(n)
for i in range(n – 1, –1, –1):
x[i] = aug_mat[i, –1] – np.dot(aug_mat[i, i+1:n], x[i+1:])
return x
|
cs |
Newton’s method (뉴턴의 방법) 구현 코드
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
|
# Newton’s method
def newton_method(F, J, initial_point, tolerance=1e–8, max_iter=100):
#Newton Method이용하기
#F: 함수, J: F의 Jacobian matrix, initial_point: 초기 설정값, tolerance: 0이와 가까운지 판단하기 위함 임계치, max_iter: 최대 iteration 횟수
x = initial_point
history = [x.copy()]
for iter in range(1, max_iter + 1):
F_value = F(x)
J_value = J(x)
#partially pivoting
delta_x = gaussian_elimination(J_value, –F_value)
x = x + delta_x
history.append(x.copy())
# convergence criterion
if np.linalg.norm(F(x)) < tolerance:
break
return x, np.array(history)
|
cs |
실제로 Newton’s method (뉴턴의 방법)을 사용하는 방법은 아래와 같이 함수를 정의하고 이함수의 gradient를 직접 입력해주면 된다.
Newton’s method (뉴턴의 방법) 예시
1
2
3
4
5
6
7
8
9
10
11
|
# 예제 1: F(x1, x2) = (x1^2 + x2^2 – 2, x1^2 – x2^2 – 1)
def func_for_3(x):
return np.array([x[0]**2 + x[1]**2 – 2, x[0]**2 – x[1]**2 – 1])
def jacob_3(x):
return np.array([[2*x[0], 2*x[1]], [2*x[0], –2*x[1]]])
# 초기 추정치
initial_guess1 = np.array([1.5, 0.5])
solution, history = newton_method(F, J, initial_guess)
|
cs |
추천글