[MMTK] Problems adding configuration to Universe Copy

vanitha@cs.wisc.edu vanitha at cs.wisc.edu
Wed Feb 8 19:33:10 CET 2006


Sorry for the copy and paste error in the code snippets. Some of the
sections of the code are repeated twice. I wouldn't want to resend this
again to the mailing list but the order is as follows:

1) Termination criteria
2) Choose initial approximation to inverse Hessian
3) Get the descent direction
4) Get the Step size
5) Update the universe
6) Update/Discard stored vectors
7) Update initial energy, gradient and configuration to current one.

Hope this is clear!

Thanks,

- Vanitha


> 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
>
>
>
> _______________________________________________
> mmtk maillist  -  mmtk at starship.python.net
> http://starship.python.net/mailman/listinfo/mmtk
>





More information about the mmtk mailing list