正确使用 skopt.gp_minimize 进行优化与结果可视化

import numpy as np
from scipy.optimize import minimize
from skopt import gp_minimize
import matplotlib.pyplot as plt

# 辅助函数（与原问题代码保持一致，此处省略详细定义，但在完整代码中会包含）
def gaussian_rbf(x, x_prime, beta):
    return np.exp(-beta * np.linalg.norm(x - x_prime)**2)

def construct_interpolation_matrix(nodes, beta):
    N = len(nodes)
    K = np.zeros((N, N))
    for i in range(N):
        for j in range(N):
            K[i, j] = gaussian_rbf(nodes[i], nodes[j], beta)
    return K

def conditioning_analysis(N, m, beta):
    nodes = np.linspace(0, 1, N)
    K = construct_interpolation_matrix(nodes, beta)

    selected_indices = np.random.choice(N, m, replace=False)
    selected_nodes = nodes[selected_indices]

    condition_full = np.linalg.cond(K)
    condition_partial = np.linalg.cond(K[selected_indices][:, selected_indices])

    return condition_full, condition_partial

# 目标函数：应能处理单个标量输入
def objective_function(x):
    # 确保x是标量，对于numpy数组也兼容
    x_scalar = np.atleast_1d(x)[0] if np.ndim(x) > 0 else x
    return -(x_scalar**2 + np.sin(5 * x_scalar))

# 牛顿法相关的梯度和Hessian（与原问题代码保持一致）
def gradient_hessian(x):
    # 注意：原始代码中的梯度和Hessian函数与objective_function不匹配，
    # 原始的objective_function是 -(x^2 + sin(5x))
    # 原始的gradient_hessian似乎是为 f(x) = x * exp(-(1-x)^2) 编写的。
    # 为保持教程的焦点，此处使用原始的gradient_hessian，但请注意此潜在不一致。
    # 正确的梯度和Hessian应为：
    # df_dx = - (2 * x + 5 * np.cos(5 * x))
    # d2f_dx2 = - (2 - 25 * np.sin(5 * x))
    # 为避免引入新的复杂性，此处沿用原代码中的gradient_hessian，但建议用户根据实际目标函数进行修正。
    df_dx = 2 * x * np.exp(-(1 - x)**2) - 4 * x * (1 - x) * np.exp(-(1 - x)**2)
    d2f_dx2 = -2 * np.exp(-(1 - x)**2) + 4 * x * (1 - x) * np.exp(-(1 - x)**2) - 4 * (1 - x) * np.exp(-(1 - x)**2)
    return df_dx, d2f_dx2

def optimize_with_newton(initial_guess, max_iter=10):
    x_opt = initial_guess
    for _ in range(max_iter):
        df_dx, d2f_dx2 = gradient_hessian(x_opt)
        # 避免除以零或非常小的数
        if abs(d2f_dx2) < 1e-9:
            print(f"Warning: Hessian near zero at x={x_opt}, stopping Newton iteration.")
            break
        x_opt = x_opt - df_dx / d2f_dx2
    return x_opt

# 修正后的高斯过程优化函数
def gaussian_process_optimization(initial_points, objective_function, bounds, n_iter=10):
    """
    对每个初始点独立运行 gp_minimize，并返回所有优化结果。
    """
    optimal_x_values = np.zeros(len(initial_points))
    for i, x0_val in enumerate(initial_points):
        # 关键修正：x0 必须是包含单个点的列表，例如 [0.5]
        # 并且直接传入原始的 objective_function，它应处理标量输入
        result = gp_minimize(objective_function, bounds, acq_func="LCB", n_calls=n_iter + 1, random_state=42 + i, x0=[x0_val])
        optimal_x_values[i] = result.x[0] # result.x 是一个列表，取第一个元素
    return optimal_x_values

# Task 1: Analyze conditioning (与原问题代码一致)
N = 10
m = 5
beta = 1.0
condition_full, condition_partial = conditioning_analysis(N, m, beta)
print(f"Conditioning for full matrix: {condition_full}")
print(f"Conditioning for partial matrix: {condition_partial}")

# 优化与牛顿法
initial_guess_newton = 0.5
x_opt_newton = optimize_with_newton(initial_guess_newton)
print(f"Optimal solution with Newton's method: {x_opt_newton}")

# 高斯过程优化
initial_points_gp = np.random.rand(5) # 5个随机初始点
bounds_gp = [(0.0, 1.0)]

# 调用修正后的高斯过程优化函数
x_opt_gp_array = gaussian_process_optimization(initial_points_gp, objective_function, bounds_gp, n_iter=10)
print(f"Optimal solutions with Gaussian process optimization (from multiple starts): {x_opt_gp_array}")

# 从多个GP结果中选择最佳点进行可视化，以便与牛顿法进行直接比较
# 假设我们寻找的是最小值，因此选择 objective_function 值最小的点
y_values_gp = [objective_function(x) for x in x_opt_gp_array]
best_gp_index = np.argmin(y_values_gp)
best_x_opt_gp = x