光滑牛顿法python程序

以下是一种实现光滑牛顿法的Python程序：

import numpy as np

def smooth_newton_method(f, x0, eps=1e-6, max_iter=100):
    x = x0
    for i in range(max_iter):
        fx = f(x)
        grad_f = gradient(f, x)
        hess_f = hessian(f, x)
        p = np.linalg.solve(hess_f, -grad_f)
        alpha = backtracking_line_search(f, x, p, grad_f)
        x_next = x + alpha * p
        if np.linalg.norm(x_next - x) < eps:
            return x_next
        x = x_next
    return x

def gradient(f, x):
    h = 1e-6
    grad_f = np.zeros_like(x)
    for i in range(len(x)):
        x_plus = x.copy()
        x_minus = x.copy()
        x_plus[i] += h
        x_minus[i] -= h
        grad_f[i] = (f(x_plus) - f(x_minus)) / (2 * h)
    return grad_f

def hessian(f, x):
    h = 1e-6
    n = len(x)
    hess_f = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            x_pp = x.copy()
            x_pp[i] += h
            x_pp[j] += h
            x_pm = x.copy()
            x_pm[i] += h
            x_pm[j] -= h
            x_mp = x.copy()
            x_mp[i] -= h
            x_mp[j] += h
            x_mm = x.copy()
            x_mm[i] -= h
            x_mm[j] -= h
            hess_f[i, j] = (f(x_pp) - f(x_pm) - f(x_mp) + f(x_mm)) / (4 * h ** 2)
    return hess_f

def backtracking_line_search(f, x, p, grad_f, alpha=1, rho=0.9, c=1e-4):
    while f(x + alpha * p) > f(x) + c * alpha * np.dot(grad_f, p):
        alpha *= rho
    return alpha

其中，`smooth_newton_method`函数是主函数，它接受一个目标函数`f`、一个初始点`x0`、一个可选的精度`eps`和一个可选的最大迭代次数`max_iter`，并返回一个近似最小值点。该函数在每次迭代中计算目标函数的梯度和海森矩阵，并使用求解线性方程组的方法计算步长，然后使用回溯线搜索确定步长大小。

`gradient`函数计算目标函数在点`x`处的梯度，使用中心差分法近似计算。

`hessian`函数计算目标函数在点`x`处的海森矩阵，同样使用中心差分法近似计算。

`backtracking_line_search`函数使用回溯线搜索确定步长大小，其中`alpha`是初始步长大小，`rho`是步长缩放因子，`c`是Armijo条件中的常数。