The problem

I have been wanting to make my code which solves an optimisation problem faster. This led to this question on here and I realised from there that the root of the problem is the very slow FindMinimize. Then, using Compile and Parallelization->True together seemed like the best and one of the last options remaining.

Now at this point, I must say that I am a total noob in this area. I don't know how (maybe because I had already executed the old code without Compile before and variables were stored in memory), but I managed to get my compiled function to run in parallel (with some help of Evaluate to avoid errors which came up otherwise) and execute without any errors and the runtime was almost 100,000 times faster. But the joy was shortlasting because restarting the kernel and running the code again threw many errors and I suspect that using Evaluate might have been the problem. If this is all getting too abstract, below is my exact code (without Evaluate since I can't make it help right now anyway) and I have mentioned errors I got as well.

So, how do I get Compile to work here?

My code

The first part of the code is simply data given so that you can run the code. You may as well ignore the entire Part 1 of code since it's just constants and variables to make the code complete. If you hopefully manage to run a modified code without errors, you may use ic={0,0} as an input to the function.

(*PART I*)

inputs1 = Table[Unique["u"], {2}];
ga={{0.0071, -0.0564, 0.00580, 0.0104, 0.00240, -0.0321}};
zminus2M = ConstantArray[0, {4, 1}];
P={{6.78614, -0.962986, 4.91036, -0.839262, 
  1.84927, -0.218886}, {-0.962986, 1.8242, -1.07876, 
  2.54396, -0.536404, 1.83934}, {4.91036, -1.07876, 3.70497, -1.16611,
   1.4347, -0.571865}, {-0.839262, 2.54396, -1.16611, 
  3.64699, -0.633915, 2.69184}, {1.84927, -0.536404, 
  1.4347, -0.633915, 0.56924, -0.362581}, {-0.218886, 
  1.83934, -0.571865, 2.69184, -0.362581, 2.02307}};
windowfunction = Array[HammingWindow, 2*M + 1, {-1/2, 1/2}];
X={{3.03185*10^-6, 1.29453*10^-6, 2.80629*10^-7}, {1.29453*10^-6, 
  7.31323*10^-7, 1.01069*10^-7}, {2.80629*10^-7, 1.01069*10^-7, 
Q = {{0, 0}, {0, 0}};
R = {{1}};
Qf={{3231.8, 4.28555}, {4.28555, 3231.8}};
statesxu = Sxu.inputs1;
outputsuz = Suz.inputs1;


optimalinputs1 = 
  Compile[{{ic, _Real}},
    Block[{zetas1 = {}, zetas, zetap1, zetap, states1, states, 
      outputs, outputs0, obj, speccons, setcons, doms, constraints, 
      shiftedstates, shiftedinputs},
     states1 = Sxp.ic + statesxu;
     states = 
      Table[states1[[i ;; i + n - 1]], {i, 1, (hor + 2*M)*n + 1, n}];
     outputs0 = Sxz.states1 + outputsuz;
     outputs3 = 
      Table[outputs0[[i ;; i + o - 1]], {i, 1, (hor + 2*M)*o+1,o}];
     outputs = 
      ArrayFlatten[{{zminus2M}, {outputs3}}];
     zetas1 = 
          windowfunction[[1 ;; 2*M]]*outputs[[p + M + 1 ;; p+3*M]], 
          1], {p, -M, hor + M}]];
     zetas = 
      Table[zetas1[[i ;; i + x - 1]], {i, 
        1, (2*M + 1)*(hor + 2*M + 1)*x - x + 1, x}];
     shiftedstates = Table[states[[i]] - xs, {i, hor}];
     shiftedinputs = inputs1 - us; 
    obj = Sum[shiftedstates[[i]].Q.shiftedstates[[i]], {i, hor}] + 
   Sum[(shiftedinputs[[i]]*R*shiftedinputs[[i]])[[1]][[1]], {i, 
     hor}] + (states[[hor + 1]] - xs).Qf.(states[[hor + 1]] - 
      xs) + (xs - xshat).Ta.(xs - xshat);
     speccons = 
      Block[{itr, outputsp, zetaoutputsp, zetaoutputsp1, cons, 
        speccon = True},
         itr = p + M + 1;
        zetap = zetas[[(itr - 1)*2*M + itr ;; itr*2*M + itr]];
        outputsp = outputs[[p + M + 1 ;; p + 3*M + 1]];
        zetaoutputsp = Transpose[{zetap, outputsp}];
        zetaoutputsp1 = 
         Table[ArrayFlatten[zetaoutputsp[[i]], 1], {i, 2*M + 1}];
        cons = Sum[
              windowfunction[[k]]*de}}]).zetaoutputsp1[[k]], {k, 
         2*M + 1}] + Last[zetap].P.Last[zetap] <= alpha;
   speccon = speccon && cons, {p, -M, hor + M}]; speccon];
     setcons = Block[{con = True},
    con=con&&And @@ Table[-15 <= states[[i]][[1]] <= 15 && -100 <= 
    states[[i]][[2]] <= 100, {i, 2, hor}]; 
    con=con&&And @@ Table[-10 <= inputs1[[j]] <= 10, {j, hor}];; 
     constraints = 
      Simplify[setcons && 
       speccons && (Join[states[[hor + 1]] - xs, {theta}]).(Inverse@X).(Join[states[[hor + 1]] - xs, {theta}]) <= 1];
     optvec = 
      inputs1 /. 
       Last[FindMinimum[{obj, constraints}, inputs1, 
         MaxIterations -> 10000, Method -> "InteriorPoint"]];
     optvec], RuntimeAttributes -> {Listable}, 
   Parallelization -> True];

The errors on running above code

Compile::part: Part specification outputs[[p+2+1;;p+3 2]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.
Compile::part: Part specification outputs[[p+2+1;;p+3 2]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.
Compile::part: Part specification outputs[[p+M+1;;p+3 M+1]] cannot be compiled since the argument is not a tensor of sufficient rank. Evaluation will use the uncompiled function.
General::stop: Further output of Compile::part will be suppressed during this calculation.
Compile::cset: Variable region of type True|False encountered in assignment of type _Real.
Compile::cif: The types of the two results in If[con,{-15,-100}\[VectorLessEqual]states[[i]]\[VectorLessEqual]{15,100},False] are incompatible because their ranks are different. Evaluation will use the uncompiled function.
Compile::cset: Variable con of type True|False encountered in assignment of type _Real.

I am grateful for your time and efforts to go through this longish piece of code.


As suggested by Henrik in this answer, in the updated code above, I got rid of domain constraints (the inputs are Real constraints), replaced Do loops in setcons with And@@ and Table and finally used Simplify before supplying constraints to FindMinimum and the results are great (about factor of 10 improvement for horizon = 10). However, my actual problem has horizon = 50 and the new code still takes an eternity to solve it, so that problem is to be solved still.

  • 3
    $\begingroup$ Compile is not a silver bullet that will magically make your code faster. In particular only certain functions can be compiled. FindMinimum is not included among those so its inclusion in a compiled function will force it to drop back to uncompiled evaluation, and ultimately your efforts will be for nothing if that is the bottleneck. $\endgroup$
    – MarcoB
    Commented Mar 31, 2021 at 12:37
  • 1
    $\begingroup$ @MarcoB I didn't know about this until now. I got excited because somehow the first time I tried this, my optimisation finished in only 0.000xx seconds instead of 12 seconds for the same parameter values but I don't know why and it never worked again after restarting the kernel. Are you aware of any other method to parallelize the evaluation of FindMinimum? $\endgroup$
    – ModCon
    Commented Mar 31, 2021 at 12:45
  • $\begingroup$ Parallelization is unlikely to help since most iterative optimization algorithms lend themselves poorly to it (the next step depends on the results of the previous one, so they cannot be run independently in parallel). On this site you can find a compiled version of the Nelder-Mead optimization algorithm; at the moment that implementation does not accept constraints, but you could include them "by hand" by generating an appropriate objective function satisfying the KKT conditions yourself. $\endgroup$
    – MarcoB
    Commented Mar 31, 2021 at 13:13
  • 2
    $\begingroup$ I think the actual problem lies somewhere else. See mathematica.stackexchange.com/a/243819. $\endgroup$ Commented Mar 31, 2021 at 14:09
  • $\begingroup$ @MarcoB I see the rather grim picture now. I will keep trying different things and hopefully some of those give good results. $\endgroup$
    – ModCon
    Commented Mar 31, 2021 at 16:22