newton.py

from sympy import *
import numpy as np


x = IndexedBase('x')
fun0 = (x[0] + x[1] - 10) ** 2 + (x[0] - 2 * x[1] + 5) ** 2 - 100


def grad(f, x0):
    n = x0.size
    grads = np.zeros(n)
    for i in range(n):
        grads[i] = lambdify(x, diff(f, x[i]), 'numpy')(x0)
    return grads


def hess(f, x0):
    n = x0.size
    hesses = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            hesses[i, j] = lambdify(x, diff(f, x[i], x[j]), 'numpy')(x0)
    return hesses


def newton(f, init, num_iters, verbose=False):
    x0 = init
    fun = lambdify(x, f, "numpy")
    y = fun(init)
    costs = [y]
    if verbose:
        print(f'开始迭代, x: {x0}, cost : {y}')
    for i in range(num_iters):
        x1 = x0 - np.dot(np.linalg.inv(hess(f, x0)), grad(f, x0))
        y = fun(x1)
        if verbose:
            print(f'第{i+1}此迭代, x: {x1}, cost: {y}')
        x0 = x1
        costs.append(y)
        change = np.linalg.norm(x1-x0)
        if change <= 0.000001:
            break
    xopt = x0
    return xopt, costs


newton(fun0, np.array([1, 1]), 10, True)