Python: метод Ньютона, Гессе, метод Якоби

Привет, я пытаюсь использовать метод Ньютона для минимизации функции, но я продолжаю получать эту ошибку при запуске кода, и я не знаю, почему. Любая помощь горячо приветствуется. Спасибо!

Ошибка:

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)

python-3.x python hessian newtons-method
person Nu2prog    schedule 24.08.2018    source источник
comment
что означает x0=([0.7, 0.3]? это неполное выражение   -  person Evgeny    schedule 24.08.2018
comment
Я думаю, что это должно быть x0 = np.array([0.7,0.3])   -  person NoorJafri    schedule 24.08.2018
comment
@Syed Noor Ali Jafri - возможно, это решает проблему! Дополнительные подсказки для OP на stackoverflow.com/questions/39608421/   -  person Evgeny    schedule 24.08.2018
comment
x0 = np.array ([0.7, 0.3]) дает размерность 2,0, а x0 = np.array ([[0.7], [0.3]]) дает размерность 2,1. OP может проверить размер, используя x0.shape. Это может дать вам ключ к разгадке размерной ошибки.   -  person NoorJafri    schedule 24.08.2018
comment
@SyedNoorAliJafri, который не решает эту проблему, возникает та же ошибка   -  person Nu2prog    schedule 24.08.2018
comment
заглядывая глубже в код - вы повторяете def JacobianFun(x): x=np.empty((2,1)), который эффективно стирает наш аргумент в каждой функции, я думаю, что это неправильное поведение.   -  person Evgeny    schedule 24.08.2018
comment
@EPo затем он использует его для инициализации. Я не знаю, что это такое, но ошибка генерируется где-то еще.   -  person NoorJafri    schedule 24.08.2018
comment
@Syed Noor Ali Jafri - нет смысла передавать x функции, а затем перезаписывать ее. Я думаю, что ОП должен указать, что ему нужно: повторно запустить матрицу функций или вычислить значение в какой-то точке x, я думаю, что код становится смесью двух способов, и неясно, какой тип JacobianFun возвращает - вектор с плавающей запятой вектора двух функций.   -  person Evgeny    schedule 24.08.2018
comment
При отладке функция минимизации использует numpy.dot(ri, ri), что приближает нас к ошибке. Проверьте эту тему: stackoverflow.com/questions/28028991/< /а>   -  person NoorJafri    schedule 24.08.2018
comment
@ Nu2prog Nu2prog, вы уверены, что у вас есть функция jacobian, верно? Пожалуйста, проверьте мое решение, хотя оптимизация достигнута: p   -  person NoorJafri    schedule 24.08.2018

Ответы (1)


arrow_upward
0
arrow_downward

Я попытался оптимизировать ваш вывод. Вам нужно изменить функции Якобиана и Гессе. Я переделал якобиан, гессиан, который вы должны сделать сами.

См. документацию здесь.

Согласно документации: jac(x) -> array_like, shape (n,) Это означает, что функция jacobian принимает x, который является ndarray, и возвращает array с (n,0 ) измерение. В вашем случае (2,0). С другой стороны, у вас было (2,2), поэтому я переключался между ними обоими один за другим, чтобы добиться оптимизации. Вы можете просмотреть документацию, если хотите. Во-вторых, hess и hessp — это две разные функции.

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

принимает аргументы и возвращает то, о чем я действительно понятия не имею. Принимая во внимание, что, увидев ваш код, я думаю, вы думали об использовании:

hessp: hessp(x, p, *args) -> форма ndarray (n,)

Который принимает аргументы и возвращает ndarray формы (n,)

Для простоты я отключил hesp.

Плюс я не вижу смысла, почему вы инициализировали массив 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)

ВЫВОД для x, когда функция Hess отключена

    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])
person NoorJafri    schedule 24.08.2018
comment
@Nu2prog Nu2prog ты пробовал решение? - person NoorJafri; 27.08.2018