6810: 1094 Session 17

Handouts: excerpt on simulated annealing

In this session, we'll first consider a toy problem comparing multidimensional minimization with standard downhill routines for local minima to an application of Monte Carlo techniques to find the global minimum: simulated annealing. Then we'll consider the bound-state Schrodinger equation in momentum space, which takes the form of an integral equation.

Simulated Annealing

Standard optimization methods are very good in finding local minima near where the minimization was started, but not good at finding the global minimum. Finding the global minimum of a function (such as the energy) is often (but not always) the goal. One strategy using conventional minimizers is to start the minimization at many different places in the parameter space (perhaps chosen at random) and keep the best minimum found.

Here we'll adapt our Monte Carlo algorithm for generating a canonical Boltzmann distribution of configurations at a temperature T to mimic how physical systems find their ground states (i.e., the energy minimum at T=0). At high temperature (compared to characteristic energy spacings), the equilibrium distribution will include many states. If the system is cooled slowly, then it will have the opportunity to settle into the lowest energy state as T goes to zero. This is called annealing. If the system is cooled quickly, it can get stuck in a state that is not the minimum ("quenching"); this is analogous to the routines we looked at in Session 10, which rapidly go "downhill" to find local minima.

The strategy here is to simulate the annealing process by treating the function to be minimized as an energy (it might actually be an energy!), introducing an artificial temperature T, generating a sequence of states in a canonical distribution via the Metropolis algorithm. Then we lower the temperature according to a "schedule" and let the system settle into (hopefully!) a configuration that minimizes the energy.

Global vs. Local Minimum

We first consider the simple example of a damped sine wave in one dimension, which has many local minima.
f(x) = e-(x-1)2sin(8x)

(This example is from the GSL documentation.) Suppose we start with an initial guess of 15.5; let's see how "downhill" minimization methods compare to simulated annealing.

  1. Plot the function in Mathematica (see test_function.nb) so we know what to expect ahead of time (remember that in general we won't be able to do this!). Use FindMinimum starting near 15.5 to find the local minimum. Then find the global minimum (using the graph to identify a good starting value). Answers:

  2. The code multimin_sa_compare.cpp (compile and link with make_multimin_sa_compare) applies the GSL routines we used in a previous session to find the minimum and then applies a simple simulated annealing algorithm. Run the program and look at the code to figure out where the simulated annealing control parameters (starting and finishing temperatures, maximum step size, cooling rate, # of iterations at each temperature) are set. How well does each method work?

  3. Now adjust the simulated annealing control parameters to give it a chance! The step size should be large enough (or small enough) so that successes are neither assured nor happen too infrequently (50% would be ideal!). The initial temperature should be considerably larger than the largest typical delta_E normally encountered. You shouldn't cool too rapidly. Make these adjustments and rerun the code. What did you do to get it to work best?

  4. Just to check that the Metropolis algorithm is doing something, once the code is working, change the line with "random < exp()" to "random > exp()" and verify that it no longer works. What happened?

  5. Bonus: Switch the code to use classes defined in previous sessions.

Bound States in Momentum Space: Delta-Shell Potential

As discussed in the background notes, the bound-state Schrodinger equation in momentum space can be solved as a matrix eigenvalue problem. We have three versions of a code to solve it: two with the matrix in symmetric form, which use GSL routines and LAPACK routines, respectively, and one that uses a LAPACK routine for non-symmetric matrices. We also dynamically allocate arrays using "new" and "delete".

The delta-shell potential is a convenient choice for a test potential because we can easily find the bound-state eigenvalue (numerically) and it takes a simple form in momentum space. It has only one bound state; if we only print out the negative-energy (E<0) eigenvalues then we should find just one. See the background notes for additional details. [Warning: It is not a well-behaved potential in momentum space, so you will be challenged to find good solutions.]

  1. What are the units of lambda in the delta-shell potential? [Note: hbar=1 here.]

  2. For mu=0.5, b=2.0, and lambda=-1.0, use Mathematica to find the correct eigenvalue for the (single) bound state in the delta-shell potential by using the appropriate equation in the background notes. [Hint: FindRoot] What is the answer?

  3. Look at deltashell_boundstates_gsl.cpp. Note how the non-gsl arrays are dynamically allocated and removed using "new" and "delete". If we had an int called size and wanted to create an array called my_array with dimension size, what command would we give?

  4. Compile and link deltashell_boundstates_gsl.cpp using the corresponding makefile. Run with lambda=-1.0 but specify on your own the number of gauss points to use and the upper bound for the integrations of momentum. The latter may be much larger than you would guess because the delta-shell potential is actually rather poorly behaved in momentum space! Based on how the potential falls off with k, what do you think is a reasonable momentum cutoff? (What should you do to discover how the potential falls off?)

  5. Now look at deltashell_boundstates_lapack.cpp. This is essentially the same code but it uses a Fortran LAPACK function instead of a GSL function to diagonalize the matrix. (Note: You have to be on a Linux machine for this to work.) We'll also use the Intel compiler and the MKL library. Load them with:
    module load intel-12.1.3-64
    Look at the corresponding makefile to see how the code is compiled and linked. Can you reproduce the results from the GSL code?

  6. Run one of the codes with different numbers of gauss points and different momentum cutoffs to try to improve the predicted eigenvalue. Estimate the precision of your best value for the bound-state energy. How could you improve the numerical accuracy of the calculation? (For example, would it help to change "job" when calling "gauss"? [Hint, hint.])

  7. What does the momentum-space wave-function look like (it is printed out to a file)?
  8. Bonus: Verify that, regardless of the potential's strength, there is only a single bound state.
  9. Bonus: Modify the code to calculate the l=1 bound state.
  10. Bonus: Make a version of the code with classes that can use either GSL or LAPACK.

6810: 1094 Session 17. Last modified: 09:53 am, April 19, 2013.