Python：牛顿法、黑森法、雅可比法

Question

您好，我正在尝试使用牛顿法来最小化一个函数，但是当我运行代码时我不断收到这个错误，我不知道为什么。任何帮助深表感谢。谢谢！

错误：

ValueError: shapes (2,1) and (2,1) not aligned: 1 (dim 1) != 2 (dim 0)  

代码：

import sympy as sy
from sympy import symbols
import numpy as np
from numpy import linalg as la
from scipy.optimize import minimize
a1=0.3
a2=0.6
a3=0.2
b1=5
b2=26
b3=3
c1=40
c2=1
c3=10
h=0.000001

def TutFunc(x):
    x=np.empty((2,1))
    u = x[0][0] - 0.8
    v = x[1][0] - ((a1+(a2*u**2))*(1-u)**0.5-(a3*u))
    alpha = -b1+b2*u**2*(1-u)**0.5+b3*u
    beta = c1*v**2*(1-c2*v)/(1+c3*u**2)
    y= alpha*np.exp(-beta)
    return y

def JacobianFun(x):
    x=np.empty((2,1))
    Jx1 = (TutFunc(x+[[h],[0]]) - TutFunc(x-[[h],[0]]))/(2*h)
    Jx2 = (TutFunc(x+[[0],[h]]) - TutFunc(x-[[0],[h]]))/(2*h)
    jaco = np.array([[Jx1],[Jx2]])
    return jaco

def HessianFun(x): 
    x=np.empty((2,1))
    Hx1x1 = (TutFunc(x+[[h],[0]]) - 2*TutFunc(x) + TutFunc(x-[[h],[0]]))/h**2
    Hx1x2 = (TutFunc(x+[[h],[h]]) - TutFunc(x+[[h],[-h]]) - TutFunc(x+[[-h],[h]]) + TutFunc(x-[[h],[h]]))/(4*h**2)
    Hx2x1 = Hx1x2
    Hx2x2 = (TutFunc(x+[[0],[h]]) - 2*TutFunc(x) + TutFunc(x-[[0],[h]]))/h**2
    Hess = np.array([[Hx1x1, Hx1x2],[ Hx2x1, Hx2x2]])
    return Hess

x0=([0.7, 0.3]
x=minimize(TutFunc,x0,method= 'Newton-CG', jac=JacobianFun, hess=HessianFun)

Answer 1

我已尝试优化您的输出。您需要更改 jacobian 和 hessian 函数。我已经改变了你需要自己做的 jacobian, hessian。

See the documentation here.

根据文档：jac(x) -> array_like, shape (n,) 这意味着 jacobian 函数采用 x 这是一个 ndarray 并返回具有 (n,0) 维度的 array。在你的情况下（2,0）。另一方面，您有 (2,2) 所以我一个一个地在它们之间切换以实现优化。如果您愿意，可以浏览文档。第二件事是 hess 和 hessp 是两个不同的函数。

hess : hess(x, *args) -> {LinearOperator, spmatrix, array}, (n, n)

接受参数并返回一些我真的不知道的东西。而通过查看您的代码，我认为您正在考虑使用：

hessp: hessp(x, p, *args) -> ndarray 形状 (n,)

接受参数并返回形状（n，）的ndarray

为简单起见，我禁用了 hessp

另外，我不明白为什么您将数组 X 初始化为一个空数组。

代码如下：

import sympy as sy
from sympy import symbols
import numpy as np
from numpy import linalg as la
from scipy.optimize import minimize
a1=0.3
a2=0.6
a3=0.2
b1=5
b2=26
b3=3
c1=40
c2=1
c3=10
h=0.000001
flag = 0

def TutFunc(x):
    u = x[0] - 0.8
    v = x[1] - ((a1+(a2*u**2))*(1-u)**0.5-(a3*u))
    alpha = -b1+b2*u**2*(1-u)**0.5+b3*u
    beta = c1*v**2*(1-c2*v)/(1+c3*u**2)
    y= alpha*np.exp(-beta)
    return y

def JacobianFun(x):
    global flag
    Jx1 = (TutFunc(x+[[h],[0]]) - TutFunc(x-[[h],[0]]))/(2*h)
    Jx2 = (TutFunc(x+[[0],[h]]) - TutFunc(x-[[0],[h]]))/(2*h)
    jaco = np.array([Jx1 , Jx2])
    if flag == 0:
        flag = 1
        return jaco[0]
    else:
        flag = 0
        return jaco[1]

def HessianFun(x): 
    x=x
    Hx1x1 = (TutFunc(x+[[h],[0]]) - 2*TutFunc(x) + TutFunc(x-[[h],[0]]))/h**2
    Hx1x2 = (TutFunc(x+[[h],[h]]) - TutFunc(x+[[h],[-h]]) - TutFunc(x+[[-h],[h]]) + TutFunc(x-[[h],[h]]))/(4*h**2)
    Hx2x1 = Hx1x2
    Hx2x2 = (TutFunc(x+[[0],[h]]) - 2*TutFunc(x) + TutFunc(x-[[0],[h]]))/h**2
    Hess = np.array([[Hx1x1, Hx1x2],[ Hx2x1, Hx2x2]])
    if flag == 0:
        flag = 1
        return Hess[0]
    else:
        flag = 0
        return Hess[1]


x0= np.array([0.7, 0.3])
x0.shape
x=minimize(TutFunc,x0,method= 'Newton-CG', jac=JacobianFun, hessp=None)
print(x)

当 hess 功能被禁用时 x 的输出

    jac: array([ 2.64884244, -2.89355671])
message: 'Optimization terminated successfully.'
   nfev: 2
   nhev: 0
    nit: 1
   njev: 6
 status: 0
success: True
      x: array([0.69999996, 0.30000004])