準ニュートン法（BFGS)を使った最小値の探索にて、線形探索の中で無限ループから抜け出せません。

Question

前提・実現したいこと

バナナ関数(rosenblock)を最適化したいです。
私のコードは、他の簡単な関数の最適化ができたのですが、どうしてもバナナ関数ができません。
BFGSの準ニュートン法を使って最適化をしています。

発生している問題・エラーメッセージ

BFGSを使い、最初の最適値aを求める線形探索で無限ループに入ってしまい、計算することができません。

該当のソースコード

def rosenbrock(x):
    global FUNCCALLS
    FUNCCALLS += 1

    funcValue = (100*((x[1]-x[0]****2)****2) + (1-x[0])**2)

    gradient = np.matrix([(400*(x[0]**3)-400*x[0]*x[1]+2*x[0]-2), (200*(x[1]-x[0]**2))]).T

    return funcValue, gradient

#この関数を、optimize に呼び出します。

def optimize(func, x0, epsilon_g,n):

    H=np.identity(n)
    I=np.identity(n)
    x=np.matrix([x0])
    x_star=np.matrix([x0])
    a_star=0.1
    g_pre=None
    s_pre=None
    k=0


    while True:

        l=x[0,0]
        m=x[0,1]
        t=[l,m]
        J,g=func(t)

        if g_pre is None:
            row= ((-g).T*(-x).T)/(-g.T*(-g))
            H=np.multiply(H,row)
            print(row)


        p=-np.linalg.inv(H)*g



        sumJ=abs(np.sum(g))
        if np.linalg.norm(g, ord = np.inf) <= 1e-6:
            print("find x_star")
            break
        print(x,p)


        a_star=bracket(x,p,func)
        print(a_star)
        print(p)

        s=np.multiply(a_star,p)
        x+=a_star*s.T


        x_star=np.vstack((x_star,np.array(x)))

        if g_pre is not None:
            y=g-g_pre
            dx=x_star[-1]-x_star[-2]
            sk=dx.T
            print(sk)
            print(y)
            print(H)
            '''dfpg'''
            H=H+(y*y.T)/(y.T*s_pre)-(H*s_pre*s_pre.T*H)/(s_pre.T*H*s_pre)
            '''dfp'''

        g_pre=g
        s_pre=(a_star*s.T).T
        k+=1
        print("loop" +str(k))

    xopt=x
    fopt=J
    outputs=k

    return xopt, fopt, outputs

#アルファの値を見つけるために、bracketing （線形探索）のコードを呼び出します。

def bracket(t,p,func):
    print("in the blacket")
    print(t[0,0])
    print(t[0,1])
    a_max=10 # alpha max
    pre_a=0
    a=3 # define arrays
    print(a)
    print(pre_a)
  
    i=1 #counter
    rho=5 # the number of step size
    mu_1=0.1 #num of mu1
    mu_2=0.9 #num of mu2

    J_0,g_0=utils.phi(0,t[0,0],t[0,1],p,func) #find initial value,
    dp_0=utils.dphi(0,t[0,0],t[0,1],p,func)
    print(dp_0)
    print(t[0,0],t[0,1],p)

    while i<=a_max:
        J,g=utils.phi(a,t[0,0],t[0,1],p,func) #(0,x[0],x[1],p)
        p=p
        print("what was the inside?")
        print(J)
        print(mu_1*a*dp_0)

        if J>(J_0+mu_1*a*dp_0):
            print("pinpoint I")
            print(pre_a)
            print(a)
            a_star = utils2.pinpoint(pre_a,a,mu_1,mu_2,t,p,func)
            print("number of iterative is "+str(i))
            return a_star


        elif ((i>1) and (J>utils.phi(pre_a,t[0,0],t[0,1],p,func)[0])):
            print("pinpoint II")
            a_star = utils2.pinpoint(pre_a,a,mu_1,mu_2,t,p,func)
            print("number of iterative is "+str(i))

            return a_star


        dp=utils.dphi(a,t[0,0],t[0,1],p,func)
        if abs(dp)<=(-mu_2*dp_0):
            a_star=a
            print("number of iterative is "+str(i))

            return a_star

        elif dp>=0:
            print("pinpoint III")
            a_star = utils2.pinpoint(pre_a,a,mu_1,mu_2,t,p,func)
            print("number of iterative is "+str(i))
            return a_star

        else:
            pre_a=a
            a=(rho*a) #add new value of a, a[i+1]=rho*a[i]
            print(a)


        i=I+1

#最小値が見つかりそうな場所に入ったら、pinpointにて値を出します。

def pinpoint(l,h,mu_1,mu_2,x,p,func):
    print("in the pinpoint")
    global high
    global low
    high=h
    low=l
    print(low)
    print(high)
    a_max=10
    a=0 # define the range of arrays
    j=0
    J_0,g_0=utils.phi(0,x[0,0],x[0,1],p,func) #calculate initial value, J,g=endurance(0)

    while True:
        print("in the pinpoint loop" +str(j))
        print(x)
        print(low)
        print(high)
        J_low,g_low=utils.phi(low,x[0,0],x[0,1],p,func)
        J_high,g_high=utils.phi(high,x[0,0],x[0,1],p,func)
        d_phi=utils.dphi(low,x[0,0],x[0,1],p,func)
        print("Are you sure to calculate?")
        print(J_low,J_high,d_phi)
        print("actual is")
        print(utils.phi(low,x[0,0],x[0,1],p,func),utils.phi(high,x[0,0],x[0,1],p,func),utils.dphi(low,x[0,0],x[0,1],p,func))
        num=2*low*(J_high-J_low)+d_phi*(low**2-high**2)
        den=2*((J_high-J_low)+d_phi*(low-high))
        a=(num)/(den)
        print("this is a")
        print(a)
        print(num, den )

        J_j,g_j=utils.phi(a,x[0,0],x[0,1],p,func)
        dp_0=utils.dphi(0,x[0,0],x[0,1],p,func)
        print("whats the p so far:1")
        print(p)
        print("J is")
        print(J_j)
        print ("the smaller one is" )
        print(J_0+mu_1*a*dp_0)
        if J_j>(J_0+mu_1*a*dp_0):
            print("I am here1")
            print(high)
            print(low)
            print(a)
            high=a
            print(high)
            print(a)

        elif J_j>J_low:
            print("I am here2")
            high=a


        else:
            print("whats the high so far:2")
            print(high)
            print("I am here3")
            J_dash_j=utils.dphi(a,x[0,0],x[0,1],p,func)
            if abs(J_dash_j)<=-mu_2*dp_0:
                print("I am here4")
                a_star=a
                print("total num of j is "+str(j))
                return a_star
            elif J_dash_j*(high-low)>=0:
                print("I am here5")
                high=low
                print("whats the high so far :unexpect")
                print(high)
            print("whe the a doesnt moove")
            print(a)
            low=a

            print("I overrride")
            print(low)
            print(high)
            print(low)
            print(high)

        j=j+1

上記のピンポイントのループから抜け出せなくなり、困ってます。
助けて頂けないでしょうか。

#よび出し関数です。

testFuncs = [rosenbrock]



for _, func in enumerate(testFuncs):

	# set the seed so every one is tested at the same points
    np.random.seed(0)

    trialScore = 0
    for _ in range(nTrails):
        FUNCCALLS = 0


        # get a random starting point
        x0 = np.random.normal(size=2)*3


        # call the optimizer
        xopt, fopt, output = uncon(func, x0, 1e-6)


        l=xopt[0,0]
        m=xopt[0,1]
        xopt=[l,m]
        _, gopt =  rosenbrock(np.array(xopt))
        #_, gopt =  matyas(np.array(xopt))
        #_, gopt =  paraboloid(np.array(xopt))
        if (not (np.linalg.norm(gopt, ord = np.inf) <= 1e-6)):
         
            
            continue


        print(xopt)
        print('fopt ', fopt, FUNCCALLS)

print(output['alias'], score)

python

試したこと

それぞれの初期値を変えたり、初期ヘッセ行列を等倍したり、値を変えてみたのですが、
どうしてもp（方向）の値が小さくなってしまい、最小値a_starを見つけるコードから抜け出せません。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Comments

[t_obara]

質問は修正可能ですので、pythonのコードは```で囲ってインデントもわかるようにみやすくしてください。

[alpha_high]

コメントありがとうございます。
また、サイトの利用方法について教えて頂き、ありがとうございます。すぐに修正します。