# Minimization by Nelder-Mead

Finding a global minimum for this problem (non-linear optimization by the Nelder-Mead downhill simplex method) may not be possible, but by finding local minimum, I am expecting the value of the function at the minimum is around 1 or (far) less than 1 (the lower the value, the better is the answer).

"Nelder-Mead" method is available for NMinimize in Mathematica. The problem I am facing is, even though I am changing the method options for Nelder-Mead (for example, "ShrinkRatio", "ContractRatio", "ReflectRatio", and so on) and also initializing the calculation for different initial values (which is chosen by the "RandomSeed"), I am getting minima which are way greater than 1.

I found people deal these kind of problems by writing codes either in C++ or Fortran and the results are good, but I know neither of them so I am trying to do it with Mathematica using the built-in functions.

Here is the what I am trying (a simplified version of the original problem that I have to deal with which has even more number of parameters):

A = ({{a11, 0},{0, a22}});
B = ({{b11, b13},{b13, b22}});
M1 = A + Exp[I phi] B;
M2 = a A + b Exp[I theta] B;
M3 = a A - 3 b Exp[I theta] B;
m1 = Abs[Eigenvalues[M1]];
m2 = Abs[Eigenvalues[M2]];
m3 = Abs[Eigenvalues[M3]];

t1 = T1[a11_?NumberQ, a22_?NumberQ, b11_?NumberQ, b22_?NumberQ,
b13_?NumberQ, phi_?NumberQ] = m1[[1]];
t2 = T2[a11_?NumberQ, a22_?NumberQ, b11_?NumberQ, b22_?NumberQ,
b13_?NumberQ, phi_?NumberQ] = m1[[2]];

t3 = T3[a11_?NumberQ, a22_?NumberQ, b11_?NumberQ, b22_?NumberQ,
b13_?NumberQ, theta_?NumberQ, a_?NumberQ, b_?NumberQ] = m2[[1]];
t4 = T4[a11_?NumberQ, a22_?NumberQ, b11_?NumberQ, b22_?NumberQ,
b13_?NumberQ, theta_?NumberQ, a_?NumberQ, b_?NumberQ] = m2[[2]];

t5 = T5[a11_?NumberQ, a22_?NumberQ, b11_?NumberQ, b22_?NumberQ,
b13_?NumberQ, theta_?NumberQ, a_?NumberQ, b_?NumberQ] = m3[[1]];
t6 = T6[a11_?NumberQ, a22_?NumberQ, b11_?NumberQ, b22_?NumberQ,
b13_?NumberQ, theta_?NumberQ, a_?NumberQ, b_?NumberQ] = m3[[2]];

O2 = 0.235; O1 = 72;
E2 = 0.031; E1 = 0.87;

O4 = 0.0222196; O3 = 0.958157;
E4 = 0.006198; E3 = 0.018862;

O6 = 0.0941698; O5 = 1.60087;
E6 = 3.6079*^-7; E5 = 0.00026;

function[a11, a22, b11, b22, b13, a, b, phi,theta] =
((t1 - O1)/E1)^2 + ((t2 - O2)/E2)^2 + ((t3 - O3)/E3)^2 + ((t4 - O4)/E4)^2 +
((t5 - O5)/E5)^2 + ((t6 - O6)/E6)^2;

Do[Print[NMinimize[
function[a11, a22, b11, b22, b13, a, b, phi, theta],
{a11, a22,b11, b22, b13, a, b, phi, theta},
Method -> {"NelderMead", "ShrinkRatio" -> 0.5,"ContractRatio" -> 0.7,
"ReflectRatio" -> 1, "ExpandRatio" -> 2,"RandomSeed" -> j}]], {j, 5}]


Five random iterations (with different initial conditions chosen randomly) find different minima but all are much greater than the expected result. Here they are:

(* {9081.63
{a11->0.0585825,a22->1.22145,b11->-0.26137,b22->0.327416,b13->-0.158818,
a->-1.2025,b->0.168387,phi->-1.01436,theta->0.698082}}

{32448.8
{a11->-0.542868,a22->0.845826,b11->1.10494,b22->0.739838,b13->-0.16513,
a->0.795721,b->0.346771,phi->-0.369352,theta->-0.077746}}

{7122.85
{a11->0.0846693,a22->-0.643122,b11->-0.0834931,b22->-1.14742,b13->-0.37268,
a->1.36398,b->-0.208855,phi->-0.367912,theta->-0.359771}}

{14919
{a11->-0.14745,a22->-0.0507363,b11->0.414601,b22->0.306933,b13->0.420321,
a->0.364225,b->0.663577,phi->-0.625049,theta->-0.33245}}

{13582.6
{a11->-0.182997,a22->0.431253,b11->0.201061,b22->0.507661,b13->0.258584,
a->-0.471931,b->0.741867,phi->-0.5953,theta->-0.0911202}} *)


I tried with different combinations of the Nelder-Mead options with no improvements. The default value of the MaxIterations is 100, and I also tried by increasing the number of iterations by using the package of Nelder-Mead from the excellent post Shaving the last 50 ms off NMinimize by Oleksandr R., but again did not get any improvement. Can anyone help me with what I am missing?

• Yes, I have seen that and at first I forgot to mention that, so I just edited my question few mins ago, where I included the fact that, I also used that package and by changing the number of Max Iteration, I did not get any improved result. May be still I will be able to get good result by using that package, but could not figure out what I am missing. @rm-rf – SAS Jun 28 '13 at 20:12
• @MichaelE2, Oh, Really ?!! I am using version 8. So you are saying version 9 is giving god results ! Btw, are you saying you ran my code without any modification ? If so then I have to try to get Mathematica 9. Thanks for your useful reply. – SAS Jun 28 '13 at 21:40
• @MichaelE2, thanks! though I am using version "8", but it is the "student" version, I am guessing it could be another reason. Let me get some other versions and then I will let you know whether I get good results. – SAS Jun 28 '13 at 21:47
• @MichaelE2 the values you quote are not proper minima; you have a copy/paste error in the definition of the function. It isn't easy to find good minima for this function--using many thousands of different initial simplexes, the lowest I found was about 1500, but values around 6500 are more common. I didn't see any less than 1. I will write an answer later giving some advice on this problem... too hung over at the moment to concentrate properly. – Oleksandr R. Jun 29 '13 at 15:43
• @OleksandrR. Thanks, I appreciate your concern and waiting for your valuable suggestions . And as you mentioned, unfortunately most of the time I am also getting min which is about 6506.35 (which is too big for my purpose) – SAS Jun 29 '13 at 15:47

It seems to me that there are two possible approaches one might use to solve this problem. Either we can follow others and use the Nelder-Mead method, or we can try to use another, better suited method, such as differential evolution. This problem is very strongly multimodal, with a huge number of deceptive local minima, and this is exactly the type of function for which the Nelder-Mead method performs least well. To overcome this difficulty, it will be necessary to try the minimization repeatedly with enormously many initial simplices in the hope of falling into a favorable part of the parameter landscape.

I will show the approach using the Nelder-Mead method first, since this is what the question principally concerns. It is also the most relevant to Mathematica, since although NMinimize does support differential evolution (Method -> "DifferentialEvolution"), I actually didn't use Mathematica here, but rather my own differential evolution minimizer written in Python.

Sampling so many initial simplices will take a long time, and if we are going to minimize the same function repeatedly, it makes sense to use my compiled Nelder-Mead implementation. This will probably not be quite as fast as MINUIT or other native C++/Fortran codes, but when compiled using a C compiler, it will not be too far away from their performance, either.

First, let's write your problem down a bit more cleanly:

A = {{a11, 0}, {0, a22}}; B = {{b11, b13}, {b13, b22}};

M1 = A + Exp[I phi] B;
M2 = a A + b Exp[I theta] B;
M3 = a A - 3 b Exp[I theta] B;

m1 = FullSimplify@Abs@Eigenvalues[M1]; t1 = m1[[1]]; t2 = m1[[2]];
m2 = FullSimplify@Abs@Eigenvalues[M2]; t3 = m2[[1]]; t4 = m2[[2]];
m3 = FullSimplify@Abs@Eigenvalues[M3]; t5 = m3[[1]]; t6 = m3[[2]];

problem = Simplify[
((t1 - O1)/E1)^2 + ((t2 - O2)/E2)^2 +
((t3 - O3)/E3)^2 + ((t4 - O4)/E4)^2 +
((t5 - O5)/E5)^2 + ((t6 - O6)/E6)^2
];
parameters = {a11, a22, b11, b22, b13, a, b, phi, theta};

function = Block[{
O2 = 0.235, O1 = 72,
E2 = 0.031, E1 = 0.87,

O4 = 0.0222196, O3 = 0.958157,
E4 = 0.006198, E3 = 0.018862,

O6 = 0.0941698, O5 = 1.60087,
E6 = 3.6079*^-7, E5 = 0.00026
},
Function @@ {parameters, problem}
];


Next, we compile the minimizer with which we will attempt to solve it. Here we use a looser tolerance than usual (Sqrt[$MachineEpsilon] vs. $MachineEpsilon), since if we have already fallen into a poor local minimum there is not much point in expending further effort to find the very lowest function value in this part of parameter space. It is better to give up and try again somewhere else.

minimizer = With[{
nelderMead = NelderMeadMinimizeDumpCompiledNelderMead[function, parameters],
bounds = {-1, 1}, dimension = Length[parameters], tolerance = Sqrt[$MachineEpsilon] }, Compile[{{dummy, _Integer, 0}}, nelderMead[RandomReal[bounds, {dimension + 1, dimension}], tolerance, -1], CompilationOptions -> {"ExpressionOptimization" -> True, "InlineCompiledFunctions" -> True}, RuntimeOptions -> {"CompareWithTolerance" -> False, "EvaluateSymbolically" -> False}, RuntimeAttributes -> Listable, "Parallelization" -> True, CompilationTarget -> "C" ] ];  Now it's just a matter of running the minimization sufficiently many times to get a decent answer. For present purposes we will make only 1 000 trials, although in practice many more will be required in order to get anywhere near a function value of 1. Beware: this code is very CPU-intensive, and can freeze your computer while it runs. {First[#], Thread[parameters -> Rest[#]]} & /@ Take[minimizer@Range[1000] // Sort, 10]  I must admit, I don't know how many attempts are likely to be needed for this, as I stopped after a few hundred thousand tries without getting close to the target range. Since others claim to use this method successfully, I assume that it will eventually succeed, but this is presumably the main reason why minimizing this type of problem is reputed to be such a tedious task. Alternatively, we might expect more rapid success using differential evolution, which is much more suited to this kind of problem. The main difficulty with this method is that it contains two tuning parameters, the values of which are critical to its success or otherwise, and (unlike the Nelder-Mead parameters) cannot generally be chosen using a priori theoretical arguments. Therefore, it is usually necessary to tune these according to some sort of meta-optimization scheme. More prosaically, NMinimize is quite slow, and only implements one mutation method (the classical rand/1 operator), although other, dramatically more efficient, possibilities have appeared in the literature. For these reasons, I found it preferable to use my own code rather than Mathematica's, and this happens to be written in Python. Your function can be programmed in Python as follows: from numpy import (exp, sqrt, abs) def function(a11, a22, b11, b22, b13, a, b, phi, theta): O2 = 0.235; O1 = 72; E2 = 0.031; E1 = 0.87; O4 = 0.0222196; O3 = 0.958157; E4 = 0.006198; E3 = 0.018862; O6 = 0.0941698; O5 = 1.60087; E6 = 3.6079e-7; E5 = 0.00026; var3 = E1**2; var4 = 1/var3; var5 = -2*O1; var6 = b11 + b22; var7 = 1j*phi; var8 = exp(var7); var9 = var6*var8; var10 = -a22; var11 = a11 + var10; var12 = var11**2; var13 = -b22; var14 = b11 + var13; var15 = 2*var11*var14*var8; var16 = b13**2; var17 = 4*var16; var18 = var14**2; var19 = var17 + var18; var20 = 2j*phi; var21 = exp(var20); var22 = var19*var21; var23 = var12 + var15 + var22; var24 = sqrt(var23); var25 = -var24; var26 = a11 + a22 + var9 + var25; var27 = abs(var26); var28 = var5 + var27; var29 = var28**2; var31 = E2**2; var32 = 1/var31; var33 = -2*O2; var34 = a11 + a22 + var9 + var24; var35 = abs(var34); var36 = var33 + var35; var37 = var36**2; var39 = E3**2; var40 = 1/var39; var41 = -2*O3; var42 = a11 + a22; var43 = a*var42; var44 = 1j*theta; var45 = exp(var44); var46 = b*var6*var45; var47 = a**2; var48 = var47*var12; var49 = 2*a*var11*b*var14*var45; var50 = b**2; var51 = 2j*theta; var52 = exp(var51); var53 = var50*var19*var52; var54 = var48 + var49 + var53; var55 = sqrt(var54); var56 = -var55; var57 = var43 + var46 + var56; var58 = abs(var57); var59 = var41 + var58; var60 = var59**2; var62 = E4**2; var63 = 1/var62; var64 = -2*O4; var65 = var43 + var46 + var55; var66 = abs(var65); var67 = var64 + var66; var68 = var67**2; var70 = E5**2; var71 = 1/var70; var72 = -2*O5; var73 = -3*b*var6*var45; var74 = -6*a*var11*b*var14*var45; var75 = 9*var50*var19*var52; var76 = var48 + var74 + var75; var77 = sqrt(var76); var78 = -var77; var79 = var43 + var73 + var78; var80 = abs(var79); var81 = var72 + var80; var82 = var81**2; var84 = E6**2; var85 = 1/var84; var86 = -2*O6; var87 = var43 + var73 + var77; var88 = abs(var87); var89 = var86 + var88; var90 = var89**2; return (var4*var29 + var32*var37 + var40*var60 + var63*var68 + var71*var82 + var85*var90)/4  As you might guess, this is not hand-written code, but rather translated from the original using Mathematica. Perhaps it isn't the most efficient way to write it, but it's not too bad, taking advantage of Mathematica's simplifications and common subexpression elimination. (It may not matter very much anyway, since Python isn't a high performance language to start with. Its performance is good enough, but I use it mainly because it's widely available and very pleasant to work with.) For obvious reasons, I don't think it's a particularly good idea to post large amounts of Python code here on this site, so I'll let the definition of the function stand on its own. You can download the remaining code, i.e. the differential evolution implementation, the meta-optimization script, and an example using meta-optimized tuning parameters, from the git repository. You'll need Python 2.5 and NumPy 1.3 (or newer) to run them. If you aren't familiar with Python, it's important to note that the language has undergone some incompatible changes in Python 3, and so, although automated conversion using 2to3 should be possible, if you want to use the files unmodified, you will need to use Python 2. I should say first of all that I did try the rand/1 mutation operator, but found it ineffective for this problem. Instead, I used the MDE5 mutator described in R. Thangaraj, M. Pant, and A. Abraham, Appl. Math. Comp. 216 (2), 532 (2010), which has often proven to be a more powerful alternative to the classical method. It's possible that other schemes might be even better, but meta-optimization takes a long time, so I didn't try any others. A variety of methods are implemented in the code that I found to have advantages for particular problems, and many more are available in the literature. Before you implement any of the latter, though, you should be aware that many publications in this area fail to employ appropriate methodology or arrive at robust conclusions, so quite often a method will be described as effective that actually doesn't work well for real problems. My advice would be just to use the ones already provided, and not waste your time with this. Meta-optimization appears to be absolutely vital for good results on your particular problem, as is frequently the case in general. Here, we choose the tuning parameters$0 < F \le 2$(mixing coefficient) and$0 < C \le 1$(crossover probability) by running differential evolution on itself, aiming to minimize the objective function as robustly as possible while also limiting the number of iterations required to do so. This is somewhat difficult due to the tendency of black-box numerical minimizers to optimize exactly what you've asked for, yet not at all what you actually meant. For example, a meta-optimization scheme designed to maximize the initial rate of decrease of the function value can often result in poor performance for problems that have different characteristics at different scales. This problem becomes significantly more difficult after reaching function values of about 6000, so it's important that we also sample results from beyond this threshold. You can run the meta-optimization code (optimization_function.py) for yourself in about a day, but the results I obtained were that$F \approx 0.850$and$C = 0.975$. This is still a difficult problem, so even using these for the actual minimization (function.py) results in success (i.e., a function value far less than unity) only about 40-50% of the time. Nonetheless, this is a lot better than the Nelder-Mead results. After a run time of about two or three minutes, the result (if it succeeds) will be around$10^{-20}$, as shown by the following example solutions: solutions = { {0.006302962181370445, 0.16279423867218476, -9.548089292946795, -62.2208122770191, 24.715755884121787, -6.7915955413380225, -0.002710346806452485, -7.127353886016603*^10, 0.43823648438746504}, {-3.7513268290908126, -0.10760731343238998, -57.860129251744794, -10.517744414643426, 25.254637364611664, 0.2936728267065106, -0.003026492897390265, 0.04259960531247162, -0.5451053025455311}, {-0.2708690475702807, -69.15868531322194, -0.11698443104717468, 4.8543712324002435, 0.9233797217848526, 0.009694619653952973, -0.10947502700209849, 8.451868538832148, -4.973378159025613} }; function @@@ solutions (* -> {1.71426*10^-21, 3.81267*10^-19, 1.80781*10^-19} *)  It's interesting to note that these results are significant up to the last place--in fact, just importing the output of the Python program into Mathematica introduced a small rounding error that increased these significantly from the original values of$\approx 2 × 10^{-21}$. Here we see that rounding them off to the nearest multiple of$10^{-12}\$ can produce an increase of around 10 orders of magnitude:

function @@@ Round[solutions, 1*^-12]
(* -> {5.20078*10^-10, 2.19968*10^-10, 1.38377*10^-13} *)

• ,Wow!Differential evolution seems far better to minimize the particular function in hand.And it is a good finding that the problem is too sensitive to the number of significant digits,so I have to be really careful about it! I appreciate your help,Thanks :) – SAS Jun 30 '13 at 14:46
• @Taarchira please see the updated answer including the links to the Python code, which you are free to use as you see fit. Sorry for the delay in finishing the job. – Oleksandr R. Jul 2 '13 at 23:41
• ,Thank you for the update. – SAS Jul 6 '13 at 20:55
• @SAS yes, it's possible. But you have to build your constraints into the objective function somehow. A simple way to do it is to return a function value of infinity if the constraints are violated--but this excludes large parts of the complete search space and may not be appropriate for soft constraints. What you seek to do is not easy in a general and abstract sense, which is why I didn't address it in the algorithm itself. But for a specific problem, modifying the objective function is normally not too difficult. – Oleksandr R. Jul 12 '14 at 17:38
• @SAS One way to deal with simple bounds as constraints is to treat the boundaries as reflective, i.e. if the value lies outside of this range then it is pushed back inside according to the amount of constraint violation. This can be appropriate for some problems and inappropriate for others; you might prefer e.g. to randomly reassign values that violate constraints. Everything in global optimization is so dependent on the problem, that I'm not sure there is any general advice on this topic. – Oleksandr R. Jul 12 '14 at 17:43