Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a function of 20 parameters, which 3 of the parameters are my physical parameters, and the others are pull terms to fix the errors. The goal is finding the global minimum of this function, to find the best fit of those three physical parameters which are minimizing my function. I am using NMinimize, with the DifferentialEvolution method; however, choosing different options of DifferentialEvolution changes my results drastically. I really don't know if I am using those options correctly or not, and I can't be sure if they are really giving my correct minimum. How can I choose these options?: "CrossProbability", "InitialPoints", "PenaltyFunction", "PostProcess", "RandomSeed", "ScalingFactor", "SearchPoints", "Tolerance".

My code is long, I can't make it shorter, as I have to make some tables and do some Interpolation before defining my final function; however, I put my code here if somebody may need it:

d11 = 667.9; d12 = 451.8; d13 = 304.8; d14 = 336.1; d15 = 513.9; d16= 739.1;
d21 = 1556.5; d22 = 1456.2; d23 = 1395.9; d24 = 1381.3; d25 = 1413.8; d26 = 1490.1;
f11 = 0.0678; f12 = 0.1493; f13 = 0.3419; f14 = 0.2701; f15 = 0.115; f16 = 0.0558;
f21 = 0.1373; f22 = 0.1574; f23 = 0.1809; f24 = 0.1856; f25 = 0.178;f26 = 0.1608;

rhodatar = {{1.70059, 1.38938}, {1.88047, 1.24779}, {2.13609, 
1.08850}, {2.39172, 0.93805}, {2.68521, 0.76991}, {2.97870, 
0.61947}, {3.42367, 0.45133}, {3.88757, 0.30973}, {4.28521, 
0.21239}, {4.68284, 0.14159}, {5.09941, 0.08850}, {5.55385, 
0.06195}, {5.88521, 0.03540}, {6.39645, 0.01770}, {6.99290, 
0.01770}, {7.68402, 0.01770}, {8.41302, 0.00885}, {9.25562, 
0.00885}, {9.89941, 0.00885}, {10.89941, 0.00885}, {12., 0.00885}};

rhor = Interpolation[rhodatar];

rhofinalr[x_] := rhor[x]/NIntegrate[rhor[x], {x, 1.8, 12}];

sterm2ofp11 =NIntegrate[Sin[1.267*2.32*10^-3*d11/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f11sin1 =Table[{w,NIntegrate[Sin[1.267*w*d11/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]},{w, 
Table[10^w, {w, -2.634512, -0.5, 0.01}]}];sterm3ofp11 = Interpolation[f11sin1];
f11sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d11/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm4ofp11 = Interpolation[f11sin2];
p11n =.;
p11n[y_, z_, w_] :=f11*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp11 - ((1 + Sqrt[1 y])/2)*z*sterm3ofp11[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp11[w]);


sterm2ofp12 =NIntegrate[Sin[1.267*2.32*10^-3*d12/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f12sin1 = Table[{w,NIntegrate[Sin[1.267*w*d12/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}{w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm3ofp12 = Interpolation[f12sin1];
f12sin2 = Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d12/enu]^2*rhofinalr[enu],{enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm4ofp12 = Interpolation[f12sin2];
p12n =.;
p12n[y_, z_, w_] :=f12*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp12 - ((1 + Sqrt[1 - y])/2)*z*sterm3ofp12[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp12[w]);


sterm2ofp13 =NIntegrate[Sin[1.267*2.32*10^-3*d13/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f13sin1=Table[{w,NIntegrate[Sin[1.267*w*d13/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm3ofp13 = Interpolation[f13sin1];
f13sin2=Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d13/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm4ofp13 = Interpolation[f13sin2];
p13n =.;
p13n[y_, z_, w_] :=f13*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp13 - ((1 + Sqrt[1-y])/2)*z*sterm3ofp13[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp13[w]);


sterm2ofp14=NIntegrate[Sin[1.267*2.32*10^-3*d14/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f14sin1=Table[{w,NIntegrate[Sin[1.267*w*d14/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm3ofp14 = Interpolation[f14sin1];
f14sin2=Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d14/enu]^2*rhofinalr[enu], {enu, 1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm4ofp14 = Interpolation[f14sin2];
p14n =.;
p14n[y_, z_, w_]:=f14*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp14 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp14[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp14[w]);


sterm2ofp15=NIntegrate[Sin[1.267*2.32*10^-3*d15/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f15sin1=Table[{w,NIntegrate[Sin[1.267*w*d15/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm3ofp15 = Interpolation[f15sin1];
f15sin2=Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d15/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp15 = Interpolation[f15sin2];
p15n =.;
p15n[y_, z_, w_] :=f15*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp15 - ((1 + Sqrt[1 - y])/2)*z*sterm3ofp15[w] - ((1 - Sqrt[y])/2)*z*sterm4ofp15[w]);


sterm2ofp16=NIntegrate[Sin[1.267*2.32*10^-3*d16/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f16sin1=Table[{w,NIntegrate[Sin[1.267*w*d16/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}{w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}];
sterm3ofp16 = Interpolation[f16sin1];
f16sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d16/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp16 = Interpolation[f16sin2];
p16n =.;
p16n[y_, z_, w_] :=f16*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp16 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp16[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp16[w]);


sterm2ofp21=NIntegrate[Sin[1.267*2.32*10^-3*d21/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f21sin1=Table[{w,NIntegrate[Sin[1.267*w*d21/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}]  
sterm3ofp21 = Interpolation[f21sin1];
f21sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d21/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp21 = Interpolation[f21sin2];
p21n =.;
p21n[y_, z_, w_] :=f21*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp21 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp21[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp21[w]);

sterm2ofp22=NIntegrate[Sin[1.267*2.32*10^-3*d22/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f22sin1=Table[{w,NIntegrate[Sin[1.267*w*d22/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}]  
sterm3ofp22 = Interpolation[f22sin1];
f22sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d22/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp22 = Interpolation[f22sin2];
p22n =.;
p22n[y_, z_, w_] :=f22*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp22 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp22[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp22[w]);

sterm2ofp23=NIntegrate[Sin[1.267*2.32*10^-3*d23/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f23sin1=Table[{w,NIntegrate[Sin[1.267*w*d23/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}]  
sterm3ofp23 = Interpolation[f23sin1];
f23sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d23/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp23 = Interpolation[f23sin2];
p23n =.;
p23n[y_, z_, w_] :=f23*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp23 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp23[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp23[w]);


sterm2ofp24=NIntegrate[Sin[1.267*2.32*10^-3*d24/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f24sin1=Table[{w,NIntegrate[Sin[1.267*w*d24/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}]  
sterm3ofp24 = Interpolation[f24sin1];
f24sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d24/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp24 = Interpolation[f24sin2];
p24n =.;
p24n[y_, z_, w_] :=f24*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp24 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp24[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp24[w]);


sterm2ofp25=NIntegrate[Sin[1.267*2.32*10^-3*d25/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f25sin1=Table[{w,NIntegrate[Sin[1.267*w*d25/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}]  
sterm3ofp25 = Interpolation[f21sin1];
f25sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d25/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp25 = Interpolation[f25sin2];
p25n =.;
p25n[y_, z_, w_] :=f25*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp25 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp25[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp25[w]);


sterm2ofp26=NIntegrate[Sin[1.267*2.32*10^-3*d26/enu]^2*rhofinalr[enu], {enu, 1.8, 12}];
f26sin1=Table[{w,NIntegrate[Sin[1.267*w*d26/enu]^2*rhofinalr[enu], {enu, 1.8, 12}]}, {w,Table[10^w, {w, -2.634512, -0.5, 0.01}]}]  
sterm3ofp26 = Interpolation[f26sin1];
f26sin2 =Table[{w,NIntegrate[Sin[1.267*(w - 2.32*10^-3)*d26/enu]^2*rhofinalr[enu], {enu,1.8,12}]}, {w, Table[10^w, {w, -2.634512, -0.5, 0.01}]}]; 
sterm4ofp26 = Interpolation[f26sin2];
p26n =.;
p26n[y_, z_, w_] :=f26*(1 - y*((1 + Sqrt[1 - z])/2)^2*sterm2ofp26 - ((1 + Sqrt[1 -y])/2)*z*sterm3ofp26[w] - ((1 - Sqrt[1 - y])/2)*z*sterm4ofp26[w]);

normfactorfar = 17566.8; normfactornear = 151725.5;

obsnear = 149904.8; obsfar = 16161.5;

chi2reno[y_, z_, w_, a_, xinear_, fnear1_, fnear2_, fnear3_, fnear4_,fnear5_, fnear6_,xifar_, ffar1_, ffar2_, ffar3_, ffar4_, ffar5_, ffar6_, bnear_,bfar_] :=
(obsnear + bnear -normfactornear*(1 + a +xinear)*((1 + fnear1)*(p11n[y, z, w]) + 
(1 + fnear2)*(p12n[y, z, w]) + (1 + fnear3)*(p13n[y, z, w]) + 
(1 +fnear4)*(p14n[y, z, w]) + (1 + fnear5)*(p15n[y, z,w]) + 
(1 + fnear6)*(p16n[y, z, w])))^2/obsnear + 
(fnear1^2 + fnear2^2 + fnear3^2 + fnear4^2 +fnear5^2 +fnear6^2)/(0.009)^2 +xinear^2/(0.002)^2 +bnear^2/(1140.93)^2 + 
(obsfar + bfar-normfactorfar*(1 + a +xifar)*((1 + ffar1)*(p21n[y, z, w]) + 
(1 + ffar2)*(p22n[y,z, w]) +(1 + ffar3)*(p23n[y, z,w]) + 
(1 + ffar4)*(p24n[y, z, w]) + (1 + ffar5)*(p25n[y, z, w]) + 
(1 +ffar6)*(p26n[y, z, w])))^2/obsfar + 
(ffar1^2 + ffar2^2 + ffar3^2 + ffar4^2 + ffar5^2 + ffar6^2)/(0.009)^2 + xifar^2/(0.002)^2 + bfar^2/(166.545)^2;

renovars = {y, z, w, a, xinear, fnear1, fnear2, fnear3, fnear4,fnear5, fnear6, xifar,ffar1, ffar2, ffar3, ffar4, ffar5, ffar6,bnear, bfar};

renobounds = {0. <= y <= 1, 0 <= z <= 1, 0.00232 <= w <= 0.1,bnear >= 0, bfar >= 0};

Changing the values of these options, give me very different results:

Do[Print[NMinimize[{chi2reno[y, z, w, a, xinear, fnear1, fnear2, 
fnear3, fnear4, fnear5, fnear6, xifar, ffar1, ffar2, ffar3, 
ffar4, ffar5, ffar6, bnear, bfar], renobounds}, renovars, 
Method -> {"DifferentialEvolution", "SearchPoints" -> Automatic, 
"ScalingFactor" -> 0.9, "CrossProbability" -> 0.1, 
"PostProcess" -> {FindMinimum, Method -> "QuasiNewton"},
"RandomSeed" -> i}]], {i, 10}]
share|improve this question
    
@OleksandrR, I used the methods that you gave me the link, however, I have this problem that I have asked in the zbove question. Do you know how can I solve it? –  ZKT Mar 14 '13 at 6:03
    
Hello and welcome to the site. First of all, you have commas missing in your definitions of f12sin1 and f16sin1. Global minimization is a difficult task, and differential evolution, being a heuristic, can't give you strong performance guarantees. Appropriate values of the scaling factor $F$ and the crossover probability $C$ are strongly problem-dependent, and choosing these is arguably beyond Mathematica's scope--you should refer to the literature on differential evolution. One way to tune the parameters is to choose an easier model problem that shares some of the characteristics ... –  Oleksandr R. Mar 14 '13 at 6:15
    
... of the real one, and numerically minimize the result of the minimization of this function with respect to the parameter values. I have done this with "SearchPoints" -> 60 for the (20-dimensional) Rastrigin's function ("ScalingFactor" -> 0.2, "CrossProbability" -> 0.6) and the (20-dimensional) Rosenbrock's function ("ScalingFactor" -> 0.6, "CrossProbability" -> 0.1). Whether these results can be of use to you I don't know. (Incidentally, I didn't use Mathematica for this meta-optimization.) –  Oleksandr R. Mar 14 '13 at 6:19
    
Sorry, I didn't recognise your username from our previous discussion, or otherwise I wouldn't have welcomed you to the site again. :) Perhaps you'd like to choose a more memorable name? As another member once commented, this gives the site less of the feel of a prison and more of a happy community! –  Oleksandr R. Mar 14 '13 at 6:23
    
Dear OleksandrR, thanks a lot for your answer. I chose DifferentialEvolution because I thought it is a better way for finding global minimum, the other methods don't seem to work well. I am dealling with real experimental data in my problem, and to be able to compare my results with the ones given by them, I have to use exactly the same fuction they have defined... –  ZKT Mar 14 '13 at 6:28
show 8 more comments

closed as too localized by Oleksandr R., m_goldberg, Sjoerd C. de Vries, halirutan, rm -rf Mar 17 '13 at 23:44

This question is unlikely to help any future visitors; it is only relevant to a small geographic area, a specific moment in time, or an extraordinarily narrow situation that is not generally applicable to the worldwide audience of the internet. For help making this question more broadly applicable, visit the help center.If this question can be reworded to fit the rules in the help center, please edit the question.

1 Answer

[Not really a solid answer but more code than I want in a comment.]

I think your function is just very difficult to work with. It may have a very flat landscape but I suspect, more likely, it jumpt quite a bit and is numerically sensitive to any number of things. I notice that Interpolation was used to construct it and that can give awkward wiggles. Possibly using Method->"Spline" might help if that is causing trouble.

I do a couple of things that seem to help slightly. One is to use large values for iterations and generation sizes. Another is to pull out all the stops on postprocessing.

This run gives results that tend to be consistently in the range of 0.1-0.3. I believe I have seen results much closer to zero so this might not indicate best possible, but at least they are not giving results like 1000.

Timing[Do[
  Print[Timing[
    NMinimize[{chi2reno[y, z, w, a, xinear, fnear1, fnear2, fnear3, 
       fnear4, fnear5, fnear6, xifar, ffar1, ffar2, ffar3, ffar4, 
       ffar5, ffar6, bnear, bfar], renobounds}, renovars, 
     MaxIterations -> 2000, 
     Method -> {"DifferentialEvolution", 
       "SearchPoints" -> 200,
       "PostProcess" -> {FindMinimum, 
         Method -> {"QuasiNewton", "InteriorPoint", "KKT"}}, 
       "RandomSeed" -> i}]]], {i, 10}]]

Certain variables tend to be near zero in most or all cases. Others vary by a fair amount. This might just be a very sensitive function.

--- edit ---

As a general remark, this example seems to show pernicious effects from the bnear and bfar variables, as results have them varying wildly whereas all others, best I can tell, seem to stay in teh same region from one result to another.

--- end edit ---

share|improve this answer
    
+1; I agree with your assessment although the question should perhaps be closed as too localized or not answerable (through no fault of the OP; this is a general problem). Some comments follow. A large population is seldom useful for differential evolution except in the case of high dimensional functions (contrary to the documentation and the early literature on the method, $NP \approx d$ is better than $NP \gg d$ except where $d$ is so low that this would lead to an impractically small number of search points), and for 20-dimensional functions I find that $NP \approx 40-60$ ... –  Oleksandr R. Mar 16 '13 at 21:42
    
... is sufficient; 80 is ample and 200 is certainly too large. The main problem is tuning $F$ and $CR$. "SearchPoints" -> 80, "ScalingFactor" -> 0.55, "CrossProbability" -> 0.05 works quite well for this problem, though these values are meta-optimized based on the Rosenbrock's function. The function value with these settings usually turns out around $10^{-10}$. I wonder though, could you elaborate on what it means to specify more than one Method option for FindMinimum here? In particular, I'm unsure what relationship "QuasiNewton", "InteriorPoint", and "KKT" have to each other? –  Oleksandr R. Mar 16 '13 at 21:45
    
@Oleksandr R. If I recall correctly, it means the postprocessing step will try all three to home in on a best result. –  Daniel Lichtblau Mar 16 '13 at 21:56
    
@Oleksandr R. Re "SearchPoints" settings: my experience has been that if one is going to do many iterations, fairly high settings for this parameter help to ward off premature convergence. –  Daniel Lichtblau Mar 16 '13 at 22:18
    
@Oleksandr R. Re "CrossProbability" and, to a lesser extent, "ScalingFactor" settings: I rarely have consistent luck with these except when some or all variables are integer valued. Your settings in this case do seem to give an improvement. That said, i will point out that it is inconsistent, and I have seen some amount of fluctuation even with those settings (for example, where using twice as many search points gives a significantly worse claimed optimum). –  Daniel Lichtblau Mar 16 '13 at 22:24
show 6 more comments

Not the answer you're looking for? Browse other questions tagged or ask your own question.