Newton's Method (Quadratic interpolation with 2nd derivative) 설명 (python 구현)

2024.02.02 - [데이터 분석/최적화 알고리즘] - Golden Section Search 설명 (python 구현)

 

Newton's Method (Quadratic interpolation with 2nd derivative)

방정식 $f(x) = 0$이 되게 하는 $x$를 찾는 알고리즘

 

 

수식

함수 $f(x)$가 해 $x*$ 근처에서 연속이고 두 번 미분가능할 때
- 임의의 점 $x_0$에서 테일러 급수를 2차 함수까지 전개
- 해당 2차 함수의 근을 $x_1$으로 삼아 반복
- 또는 2차 함수가 근을 갖지 않는 경우, local minimum이 되게 하는 $x$를 $x_1$으로 삼아 반복
- $|x_{n+1}-x_{n}| \le \text{tol}$일 때까지 반복, `tol`은 임의의 작은 수

 

 

장단점

- 테일러 급수를 더 높은 차수까지 전개하면 더 빠르게 해를 구할 수 있다.
- 5차함수부터는 일반해를 구하는 공식이 없기 때문에 원래 해를 구하기 위해 이 해도 구해야 한다.
- 미분이 가능한 함수에만 사용할 수 있다.

 

 

구현 소스 코드

def derivative(func, x):
    delta_x = 1e-4
    fx_plus_delta = func(x + delta_x)
    fx_minus_delta = func(x - delta_x)
    
    return (fx_plus_delta - fx_minus_delta) / (2 * delta_x)

def twoderivative(func, x):
    delta_x = 1e-4
    return (func(x+2*delta_x) - 2*func(x+delta_x) + func(x)) / (delta_x**2)

class newton:
    
    def __init__(self, xn, func, tol=1e-6):
        self.xn = xn
        self.func = func
        self.tol = tol
        
        
    def solve(self):
        def D(a, b, c):
            return b**2 - 4*a*c

        def solve_quadratic(a, b, c): 
            d = D(a, b, c)
            xu = (-b + d**0.5) / (2*a)
            xl = (-b - d**0.5) / (2*a)
            return xl, xu
        
        history = [self.xn]
        
        while True:
            f_prime_xn = derivative(self.func, self.xn)
            f_twoprime_xn = twoderivative(self.func, self.xn)
            # 이계도함수가 0에 가까운 경우, noise 더하기
            if abs(f_twoprime_xn) <= 1e-1:
                f_twoprime_xn += 1e-1
            
            a = f_twoprime_xn / 2
            b = f_prime_xn - self.xn * f_twoprime_xn
            c = self.xn ** 2 * f_twoprime_xn / 2 - self.xn * f_prime_xn + self.func(self.xn)
            
            d = D(a, b, c)
            
            # 이차함수의 근이 존재하는 경우
            if d>0 :
                # xl, xu 중 xn에 더 가까운 값을 xn+1로 채택
                xl, xu = solve_quadratic(a, b, c)
                
                if abs(xl - self.xn) < abs(xu - self.xn):
                    xn1 = xl
                else:
                    xn1 = xu
            
            # 이차함수의 근이 없는 경우 = local minimum = 기존 newton's method 방식대로 구하기
            else:
                # 미분값이 0에 매우 가까운 경우, xn에 noise를 더하여 빠져나오기
                if abs(f_prime_xn) <= 1e-1:
                    self.xn += self.tol
                    continue

                xn1 = self.xn - self.func(self.xn) / f_prime_xn
            
            # fxn이 tol 이하인 경우
            if abs(self.func(xn1)) <= self.tol:
                self.solution = xn1
                self.history = history
                return xn1
            
            self.xn = xn1
            history.append(xn1)

 

 

 

예제

$f(x) = \sin ({ \cos ({ \rm{e}^x }) })$

NM = newton(func=f, xn = 0.2, tol=1e-10)
NM.solve() # 0.45158270529849714

f(NM.solution) # -1.4203688055504446e-11

 

$g(x) = x^3 - x$

derivative(g, 1/math.sqrt(3)) # 1.0000056338554941e-08

NM = newton(func=g, xn = 1/math.sqrt(3), tol=1e-10)
NM.solve() # 0.9999999999992966

NM.history # [0.5773502691896258, 1.04871396685623, 0.9999416299477963]

 

전체 소스 코드 → github