[MMTK] Problems adding configuration to Universe Copy

vanitha@cs.wisc.edu vanitha at cs.wisc.edu
Wed Feb 8 19:18:47 CET 2006


I've written an implementation of the Limited Memory BFGS minimization
algorithm in Python within the MMTK framework. I am running into a problem
when I create a universe, figure out a "search direction" and then I add
this to a copy of my universe. I get the following problem when I try to
do this. (see traceback and code below). How can I resolve this? I can
repeat code all over again, but that really won't help. I'm trying to see
how I can avoid computing the search direction twice, once with the
original universe and then with a copy of the universe.

Any help will be greatly appreciated!

Thanks,

- Vanitha


Traceback (most recent call last):
  File "testlbfgs.py", line 13, in ?
    Minimization.LBFGSMinimize(universe, 7, 1000, 0.0001)
  File
"/afs/cs.wisc.edu/u/v/a/vanitha/private/MMTK-2.5.11/Examples/Minimization/Test3/Minimization.py",
line 105, in LBFGSMinimize
    alfa = backtrackingLineSearch( universeCopy, initialAlfa, direction,
maxIter)
  File
"/afs/cs.wisc.edu/u/v/a/vanitha/private/MMTK-2.5.11/Examples/Minimization/Test3/Minimization.py",
line 147, in backtrackingLineSearch
    Dphi0 = gradient.dotProduct(direction)
  File
"/u/v/a/vanitha/lib/python2.4/site-packages/MMTK/ParticleProperties.py",
line 286, in dotProduct
    raise ValueError('Variables are for different universes')
ValueError: Variables are for different universes



def LBFGSMinimize( universe, m, maxIter, tolerance):

        # Iteration counter
        k = 0

        # The current configuration of the universe is the initial guess
        energy, gradient = universe.energyAndGradients()
        configuration = universe.copyConfiguration()

        while 1:
                print energy, gradient.norm()

                ###############################################################
                #               Termination Criteria
                ###############################################################
                if( k >= maxIter):
                        break

                if( gradient.norm() <= tolerance):
                        break

                ################################################################
                # Choose the initial approximation to the inverse of the
Hessian
                # This is given by Hk0 = gamma*I
                # The following steps compute the value of gamma
                ################################################################
                if( k == 0):
                        # Initialize the rho, s and y vectors
                        rho = []
                        s = []
                        y = []
                        gamma = 1
                else:
                        skMinusm = s[len(s)-1]
                        ykMinusm = y[len(y)-1]

                        # Skip the update if negative curvature
                        # was detected in the previous step.
                        #if( skMinusm.dotProduct(ykMinusm) > 0):

                        gammaTop = skMinusm.dotProduct(ykMinusm)
                        gammaBottom = ykMinusm.dotProduct(ykMinusm)
                        gamma = gammaTop/gammaBottom

                ####################################################
                #       Get the descent direction
                ####################################################
                universeCopy = copy(universe)
                energy, gradient = universeCopy.energyAndGradients()
                q = gradient
                alpha = []

                # Start from len(s)-1 in steps of -1 until 0
                for i in range(len(s)-1, -1, -1):
                        alpha.append(rho[i]*(s[i].dotProduct(q)))
                        q = q - alpha[len(s)-1-i]*y[i];

                alpha.reverse()
                r = gamma*q

        while 1:
                print energy, gradient.norm()

                ###############################################################
                #               Termination Criteria
                ###############################################################
                if( k >= maxIter):
                        break

                if( gradient.norm() <= tolerance):
                        break

                ################################################################
                # Choose the initial approximation to the inverse of the
Hessian
                # This is given by Hk0 = gamma*I
                # The following steps compute the value of gamma
                ################################################################
                if( k == 0):
                        # Initialize the rho, s and y vectors
                        rho = []
                        s = []
                        y = []
                        gamma = 1
                else:
                        skMinusm = s[len(s)-1]
                        ykMinusm = y[len(y)-1]

                        # Skip the update if negative curvature
                        # was detected in the previous step.
                        #if( skMinusm.dotProduct(ykMinusm) > 0):

                        gammaTop = skMinusm.dotProduct(ykMinusm)
                        gammaBottom = ykMinusm.dotProduct(ykMinusm)
                        gamma = gammaTop/gammaBottom

                ####################################################
                #       Get the descent direction
                ####################################################
                universeCopy = copy(universe)
                energy, gradient = universeCopy.energyAndGradients()
                q = gradient
                alpha = []

                # Start from len(s)-1 in steps of -1 until 0
                for i in range(len(s)-1, -1, -1):
                        alpha.append(rho[i]*(s[i].dotProduct(q)))
                        #if( skMinusm.dotProduct(ykMinusm) > 0):

                        gammaTop = skMinusm.dotProduct(ykMinusm)
                        gammaBottom = ykMinusm.dotProduct(ykMinusm)
                        gamma = gammaTop/gammaBottom

                ####################################################
                #       Get the descent direction
                ####################################################
                universeCopy = copy(universe)
                energy, gradient = universeCopy.energyAndGradients()
                q = gradient
                alpha = []

                # Start from len(s)-1 in steps of -1 until 0
                for i in range(len(s)-1, -1, -1):
                        alpha.append(rho[i]*(s[i].dotProduct(q)))
                        q = q - alpha[len(s)-1-i]*y[i];

                alpha.reverse()
                r = gamma*q

                # Start from 0 in steps of 1 until len(s)-1
                for i in range(0, len(s)):
                        beta = rho[i]*(y[i].dotProduct(r))
                        r = r + s[i]*(alpha[i] - beta)

                ####################################################
                # Get the step size
                ####################################################
                initialAlfa = 1
                alfa = initialAlfa
                direction = -r
                alfa = backtrackingLineSearch( universeCopy, initialAlfa,
direction, maxIter)
                ####################################################
                # Update the universe
                ####################################################
                universe.addToConfiguration(alfa*direction)

                newEnergy, newGradient = universe.energyAndGradients()
                newConfiguration = universe.copyConfiguration()

                ####################################################
                # Update/discard the stored vectors
                ####################################################
                # Discard the first vector pair from storage
                if( k >= m):
                        del s[0]
                        del y[0]
                        del rho[0]

                s.append(newConfiguration - configuration)
                y.append(newGradient - gradient)

                totalStoredVectors = len(y)
                val =
y[totalStoredVectors-1].dotProduct(s[totalStoredVectors-1])
                rho.append(1/val)

                #######################################################
                # Update the initial gradient, energy and configuration
                # to the current one!
                #######################################################
                gradient = newGradient
                energy = newEnergy
                configuration = newConfiguration

                k=k+1

def backtrackingLineSearch( universeCopy, initialAlfa, direction, maxIter):
        # Initialization
        rho = 0.05
        c = 0.9
        maxit = 100

        energy, gradient = universeCopy.energyAndGradients()
        Dphi0 = gradient.dotProduct(direction)
        phi0 = energy
        alfa = initialAlfa

        if ((alfa <= 0) or (Dphi0 > 0)):
                print "Initialization of step incorrect!"

        for i in range(maxit):
                universeCopy.addToConfiguration(alfa*direction)
                new_energy, new_gradient = universeCopy.energyAndGradients()

                #print "New Energy=", new_energy
                #print "New gradient=", new_gradient.norm()

                if(new_energy >= phi0 + alfa*c*Dphi0):
                        universeCopy.addToConfiguration(-alfa*direction)
                        alfa = rho*alfa
                else:
                        break

        # Reset the value of alfa back to its original value
        # in case the max iterations have been reached and if
        # no decrease in function value could be obtained with alfa.
        if( (new_energy >= phi0 + alfa*c*Dphi0) & (i >= maxit) ):
                alfa = initialAlfa

        #print "alfa =", alfa
        return alfa





More information about the mmtk mailing list