如何从Powell局部优化的盆地跳跃结果中获取解的误差

问题背景 #

我习惯用BFGS或LL-BFGS优化器，因为它们会提供用于误差计算的逆海森矩阵，但Nelder-Mead、Powell等优化器以及盆地跳跃算法（Basin Hopping）都不返回这个矩阵。现在我用Powell作为局部极小值求解器做盆地跳跃优化，想知道怎么获取解的误差范围（比如得到的9.349e3 ± ?）。

示例代码与输出 #

输入代码 #

from scipy.optimize import basinhopping
import numpy as np

# 假设fun是目标函数，input1、input2是函数参数
basinhopping(fun, minimizer_kwargs={"args":(input1,input2),"method" : 'Powell',"bounds":(0,np.inf)}, x0=1.0)

输出结果 #

message: ['requested number of basinhopping iterations completed successfully']
                    success: True
                        fun: 2.9445950465323385
                          x: [ 9.349e+03]
                        nit: 100
      minimization_failures: 0
                       nfev: 10141
 lowest_optimization_result: message: Optimization terminated successfully.
                             success: True
                              status: 0
                                 fun: 2.9445950465323385
                                   x: [ 9.349e+03]
                                 nit: 5
                               direc: [[ 1.805e+03]
                                       [ 2.962e+03]]
                                nfev: 175

得到的解为9.349e3，卡方值(chi²)为2.944，但不知道对应的误差范围。

我的尝试方案 #

我尝试用有限差分计算海森矩阵，再取逆矩阵的对角线元素作为误差，代码如下，不确定是否正确：

from statsmodels.base.optimizer import _fit_powell 
from statsmodels.tools.numdiff import approx_hess1
from numpy.linalg import inv
import numpy as np

solution = basinhopping(fun, minimizer_kwargs={"args":(input1,input2),"method" : 'Powell',"bounds":(0,np.inf)}, x0=1.0)
hessian = approx_hess1(solution.x, fun, args=(input1,input2))
print(np.diag(inv(hessian)))

疑问点 #

1. 上述代码是否正确？
2. 这种方法是否仅用于求微分？误差传播是否还需要协方差矩阵？
3. 已知误差传播公式：[df/dx df/dy][covariance matrix][df/dx df/dy]^T，如何计算参数本身的方差？逆海森矩阵是否就是方差？

解答 #

1. 代码修正与正确性 #

你的核心思路是对的，但代码有几处细节需要修正：

• 导入的是approx_hess1，但代码里误写为approx_hess；
• 最后一行缺少闭合括号；
• 应调用scipy.optimize.basinhopping而非内部函数_basinhopping。

修正后的代码：

from statsmodels.tools.numdiff import approx_hess1
from numpy.linalg import inv
import numpy as np
from scipy.optimize import basinhopping

solution = basinhopping(fun, minimizer_kwargs={"args":(input1,input2),"method" : 'Powell',"bounds":(0,np.inf)}, x0=1.0)
# 计算极小值点处的海森矩阵，传入边界避免试探点越界
hessian = approx_hess1(solution.x, fun, args=(input1, input2), bounds=(0, np.inf))
# 求逆得到协方差矩阵
cov_matrix = inv(hessian)
# 参数的标准差是协方差矩阵对角线元素的平方根
std_errors = np.sqrt(np.diag(cov_matrix))
print(f"解的误差范围：{solution.x[0]:.3e} ± {std_errors[0]:.3e}")

2. 协方差矩阵与误差的关系 #

• 参数本身的误差：当目标函数是卡方(chi²)拟合时，极小值点处的逆海森矩阵（除以2）就是参数的协方差矩阵。参数的方差是协方差矩阵的对角线元素，标准差则是方差的平方根——这就是你要的±值。
• 误差传播公式：这个公式是用于计算由参数导出的函数的误差时才需要的。比如你有函数f(x,y)，想求f的误差，就用该二次型计算f的方差。但如果只是求参数x、y本身的误差，直接用协方差矩阵的对角线元素即可。

3. 额外注意事项 #

• 有限差分计算海森矩阵时，approx_hess1会自动调整步长，但如果函数有噪声，可手动指定epsilon参数优化精度；
• 必须传入bounds参数，避免有限差分的试探点超出边界导致计算错误；
• 务必使用盆地跳跃返回的全局最优解solution.x计算海森矩阵，不要用中间局部极小值点。

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