python Scipy Minimize,逆hessian可以用来消除具有类似“乐趣”的多个解决方案吗?

zdwk9cvp  于 2023-04-10  发布在  Python
关注(0)|答案(1)|浏览(81)

假设你有一个最小化问题:

data_set_1=[[24.714, 24.713, 24.605, 24.607], [25.708, 25.59, 25.557, 25.753], [22.713, 22.654, 22.655, 22.812], [24.231, 24.233, 24.289, 24.109], [24.401, 24.385, 24.396, 24.408], [26.124, 26.107, 25.915, 25.847], [15.341, 15.242, 15.049, 14.847], [22.332, 22.282, 22.324, 22.314], [25.843, 25.835, 25.683, 25.983], [24.4, 24.402, 24.522, 24.409], [20.218, 20.23, 20.212, 20.254], [21.927, 21.955, 21.94, 21.917], [20.697, 20.705, 20.731, 20.692], [21.13, 21.134, 21.196, 21.266], [24.619, 24.55, 24.556, 24.784], [25.294, 25.29, 25.322, 25.337], [24.237, 24.124, 24.06, 24.143], [18.533, 18.603, 18.814, 19.229], [26.235, 26.183, 26.226, 26.137], [25.754, 25.72, 25.658, 25.639], [26.786, 26.716, 26.646, 26.668], [23.863, 23.838, 23.946, 24.034], [21.402, 21.395, 21.401, 21.317], [22.382, 22.386, 22.363, 22.244], [22.435, 22.454, 22.434, 22.579], [20.416, 20.444, 20.439, 20.583], [24.929, 24.943, 24.819, 24.916], [25.839, 25.811, 25.716, 25.752]]

data_set_2=[[0.786, 0.781, 0.766, 0.773], [0.78, 0.773, 0.782, 0.798], [0.983, 0.993, 0.991, 1.031], [0.74, 0.741, 0.756, 0.695], [0.614, 0.608, 0.642, 0.646], [0.659, 0.67, 0.666, 0.675], [0.488, 0.456, 0.438, 0.571], [0.99, 0.979, 0.97, 0.959], [0.806, 0.815, 0.801, 0.818], [0.772, 0.772, 0.756, 0.752], [0.364, 0.373, 0.357, 0.394], [0.735, 0.727, 0.729, 0.705], [0.489, 0.484, 0.495, 0.462], [0.803, 0.806, 0.817, 0.831], [1.018, 1.02, 1.021, 0.993], [0.589, 0.606, 0.599, 0.622], [0.61, 0.612, 0.591, 0.615], [0.973, 0.955, 0.94, 0.956], [0.628, 0.634, 0.643, 0.669], [0.64, 0.656, 0.637, 0.619], [0.732, 0.738, 0.742, 0.743], [0.872, 0.865, 0.859, 0.838], [0.199, 0.205, 0.21, 0.221], [0.783, 0.771, 0.775, 0.727], [1.069, 1.075, 1.064, 1.1], [0.707, 0.708, 0.705, 0.741], [1.07, 1.061, 1.058, 1.034], [0.911, 0.908, 0.91, 0.888]]

def get_populations(initial,io):
 k,k1,x,y=initial[0],initial[1],initial[2],initial[3]
 kx,k1x=k*x,k1*x
 ky,k1y=k*y,k1*y
 kxy,k1xy=k*x*y,k1*x*y
 pF=(np.sqrt((k*k1)**2+(8*io*k*k1)+(8*io*k*k1**2))-(k*k1))/(4*(1+k1))
 pC=(((4*io)+(4*io*k1))/(4*(1+k1)**2))-(pF/(1+k1))
 pO=k1*pC
 pF2=(np.sqrt((kx*k1x)**2+(8*io*kx*k1x)+(8*io*kx*k1x**2))-(kx*k1x))/(4*(1+k1x))
 pC2=(((4*io)+(4*io*k1x))/(4*(1+k1x)**2))-(pF2/(1+k1x))
 pO2=k1x*pC2
 pF3=(np.sqrt((ky*k1y)**2+(8*io*ky*k1y)+(8*io*ky*k1y**2))-(ky*k1y))/(4*(1+k1y))
 pC3=(((4*io)+(4*io*k1y))/(4*(1+k1y)**2))-(pF3/(1+k1y))
 pO3=k1y*pC3
 pF4=(np.sqrt((kxy*k1xy)**2+(8*io*kxy*k1xy)+(8*io*kxy*k1xy**2))-(kxy*k1xy))/(4*(1+k1xy))
 pC4=(((4*io)+(4*io*k1xy))/(4*(1+k1xy)**2))-(pF4/(1+k1xy))
 pO4=k1xy*pC4
 local_chi2=0
 for data1,data2 in zip(data_set_1,data_set_2):
   populations=np.array([[pF,pO,pC],[pF2,pO2,pC2],[pF3,pO3,pC3],[pF4,pO4,pC4]])
   least_squared_fit_1=lsmr(populations/io,np.array([data1])/4*800,maxiter=10)
   least_squared_fit_2=lsmr(populations/io,np.array([data2])*800,maxiter=10)
   local_chi2+=least_squared_fit_1[3]**2+least_squared_fit_2[3]**2
 return local_chi2

global_parameter_solution=minimize(get_populations,args=io, x0=np.array([2e3,2e-3,9,2e1]),bounds=np.array([[0,np.inf],[0,np.inf],[0,np.inf],[0,np.inf]]),options={'maxiter':100000})

不幸的是,这个问题的前景是相当平坦的,所以有多种解决方案。我正在寻找2个特定的解决方案,使用2个初始猜测。
使用5e 2、2 e-2、7、3el的溶液

message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
  success: True
   status: 0
      fun: 1568.1943235398385
        x: [ 5.000e+02  1.226e-02  6.638e+00  3.039e+01]
      nit: 8
      jac: [-2.683e-03 -1.183e-01 -2.933e-03 -1.796e-03]
     nfev: 50
     njev: 10
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>

的逆hessian对角线

[0.99964765 0.00841953 1.41597902 3.25752296]

使用2 e3、2 e-3、9、2 el的溶液2

message: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH
  success: True
   status: 0
      fun: 1568.1943236577647
        x: [ 2.000e+03  9.445e-03  6.398e+00  2.838e+01]
      nit: 24
      jac: [-1.137e-04  1.100e-01 -1.000e-03 -1.523e-03]
     nfev: 370
     njev: 74
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>

带逆粗麻布的

[1.00008883e+00 2.07424141e-03 5.98621521e+01 6.83255242e+01]

两个解的“fun”或残差之和基本相同。然而,一个解的方差(逆hessian的对角线)明显比另一个差。我可以根据其中一个参数的方差非常高而另一个解是有效的来排除解2吗?

kq0g1dla

kq0g1dla1#

对不起,这里没有免费的午餐。除非你能提供更多的细节,否则两种解决方案都是同样可行的。例如,一种解决方案在物理上是首选的,或者你可以对参数设置界限或约束。
除了第一个参数外,这两个解看起来并没有太大的不同,也许你应该画出get_populations作为第一个参数的函数是如何变化的?
〈免责声明-我没有详细查看代码〉

相关问题