光滑牛顿法求解信赖域子问题python

以下是使用Python实现光滑牛顿法求解信赖域子问题的代码示例：

import numpy as np
from scipy.optimize import minimize_scalar

def solve_subproblem(B, g, delta):
    n = B.shape[0]
    p = np.linalg.solve(B, -g)
    if np.linalg.norm(p) <= delta:
        return p
    else:
        def fun(alpha):
            return np.dot(g, np.linalg.solve(B, g + alpha * p))
        alpha = minimize_scalar(fun, bounds=(0, np.inf), method='bounded').x
        return alpha * p / np.linalg.norm(p)

def smooth_newton_method(f, x0, grad_f, hess_f, delta=1.0, eta=0.1, max_iter=100):
    x = x0
    for i in range(max_iter):
        g = grad_f(x)
        B = hess_f(x)
        p = solve_subproblem(B, g, delta)
        q = f(x) - f(x + p)
        pred = -np.dot(g, p) - 0.5 * np.dot(p, np.dot(B, p))
        rho = q / pred
        if rho < eta:
            delta *= 0.5
        else:
            delta = min(2 * delta, 1.0)
        if rho > 0:
            x += p
        if np.linalg.norm(g) < 1e-6:
            break
    return x

其中，`B`是当前点的黑塞矩阵，`g`是当前点的梯度向量，`delta`是信赖域半径，`p`是求解的步长向量，`f`是目标函数，`x0`是初始点，`grad_f`是目标函数的梯度向量函数，`hess_f`是目标函数的黑塞矩阵函数，`eta`是控制信赖域半径的缩放因子，`max_iter`是最大迭代次数。

在求解过程中，首先根据当前点的黑塞矩阵和梯度向量求解信赖域子问题，得到步长向量`p`。然后根据预测准确率和实际准确率的比值更新信赖域半径。如果预测准确率比实际准确率低，则缩小信赖域半径；如果预测准确率比实际准确率高，则扩大信赖域半径。最后根据步长向量更新当前点。如果梯度向量的范数小于某个阈值，则认为已经达到最优解，停止迭代。