带bounds的scipy.curve_fit拟合mpmath.zeta时遇NaN错误

问题背景 #

使用scipy.optimize.curve_fit拟合基于mpmath.zeta的函数到实验数据时，不传bounds参数代码可正常执行并返回合理结果；但传入bounds后触发错误：

ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

追踪错误到least_squares源码第836行：

if not np.all(np.isfinite(f0)):
    raise ValueError("Residuals are not finite in the initial point.")

原因是f0为numpy数组，但元素类型是mpmath.ctx_mp_python.mpf，np.isfinite无法处理该类型。修改源码将f0转为float可解决问题，但不想改动源码。

疑问与解答 #

1. _wrap_func(func, xdata, ydata, transform)的实际作用是什么？ #

_wrap_func是curve_fit内部的工具函数，核心作用是将自定义拟合函数转换为优化器所需的残差计算函数：

• 无transform时，直接返回f(xdata, *params) - ydata（拟合值与实验值的原始残差）
• 若transform是一维数组，返回加权残差（transform作为权重数组）
• 若transform是二维数组，实现基于协方差矩阵的加权残差计算（对应卡方最小化的加权形式）
  本质是把拟合逻辑包装成优化器（least_squares）能直接用于最小化的目标函数。

2. 为什么f0包含mpmath类型对象而非浮点数？ #

你的拟合函数f中用np.frompyfunc(zeta,2,1)将mpmath.zeta转换为numpy兼容函数，但frompyfunc不会自动把mpmath返回的mpf类型转为Python原生float或numpy float。zeta计算结果为mpf对象，后续算术运算会保留该类型，最终返回的数组元素均为mpf类型，导致f0（初始参数下的残差数组）成为含mpf元素的numpy数组。

3. 能否不修改least_squares源码解决该问题？ #

有两种可靠方案：

• 方案1：在拟合函数内显式转换类型
  计算完zeta1/zeta2/zeta3后，将结果转为numpy float类型：

zeta1 = npzeta(0.5, 0.5 + (lb**2)/(l**2)).astype(np.float64)
zeta2 = npzeta(0.5, 0.5 + (lb**2)/(lfi**2)).astype(np.float64)
zeta3 = npzeta(0.5, 0.5 + 4*(lb**2)/(lso**2)+(lb**2)/(lfi**2)).astype(np.float64)

这样函数返回标准numpy float数组，np.isfinite可正常处理。

• 方案2：包装拟合函数统一转换类型
  不修改原拟合函数，在外层加包装函数转换结果类型：

def f_wrapper(B, l, lfi, lso, N):
    result = f(B, l, lfi, lso, N)
    return np.asarray(result, dtype=np.float64)

将f_wrapper传入curve_fit即可。

4. 为何仅传入bounds时才出现该错误？ #

不传bounds时，curve_fit默认使用lm（Levenberg-Marquardt）方法，该方法内部流程不会在初始阶段调用np.isfinite检查残差类型；传入bounds时，curve_fit自动切换到支持边界约束的默认方法trf（Trust Region Reflective），该方法在least_squares中会先检查初始残差是否为有限值，触发了np.all(np.isfinite(f0))代码，从而暴露了mpf类型不兼容的问题。

内容的提问来源于stack exchange，提问作者kuba_pol