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

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

import numpy as np
from numpy.linalg import norm
from scipy.optimize import minimize_scalar

def smooth_newton(trust_radius, x0, f, g, h):
    x = x0
    while True:
        # Compute the gradient and Hessian at the current point
        grad = g(x)
        H = h(x)
        
        # Solve the trust region subproblem
        p = compute_step(grad, H, trust_radius)
        x_new = x + p
        
        # Compute the actual reduction and predicted reduction
        actual_reduction = f(x) - f(x_new)
        predicted_reduction = -grad.dot(p) - 0.5 * p.dot(H.dot(p))
        
        # Compute the ratio of actual reduction to predicted reduction
        rho = actual_reduction / predicted_reduction
        
        # Update the trust region radius
        if rho < 0.25:
            trust_radius /= 2
        elif rho > 0.75 and norm(p) == trust_radius:
            trust_radius = min(2*trust_radius, 1e8)
        
        # Update the iterate
        if rho > 0:
            x = x_new
        
        # Check for convergence
        if norm(grad) < 1e-6:
            break
    
    return x

def compute_step(grad, H, trust_radius):
    # Solve the trust region subproblem using the dogleg method
    p_u = -np.linalg.solve(H, grad)
    if norm(p_u) <= trust_radius:
        return p_u
    else:
        p_b = -grad.dot(grad) / grad.dot(H.dot(grad)) * grad
        if norm(p_b) >= trust_radius:
            return trust_radius / norm(grad) * grad
        else:
            a = norm(p_u - p_b)**2
            b = 2 * p_b.dot(p_u - p_b)
            c = norm(p_b)**2 - trust_radius**2
            tau = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a)
            return p_b + tau * (p_u - p_b)

# Example usage
def f(x):
    return x[0]**2 + x[1]**2

def g(x):
    return np.array([2*x[0], 2*x[1]])

def h(x):
    return np.array([[2, 0], [0, 2]])

x0 = np.array([1, 1])
trust_radius = 1.0
x = smooth_newton(trust_radius, x0, f, g, h)
print("Optimal solution:", x)

在上面的代码中，`smooth_newton`函数实现了光滑牛顿法的迭代过程，包括求解信赖域子问题、更新信赖域半径、更新迭代点等步骤。`compute_step`函数实现了用狗腿方法求解信赖域子问题的过程。在示例中，我们使用了一个简单的二次函数作为目标函数，以演示光滑牛顿法的使用方法。