大数值场景下使用lmfit进行最小二乘优化溢出报错如何解决

问题根因 #

当输入的x为绝对值较大的负数、且参数b较小时，-x/b会得到极大的正数，超出float64类型下np.exp()的可计算范围（上限约为709.78），最终生成inf值导致拟合报错。

可直接落地的修改方案 #

1. 指数项数值截断 #

在计算指数前对输入值做范围裁剪，避免溢出，同时不影响正常区间的计算结果。

2. 残差异常值兜底 #

对残差中的NaN/Inf值做替换，返回超大惩罚值，引导拟合算法避开会产生数值异常的参数组合。

3. 优化参数初始值 #

结合输入profitlist的取值范围调整b的初始值，减少迭代初期就出现数值溢出的概率。

修改后完整代码 #

import numpy as np
from lmfit import Parameters, minimize, fit_report

class LeastSquares:
    def __init__(self, profitList, utilityList):
        self.profitList = np.asarray(profitList)
        self.utilityList = np.asanyarray(utilityList)

    def function(self, params, x):
        a = params["a"]
        b = params["b"]
        c = params["c"]
        # 裁剪指数输入范围，避免溢出
        exp_input = np.clip(-x / b, -709, 709)
        return a * np.exp(exp_input) + c

    def residual(self, params, x, y):
        residual_val = (y - self.function(params, x))**2
        # 异常值兜底，替换NaN/Inf为超大惩罚值
        return np.nan_to_num(residual_val, nan=1e18, posinf=1e18, neginf=1e18)

    def setParameters(self, a_start, b_start, c_start):
        parameters = Parameters()
        parameters.add(name="a", value=a_start, min=None, max=0, vary=True)
        parameters.add(name="b", value=b_start, vary=True, min=0.1, max=None)
        parameters.add(name="c", value=c_start, vary=True)
        return parameters 

    def startOptimalization(self):
        # 自动计算适配输入范围的b初始值，避免初始迭代就溢出
        x_range = np.max(self.profitList) - np.min(self.profitList)
        init_b = max(1, x_range / 10)
        parameters = self.setParameters(-1, init_b, 1)    
        result = minimize(self.residual, parameters, args=(self.profitList, self.utilityList), method="leastsq")
        result.params.pretty_print()
        print(fit_report(result))
        print("SSE")
        print(np.sum(result.residual))

# 测试示例
if __name__ == "__main__":
    profitlist = [-10000, 100.00, 1000.00, 100000.00, 1000000.00]
    utilitylist = [0, 0.2, 0.4, 0.6, 1]
    ls = LeastSquares(profitlist, utilitylist)
    ls.startOptimalization()

补充说明 #

• 裁剪阈值709是float64类型下np.exp()能计算的最大输入值，裁剪后只会对超出合理范围的指数做饱和处理，不会影响正常拟合效果
• 残差的惩罚值可以根据实际业务的数值范围调整，只要远大于正常残差即可

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