quasiNewtonSearchDirection()

quasiNewtonSearchDirection.py

# ------------------------------------------------------------------------
# The following function is implemented in Python by Professor Terje Haukaas
# at the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without any form of warranty.
# ------------------------------------------------------------------------

def quasiNewtonSearchDirection(x, oldPoint, gradient, oldGradient, HinvOld, type):

    # Initial declarations
    numVars = len(gradient)

    # If it is the first step, use steepest descent
    if np.linalg.norm(oldGradient) == 0.0:

        Hinv = np.multiply(np.identity(numVars), 0.001)
        searchDirection = np.multiply(gradient, -0.001)

        return [searchDirection, Hinv]

    else:

        deltaX = x - oldPoint
        deltaGrad = gradient - oldGradient

        if type == "DFP":

            deltaXOuterDeltaX = np.outer(deltaX, deltaX)
            deltaXDotDeltaGrad = deltaX.dot(deltaGrad)

            # Hinv = HinvOld + deltaXOuterDeltaX / deltaXDotDeltaGrad
            #      - (HinvOld * deltaGrad *outer* deltaGrad * HinvOld) / (deltaGrad * HinvOld * deltaGrad)
            Hinv = HinvOld + np.divide(deltaXOuterDeltaX , deltaXDotDeltaGrad) \
                   - np.divide(   (np.outer( HinvOld.dot(deltaGrad) , deltaGrad)).dot(HinvOld) , \
                                  (deltaGrad.dot(HinvOld)).dot(deltaGrad)   )

        elif type == "Broyden":

            # Hinv = HinvOld + ((deltaX - HinvOld * deltaGrad) *outer* deltaX * HinvOld)
            #      / (deltaX * HinvOld * deltaGrad)
            Hinv = HinvOld + np.divide(( np.outer( deltaX - (HinvOld.dot(deltaGrad)), deltaX).dot(HinvOld)) , \
                                       (deltaX.dot(HinvOld)).dot(deltaGrad)  )

        elif type == "SR1":

            # Hinv = HinvOld + ((deltaX - HinvOld * deltaGrad) *outer* (deltaX - HinvOld * deltaGrad))
            #      / ((deltaX - HinvOld * deltaGrad) * deltaGrad)
            Hinv = HinvOld + np.divide( np.outer( (deltaX - HinvOld.dot(deltaGrad)) , \
                 ((deltaX - HinvOld.dot(deltaGrad)))) , (deltaX - (HinvOld.dot(deltaGrad))).dot(deltaGrad) )

        elif type == "BFGS":

            Eye = np.identity(numVars)
            deltaXOuterDeltaGrad = np.outer(deltaX, deltaGrad)
            deltaXDotDeltaGrad = deltaX.dot(deltaGrad)
            deltaXOuterDeltaX = np.outer(deltaX, deltaX)

            # Hinv = (Eye - deltaXOuterDeltaGrad / deltaXDotDeltaGrad) * HinvOld * (Eye - deltaXOuterDeltaGrad
            #      / deltaXDotDeltaGrad) + deltaXOuterDeltaX / deltaXDotDeltaGrad
            Hinv = (Eye - (np.divide(deltaXOuterDeltaGrad , deltaXDotDeltaGrad)).dot(HinvOld)).dot(Eye \
                 - np.divide(deltaXOuterDeltaGrad , deltaXDotDeltaGrad)) + deltaXOuterDeltaX / deltaXDotDeltaGrad

        else:
            print("Wrong type given to quasi-Newton search direction algorithm")

    # Calculate the search direction -Hinv*grad
    searchDirection = np.multiply(-1, Hinv.dot(gradient))

    return [searchDirection, Hinv]