8
$\begingroup$

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[...]: uses NelderMead to determine optimal parameter set for problem stated in question above
  • hamQ1,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...

|improve this question
$\endgroup$
  • $\begingroup$ Providing your code will definitely help others reproduce / solve this problem. $\endgroup$ – dr.blochwave Feb 2 '15 at 15:23
  • $\begingroup$ @blochwave I will try to reduce the code to the least possible amount and attach it within the next two hours. $\endgroup$ – Lukas Feb 2 '15 at 15:27
  • 1
    $\begingroup$ I don't know whether this comprises a natural consequence, but the equation preprocessing, and in particular event detection, ties itself in knots trying to make sense of (non)events due to (I think) imaginary parts becoming negative. Why this happens is outside of my understanding so I pointed it out to someone more familiar with those workings. Also there was some oddity with packed arrays that actually caused the crash and I think that goes away when one uses explicitly real input and variables. And there is apparently yet another bad spot which remains under investigation. $\endgroup$ – Daniel Lichtblau Feb 3 '15 at 16:24
  • 1
    $\begingroup$ Here are a couple more observations. (1) The eventual crash I've been seeing might be self inflicted (by you, not me). Setting the stiffness detector to False is problematic because at optimization step 98, 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 full NDSolve interpolation, just the end result. So dense output is not needed. Can avoid it using NDSolveValue. $\endgroup$ – Daniel Lichtblau Feb 3 '15 at 17:19
  • 1
    $\begingroup$ The "other bad spot" was a reference to what I later said, that with some fixes in place we still have trouble at step 98. Was due (I think) to the very small steps causing need for considerable memory to store the DE result. $\endgroup$ – Daniel Lichtblau Feb 3 '15 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.