Optimization algorithm the fastest descent method, Newton method, quasi-Newton method Python implementation

— — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — 2020.9.23 update — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — —

I’m rewriting BFGS of x to make it cleaner

def BFGS(x) :   # Quasi Newtonian method
    epsilon, h, maxiter = 10* * -5.10* * -5.10支那4
    Bk = np.eye(x.size)
    for iter1 in range(maxiter):
        grad = num_grad(x, h)
        if np.linalg.norm(grad) < epsilon:
            return x
        dk = -np.dot((np.linalg.inv(Bk)), grad)
        ak = linesearch(x, dk)
        x = x + dk*ak
        yk = num_grad(x, h) -grad
        sk = ak*dk
        if np.dot(yk, sk) > 0:
            Bs = np.dot(Bk,sk)
            ys = np.dot(yk,sk)
            sBs = np.dot(np.dot(sk,Bk),sk) 
            Bk = Bk - 1.0*Bs.reshape((n,1))*Bs/sBs + 1.0*yk.reshape((n,1))*yk/ys   
    return x
Copy the code

— — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — 2020.9.23 update — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — —

Numpy is the only library used, as long as the installation of this library should be able to run directly

import numpy as np

def f(x) :   # Objective function
    x1 = x[0]
    x2 = x[1]
    y = 100*((x2 - x1**2) * *2) + (x1-1) * *2
    return y

def num_grad(x, h) :     Gradient # o
    df = np.zeros(x.size)
    for i in range(x.size):
        x1, x2 = x.copy(), x.copy()  # use copy instead of assignment (=) because in Python = is an alias, not a copy (c/ C ++)
        x1[i] = x[i] - h
        x2[i] = x[i] + h
        y1, y2 = f(x1), f(x2)
        df[i] = (y2-y1)/(2*h)
    return df

def num_hess(x, h) :    # Find Hess matrix
    hess = np.zeros((x.size, x.size))
    for i in range(x.size):
        x1 = x.copy()
        x1[i] = x[i] - h
        df1 = num_grad(x1, h)
        x2 = x.copy()
        x2[i] = x[i] + h
        df2 = num_grad(x2, h)
        d2f = (df2 - df1) / (2 * h)   
        hess[i] = d2f     
    return hess
    
def linesearch(x, dk) :   Step # for long
    ak = 1
    for i in range(20):
        newf, oldf = f(x + ak * dk), f(x)
        if newf < oldf:
            return ak
        else:
            ak = ak / 4  Newf is smaller than oldf (ak=ak/2).
    return ak

def steepest(x) :   # Fastest descent method
    epsilon, h, maxiter = 10* * -5.10* * -5.10支那4
    for iter1 in range(maxiter):
        grad = num_grad(x, h)
        if np.linalg.norm(grad) < epsilon:
            return x
        dk = -grad
        ak = linesearch(x, dk)
        x = x + ak * dk
    return x

def newTonFuction(x) :   # Newton's method
    epsilon, h1, h2, maxiter = 10* * -5.10* * -5.10* * -5.10支那4
    for iter1 in range(maxiter):
        grad = num_grad(x, h1)
        if np.linalg.norm(grad) < epsilon:
            return x
        hess = num_hess(x, h2)
        dk = -np.dot((np.linalg.inv(hess)), grad)
        x = x + dk
    return x

def BFGS(x) :   # Quasi Newtonian method
    epsilon, h, maxiter = 10* * -5.10* * -5.10支那4
    Bk = np.eye(x.size)
    for iter1 in range(maxiter):
        grad = num_grad(x, h)
        if np.linalg.norm(grad) < epsilon:
            return x
        dk = -np.dot((np.linalg.inv(Bk)), grad)
        ak = linesearch(x, dk)
        x = x + dk*ak
        yk = num_grad(x, h) -grad
        sk = ak*dk
        if np.dot(yk.reshape(1, grad.shape[0]), sk) > 0:
            T0 = np.dot(Bk, SK) t1 = np.dot(T0. Shape [0], 1) sk.shape (1, SK.shape [0]) 0 0 = np.dot(T1, 0 Bk) temp1 = np.dot(np.dot(sk.reshape(1, sk.shape[0]), Bk), sk) tmp0 = np.dot(yk.reshape(yk.shape[0], 1), yk.reshape(1, yk.shape[0])) tmp1 = np.dot(yk.reshape(1, yk.shape[0]), sk) Bk = Bk - temp0 / temp1 + tmp0 / tmp1 '''
            # The second way to directly write the formula implementation
            Bk = Bk - np.dot(np.dot(np.dot(Bk, sk).reshape(sk.shape[0].1), sk.reshape(1, sk.shape[0])), Bk)/np.dot(np.dot(sk.reshape(1, sk.shape[0]), Bk), sk) + np.dot(yk.reshape(yk.shape[0].1), yk.reshape(1, yk.shape[0])) / np.dot(yk.reshape(1, yk.shape[0]), sk)
    return x


Array ([0.999960983973235, 0.999921911551354]) #
x0 = np.array([0.7.0.9])    # initial solution
x = steepest(x0)     Call the fastest descent method
print("The final solution vector of the fastest descent method:",x)
print("The final solution of the fastest descent method:",f(x))
print(' ')
x = newTonFuction(x0)     Call Newton's method
print("Final solution vector of Newton's method:", x)
print("Final solution of Newton's method:", f(x))
print(' ')
x = BFGS(x0)     Call the quasi-Newton method
print("Final solution vector of quasi Newtonian method:", x)
print("Final solution of quasi-Newton method:", f(x))
print(' ')
Copy the code

The results are as follows

The quasi Newtonian method feels troublesome, and I don’t have any ideas to change it for now, so let’s leave it at that