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:
Then we'll consider the bound-state
Schrodinger equation in momentum
space, which takes the form of an integral equation.
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
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
- 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).
- 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?
- 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?
- 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.
- Bonus: Switch the code to use classes defined in previous
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.]
- What are the units of lambda in the delta-shell potential?
[Note: hbar=1 here.]
- 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?
- 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?
- 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?)
- 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?
- 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
How could you improve the numerical accuracy of the calculation?
(For example, would it help to change "job" when calling "gauss"? [Hint,
- What does the momentum-space wave-function look like (it is
printed out to a file)?
- Does it fall off at large k? Does it oscillate? Is it well behaved
at the origin? Can you account for its characteristics by some simple
- The analytic wave function in coordinate space is given in
the background notes. Try to calculate the momentum
space wave function from this (as discussed in the notes)
- Is the momentum wave function normalized? What is the
- Bonus: Verify that, regardless of the potential's strength, there
is only a single bound state.
- Bonus: Modify the code to calculate the l=1 bound state.
- 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.