I am making use of the Nelder Mead algorithm inside Mathematica's NMaximize
to find a parameter set by solving the Schrodinger equation for a given system. (see below) With a 3D system and a set of 5-9 parameters everything works well but for a 9D system the kernel crashes during the optimisation without throwing an error.
Here is what I know so far:
- I could test it on Mathematica 8.0.1 (only Linux) and 10.0.1 (Linux and Windows). Each resulting in kernel crash at the same number of iterations
- Looking at the memory usage of the MathKernel in Linux's list of processes showed a memory consumption of ~250MiB which (to me) does not appear to be extremely high.
- The time evolution of my 9D compound system takes ~1s, so it should not be an issue I assume. At first I thought that solving the differential equations might cause some trouble...
I am completely clueless how to approach that problem, so I'd appreciate any help. If you need code or some output of "debugging commands", just let me know. The code is a bit longish, so I wanted to stay as general as possible in the beginning.
Some explanation about physical background
I have two Hamiltonians $\hat{\rm H}_1$ and $\hat{\rm H}_2$ (from a 3D Hilbert space), both depending on the same parameters $a_{1\ldots n}$. If I use NelderMead
to find optimal values for $a_{1\ldots n}$ in the sense that either $\hat{\rm H}_1$ or $\hat{\rm H}_2$ causes a given time evolution everything works fine. But if I try to do the same for the compound system $\hat{\rm H}=\hat{\rm H}_1 \otimes \hat{\rm H}_2$ the kernel crashes after a few iterations. That means I want to find the parameters $a_{1\ldots n}$ in $\hat{\rm H}(t)$ such that:
$$i\frac{\mathrm{d}}{\mathrm{d}t}\hat{\rm U}(t)=\hat{\rm H}(t)\hat{\rm U}(t),\quad t \in [0,T] \quad (1)$$
yields a certain evolution $\hat{\rm U}(T) = \hat{\rm U}_{\rm target}$.
Edit 1 - Link to .nb file and brief description
In order to reproduce the problem, please find the relevant code here. I tried to further reduce it but less complexity seemed to resolve the issue with kernel crashing. So I reduced it so far that 3D problem works well and 9D crashes after 20 iteration steps. This is a link to virustotal analysis of the Dropbox link above in case someone is suspicious: http://preview.tinyurl.com/pz9jmcd
What the code does Mainly every function is briefly documented inside the nb file, but these are the most important ones
solver[H_,time_,index_]
: computes time evolution according to Eq.(1) as replacement rules for matrix elements $u_{i,j}(t)$ of $\hat{\rm U}(t)$optimizer[...]
: usesNelderMead
to determine optimal parameter set for problem stated in question abovehamQ1,hamQ2,ham
: Hamiltonians representing $\hat{\rm H}_{1,2}$ and $\hat{\rm H}$ of the question. They are composed of functions depending on certain parameters
Edit 2 - Link to new .nb file with ASCII only and no global variables
As recommended by DanielLichtblau I rewrote the code with ASCII symbols and removed the global variables. However the behaviour is exactly the same. Please find the new notebook file here (The old one from Edit1 is still accessible).
Edit 3 - Possible fix by using Daniel Lichtblau's suggestion
Daniel Lichtblau mentioned here that a fix might be to remove Re[]
, Im[]
and so forth from the input given to my ODE solver. As a test I changed the following two lines in the .nb file linked in "Edit2 ":
ex1[t_, tg_, exp_] := Ox1[t, tg, exp] - b2*D[Ox1[x, tg, exp], {x, 2}] /. {x -> t};
ey1[t_, tg_, exp_] := -b1*D[Ox1[x, tg, exp], x] /. {x -> t};
Indeed, that seems to solve the problem. Also replacing the SetDelayed
in the original definitions (containing Re,Im
) of ex1,ey1
with Set
seems to solve the issue as no Re,Im
are passed to the ODE solver.
The question is: Why does Re,Im
etc cause a crash somewhere during computation? If it did at the beginning I would still not understand but it would be more reasonable to me than this...
NDSolve
gets into very small steps (order of 10^(-6)) and thus is creating a very dense result. When I remove that setting, step 98 is slow but runs to completion along with the computation. (2) You are not using the fullNDSolve
interpolation, just the end result. So dense output is not needed. Can avoid it usingNDSolveValue
. $\endgroup$ – Daniel Lichtblau Feb 3 '15 at 17:19