Calcula da hessiana e atualização - Broyden-Fletcher-Goldfarb-Shanno (BFGS)

gistfile1.txt

import numpy as np
import scipy.linalg as la
import numpy.linalg as ln

#Busca valor de alpha
def LineSearch(f,g,x,p):
    a,b = 1 - 2/(1 + np.sqrt(5)), 2/(1+np.sqrt(5))
    alpha = 1.0
    while (f(x+alpha*p) > f(x)) + a * alpha*np.dot(g(x),p) :
        alpha *= b
    return alpha
    
#Desenvolvimento do método BFGS
def BFGS(f,g,x,eps,maxiter=10, alfa=0.1,inv_hessian=None):
    
    #Se não houver valor para Hessiana
    if inv_hessian is None:
        H = scipy.eye(g(x).shape[0])
        print('Hessiana: ', H)
       
    #Inicializa a variável
    i = 0
    
    #Lopping do método onde interrempo a iteração ao atingir os critérios de parada
    while ln.norm(g(x),2) >= eps and i < maxiter:
        #Direção
        p = -np.dot(H, g(x))
        #Tamanho do passo
        alpha = LineSearch(f,g,x,p)
        s = alpha * p 
        #Calcula gradiente em x
        grad_m1 = scipy.zeros(g(x).shape)
        grad_diff = g(x) - grad_m1
        #Atualizando o valor da Hessiana
        ys = np.inner(grad_diff, s)
        Hy = np.dot(H, grad_diff)
        yHy = np.inner(grad_diff, Hy)
        print('Passo 1: ', (ys + yHy) )
        print('Passo 2: ', np.outer(s, s))
        print('Passo 3: ', ys ** 2)
        print('Juntando 1: ', (ys + yHy) * np.outer(s, s))
        print('Juntando 3: ', (ys + yHy) * np.outer(s, s)/ ys ** 2)
        H += (ys + yHy) * np.outer(s, s) / ys ** 2
        print('atualiza H: ', H)
        H -= (np.outer(Hy, s) + np.outer(s, Hy)) / ys        
#         p = -np.linalg.solve(H,g(x))
#         print('Direcao: ', p)
#         alpha = LineSearch(f,g,x,p)
#         print('Alfa:', alpha)
#         s = alpha * p 
#         print('Valor gradiente :', s)
#         y = g(x+alpha*s) - g(x)
#         x = x + s
#         H = H + np.outer(y,y)/np.inner(y,s) - (H@np.outer(s,s)@H.T)/(s.T@H@s)
#         print('Atualiza Hessiana: ', H)
        print('{0:5.0f}'.format(i),'\t','\t','{0:12.10f}'.format(f(x)),'\t','{0:10.8e}'.format(np.linalg.norm(g(x))))
        print('-------------------------------------------------------------')
        i += 1
        
    return x, i
    
#Função a ser otimizada e minimizada
def f(x):
    return x[0]**2 + 2*x[1]**2 + (-2*x[0]*x[1] - 2*x[1])
# Encontra a matriz da função a partir da derivada primeira
def f1(x):
    return np.array([2 * x[0] - 2, 4*x[1] - 4])

x0 = np.array([1,2])
eps = 1e-8
x, iteracao = BFGS(f,f1,x0,eps)
print('Iteracao = ', iteracao)
print('valor de x: ', x)

def LineSearch(f,g,x,p):
a,b = 1 - 2/(1 + np.sqrt(5)), 2/(1+np.sqrt(5))
alpha = 1.0
while (f(x+alphap) > f(x) + a * alphanp.dot(g(x),p)) : #
alpha *= b
return alpha