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 = NelderMeadMinimize`Dump`CompiledNelderMead[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} *)
