errorsk-4nQwBFp0LP8NpyK9WwwoT3BlbkFJyH23SKdFWR1P9Sr63sF8

#coding:UTF-8

from numpy import *
from function import *

def lbfgs(fun, gfun, x0):
    result = []#保留最终的结果
    maxk = 500#最大的迭代次数
    rho = 0.55
    sigma = 0.4
    
    H0 = eye(shape(x0)[0])
    
    #s和y用于保存最近m个，这里m取6
    s = []
    y = []
    m = 6
    
    k = 1
    gk = mat(gfun(x0))#计算梯度
    dk = -H0 * gk
    while (k < maxk):             
        n = 0
        mk = 0
        gk = mat(gfun(x0))#计算梯度
        while (n < 20):
            newf = fun(x0 + rho ** n * dk)
            oldf = fun(x0)
            if (newf < oldf + sigma * (rho ** n) * (gk.T * dk)[0, 0]):
                mk = n
                break
            n = n + 1
        
        #LBFGS校正
        x = x0 + rho ** mk * dk
        #print x
        
        #保留m个
        if k > m:
            s.pop(0)
            y.pop(0)
            
        #计算最新的
        sk = x - x0
        yk = gfun(x) - gk
        
        s.append(sk)
        y.append(yk)
        
        #two-loop的过程
        t = len(s)
        qk = gfun(x)
        a = []
        for i in xrange(t):
            alpha = (s[t - i - 1].T * qk) / (y[t - i - 1].T * s[t - i - 1])
            qk = qk - alpha[0, 0] * y[t - i - 1]
            a.append(alpha[0, 0])
        r = H0 * qk
            
        for i in xrange(t):
            beta = (y[i].T * r) / (y[i].T * s[i])
            r = r + s[i] * (a[t - i - 1] - beta[0, 0])

            
        if (yk.T * sk > 0):
            dk = -r            
        
        k = k + 1
        x0 = x
        result.append(fun(x0))
    
    return result

• 1.
• 2.
• 3.
• 4.
• 5.
• 6.
• 7.
• 8.
• 9.
• 10.
• 11.
• 12.
• 13.
• 14.
• 15.
• 16.
• 17.
• 18.
• 19.
• 20.
• 21.
• 22.
• 23.
• 24.
• 25.
• 26.
• 27.
• 28.
• 29.
• 30.
• 31.
• 32.
• 33.
• 34.
• 35.
• 36.
• 37.
• 38.
• 39.
• 40.
• 41.
• 42.
• 43.
• 44.
• 45.
• 46.
• 47.
• 48.
• 49.
• 50.
• 51.
• 52.
• 53.
• 54.
• 55.
• 56.
• 57.
• 58.
• 59.
• 60.
• 61.
• 62.
• 63.
• 64.
• 65.
• 66.
• 67.
• 68.
• 69.
• 70.
• 71.
• 72.