Python目标函数优化异常求助：参数未按预期收敛

现有长度为600的向量V，需最小化包含两个未知参数的目标函数：
$$\text{objective}(x) = \sum_{i=0}^{599} \left(V[i] - \left(x_0 + \frac{x_1}{\sqrt{\text{horizontalsteps}[i]}}\right)\right)^2$$
尝试过trust-constr、Nelder-Mead、SLSQP等多种优化方法，但所有方法表现相近；即使以预期最优参数作为初始猜测，优化结果也未收敛到该值。核心需求是找到使目标函数值最小的参数，不限方法。

你的目标函数本质是线性回归模型：令$z_i = \frac{1}{\sqrt{\text{horizontalsteps}[i]}}$，模型可改写为$V[i] = x_0 + x_1 \cdot z_i$，这是标准的一元线性回归问题，存在解析解，完全不需要使用通用非线性优化方法。通用优化方法可能因迭代收敛条件、数值精度偏差等问题，无法精准触达全局最优解（尽管线性问题的损失函数是凸的，但数值迭代仍可能存在误差）。

方法1：直接计算解析解 #

线性回归的参数可通过最小二乘解析解直接求解，公式为：
$$\begin{bmatrix} x_0 \ x_1 \end{bmatrix} = (X^T X)^{-1} X^T V$$
其中$X$是设计矩阵，每一行是$[1, z_i]$。

代码实现：

import numpy as np

# 假设V已加载完成
batch_size = 100
horizontalsteps = np.arange(1, 60000, batch_size)
z = 1 / np.sqrt(horizontalsteps)

# 构造设计矩阵
X = np.column_stack((np.ones(len(z)), z))

# 计算解析解
params = np.linalg.inv(X.T @ X) @ X.T @ V
x0_opt, x1_opt = params
print("解析解优化参数：", x0_opt, x1_opt)

方法2：使用专用线性最小二乘函数 #

scipy.linalg.lstsq专为最小二乘问题设计，比通用优化方法更高效、准确：

import numpy as np
from scipy.linalg import lstsq

z = 1 / np.sqrt(horizontalsteps)
X = np.column_stack((np.ones(len(z)), z))

# 求解最小二乘
params, residuals, rank, singular_values = lstsq(X, V)
x0_opt, x1_opt = params
print("线性最小二乘优化参数：", x0_opt, x1_opt)
print("残差平方和：", residuals[0])

方法3：若坚持使用minimize（不推荐） #

如果一定要用scipy.optimize.minimize，可给目标函数添加梯度信息（线性问题的梯度极易计算），大幅提升收敛速度与精度：

import numpy as np
from scipy.optimize import minimize

def objective(x):
    z = 1 / np.sqrt(horizontalsteps)
    preds = x[0] + x[1] * z
    return np.sum((V - preds) ** 2)

def gradient(x):
    z = 1 / np.sqrt(horizontalsteps)
    preds = x[0] + x[1] * z
    res = V - preds
    grad0 = -2 * np.sum(res)
    grad1 = -2 * np.sum(res * z)
    return np.array([grad0, grad1])

initial_guess = [0.003, 6]
param_bounds = [(-10000, 10000), (-10000, 10000)]

# 使用带梯度的L-BFGS-B方法
result = minimize(objective, initial_guess, method='L-BFGS-B', bounds=param_bounds, jac=gradient)
print("带梯度的优化参数：", result.x)

计算出最优参数后，可绘制拟合曲线验证效果：

import matplotlib.pyplot as plt

preds = x0_opt + x1_opt / np.sqrt(horizontalsteps)
plt.plot(horizontalsteps, V, label="原始数据")
plt.plot(horizontalsteps, preds, label="拟合曲线", linestyle="--")
plt.xlabel("horizontalsteps")
plt.ylabel("V")
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者p.golestaneh