Skip to content

Instantly share code, notes, and snippets.

@eccstartup
Last active March 23, 2022 09:17
Show Gist options
  • Save eccstartup/9c141d0f54dfb59a573331555f7feb39 to your computer and use it in GitHub Desktop.
Save eccstartup/9c141d0f54dfb59a573331555f7feb39 to your computer and use it in GitHub Desktop.
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment