# How can I speed up this NM optimization? Most delay seems to be due to long analytic expression

I am struggling with computational speed of a NelderMead optimization. Nevertheless my impression is that I lose lots of speed not only due to the NM optimization itself but rather to for instance inefficient memory usage. Let me briefly explain what I intend to do and what the notebook looks like (the notebook file can be downloaded from here. Alternatively you can use the MWE in "Edit 3" section below):

• I have two 3x3 matrices (Hamiltonians) hamQ1,hamQ2 that depend on some parameters $\Lambda_2$ and $c_{ij}$ (see sections 1.2.2 & 1.2.3 in my .nb)
• By solving the Schroedinger equation in a time interval $[0,t_g]$ (let us focus on $t_g=20$) and evaluating a fidelity function, I want to optimize for the parameters $\Lambda_2$ and $c_{ij}$. The method that solves Schroedinger equation is to be found in section 2.2, the fidelity function in 2.3 (it mainly computes the "overlap" between two unitaries)
• The optimization function is located in section 3.1. A call of that function together with constraints is in section 3.2
• The notebook contains initialization cells, which need roughly 60s on my machine since they involve some longish analytical integrands (sec. 1.2.2)

All of the above works but is not as fast as I want it to be. From the second cell in sec. 3.2 one can see that only evaluating one of my Hamiltonians takes ~4.5s if numerical values for $c_{ij}$ and $\Lambda_2$ are given. If they are still symbolic, hamQ1[t,tg] needs ~7.5s to be evaluated. Solving the Schroedinger equation is quite efficient, I would say. It only takes ~0.1s. If you now execute the first cell of sec. 3.2 you will notice that all works but each step of the NelderMead optimization takes ~1.5min on my machine.

# Questions

• Is it possible to somehow monitor where most time is spent during each NealderMead step? I need to do that optimization for many times $t_g$ and maybe for more parameters $c_{ij}$. In the notebook file there are only 5 parameters I am optimizing for, but that will probably increase to ~10 which will significantly increase computational time
• The programme is quite memory expansive since each Hamiltonian has a ByteCount of: ByteCount[hamQ1[t, tg]] (* 2 166 116 960 *). Maybe it could be sped up by optimizing memory consumption (?)

# Edit 1/Comment

From MichaelE2's comment I guess it was not clear enough that I linked the notebook file. This is simply because I did not want to post all the code lines here. The notebook is structured with sections and the most relevant passages should be easy to find/understand with the help of the "summary" above. Again, please find the notebook here: https://www.dropbox.com/s/3ei8r69mf1mb2up/Fidelity.nb?dl=0

# Edit 2:

As suggested by DanielLichtblau I rewrote the code in a way that NIntegrate instead of Integrate is used which in the end avoids memory consuming (long) symbolic expressions. Indeed, that gave a tremendous speed up. Nevertheless I will try to work out a MWE containing symbolic expressions since I am still interested if there is a solution that still involves symbolic computations.

# Edit 3 (minimal working example)

As requested here is a MWE that still shows similar behavior of the whole notebook I linked to above. This code is still slow but by far not as slow as the original one. As stated in my second edit I was abloe to work around the issue by following DanielLichtblau's remark about using only numerical integration. But still, I am interested in a way to speed up the analytic approach.

(* constant parameters *)
capDel = 2*\[Pi]*(-35/100);
del = 2*\[Pi]*(45/1000);
capLam1 = 0;

(* g depends on capLam2, which will be one of the optimized parameters *)
g = capDel - del + capLam1 - capLam2;

(* Define functions inside Hamiltonian (epsX,epsY) *)

hann[t_, tg_] := hann[t, tg] = (c11*(1 - Cos[(2*\[Pi]*t)/tg]) +  c13*(1 - Cos[(2*3*\[Pi]*t)/tg]));
omega[t_, tg_] := omega[t, tg] = hann[t, tg] - I*D[hann[x, tg], {x, 2}]/. x -> t;

A2g[tg_] := A2g[tg] = Integrate[Expand@TrigToExp@(Exp[I*g*t]*omega[t,tg]), {t, 0, tg}];
a2[tg_] := a2[tg] = \[Pi]*(1 -A2g[tg]/(capLam2 - del)*Integrate[Expand@TrigToExp@(Exp[I*(del - capDel - 2*capLam1 + capLam2)*t]*omega[t, tg]), {t, 0, tg}])^-1;
epsX[t_, tg_] := epsX[t, tg] = ComplexExpand[Re[a2[tg]*omega[t, tg]]];
epsY[t_, tg_] := epsY[t, tg] = ComplexExpand[Im[a2[tg]*omega[t, tg]]];

(* demonstrational Hamiltonian *)

hamilton[t_, tg_] := hamilton[t, tg] = {{epsX[t, tg], epsX[t, tg] - I*epsY[t, tg]}, {epsX[t, tg] + I*epsY[t, tg], 0}};

(* solve Schrodinger equation*)
(* USAGE: {foo,{solutions}} = solverSample[...], where {solutions} is structured as follows: {time poisition,{time evolution op. at particular time}} *)

solver[H_, time_, sampleStep_] :=
Module[{d, init, eqs, vars, solargs, t = time[], t0 = time[], lastt = -sampleStep, dt = sampleStep},
d = Dimensions[H][];
u[t_] := Table[Subscript[u, i, j][t], {i, 1, d}, {j, 1, d}];
init = Thread[Flatten /@ (u[t0] == IdentityMatrix[d])];
eqs = Thread[Flatten /@ (I*u'[t] == H.u[t])];
vars = Flatten[Table[Subscript[u, i, j], {i, 1, d}, {j, 1, d}]];
solargs = Join[eqs, init];
Return[
Reap@NDSolve[solargs, {}, time,
InterpolationOrder -> All,
AccuracyGoal -> 12,
PrecisionGoal -> 12,
MaxSteps -> \[Infinity],
EvaluationMonitor :> If[t >= lastt + dt, lastt = t; Sow[{t, Table[ Subscript[u, i, j][t], {i, 1, d}, {j, 1, d}]}]]]
]
];

(* demonstrational fidelity function *)
fidel[evol_, uIdeal_] := fidel[evol, uIdeal] = Abs[Tr[uIdeal.evol]];

(* function to optimize for certain parameters *)
optimizer[gateTime_, ham_, ideal_, vars_, steps_, cons___] := Block[{obj, nbr = 0, tstart},
obj[v_ /; VectorQ[v, NumericQ]] := Module[{foo, sol, final},
{foo, {sol}} = solver[ham /. Thread[vars -> v] /. tg -> gateTime, {t, 0, gateTime}, 10^(-3)];
final = Last@Last@sol;
Return[fidel[final, ideal]];];
solsOpt = NMinimize[
{obj[vars], cons},
vars,
Method -> "NelderMead",
AccuracyGoal -> 5,
EvaluationMonitor :> {nbr += 1; If[Mod[nbr, steps] == 0, Print["Step: ", nbr, " ; Current magnitude: ", obj[vars]], Print]}];
Print["Final parameter set: ", solsOpt[]];
Return[solsOpt];]

(* perform optimization *)
constraint = -2 <= capLam2 <= 8 && -4 <= c11 <= 4 && -4 <= c13 <= 4;
optimizer[20, hamilton[t, tg], IdentityMatrix, {capLam2, c11, c13}, 1, constraint];

• If people are voting to close because of a lack of code, I would point out the code is in the linked file. – Michael E2 Apr 9 '15 at 14:33
• If people are voting to close because of a lack of code, I would point out the code is in the linked file. – Michael E2 Apr 9 '15 at 14:33
• Maybe define the various functions so that they only operate if given explicitly numeric inputs, and use NIntegrate instead of Integrate (if you are recomputing integrals, I wasn't sure). Also you might want to memoize objective function values e.g. obj[vars_] := obj[vars] = ... This will avoid reevaluations (though at a memory cost, so you may sometimes need to flush them, but that's another issue). – Daniel Lichtblau Apr 9 '15 at 16:31
• Is there any way you could post a minimal working example as code in your question? Otherwise this post will be of very limited use for other users browsing the site, especially if the linked file at some point gets deleted. – Jens Apr 9 '15 at 17:49
• @DanielLichtblau The main idea about integrals is that all those factors ´a´ and ´A correspond to amplitudes needed to satisfy some area condition of my ´[Epsilon] pulses. Initially I thought it would be best to pass a full analytical expression to NMinimize in order to achieve best results - however it is very slow. Thanks for your tip about memoizing the function values. Concerning your first tip about redefining functions: would that involve that I need to e.g. turn my hamQ1[t,tg] to hamQ1[t,tg,c11_?NumericQ,...]? Just to ensure I got your point correctly. – Lukas Apr 9 '15 at 18:28