Hello my problem is simple but I am not an expert in statistics to figure out the solution to this. Basically I have a model which takes in 3 parameters (A,T,D) and then spits out a SCORE. My goal is to find the best values for the 3 parameters which gives the highest score. The difficulty of course is that I don't know the model and it's likely not linear, so derivatives and partial derivatives based solutions will not work. There can also be many local maximums so we could easily get ambushed there, so it can get very tricky, I want the global maximum for SCORE or at least an estimation of it which can be computed in an affordable time.

So I treat the model as a black box, 3 parameters go in (A,T,D) 1 comes out (SCORE) and the task is to find the 3 parameters which give the best score. The parameters luckily are constrained and not infinite, I have normalized them, so they all take values between (0.0,1.0] and the output score can be between 0 and 6200, it doesn't have an upper limit but it's extremely unlikely to go above 6200.

What I did so far is just draw random input values, basically random uniform distribution values for A,T,D and then manually check how they change with the SCORE value.

One solution here is just to draw as many samples and then pick the one with the highest SCORE and call it a day but this is not good and computationally very inefficient because if 99% of the sample is rubbish and the good values are all concentrated around the other 1% then we waste a lot of resources computing bad values.

So what I did is draw a moderate amount, sort the output by the SCORE in ascending order and then analyze both visually and analytically.

For example here is a random sample of n=10, ordered in ascending order:

enter image description here

(of course normally I would draw at least 10,000 samples)

And then I cutoff the parameter's range based on where it's likely to give rubbish results. I just take the average of the SCORE column and compute the Pearson correlation for each variable column and based on these 2 metrics I try to find the range that likely contains the best values.

For example the correlations here are:




This one suggests that the A variable range has to be increased to (0.0+x,1.0], the T variable is indifferent so then I just take the average of the column which is 0.682599909364871 and then the T parameter can be fixed to this and only work with A and B, and the B parameter also increased by (0.0+y,1.0], if the correl had been negative then obviously it would have been (0.0,1.0-y]. The X and Y value I try to gauge manually based on how big of a value I need to increase the average of the SCORE column. So if I go carefully with 0.01 increments and until the average of the SCORE column increases as I constrain the ranges, that means that I am heading in the right direction and when it doesn't then I stop. And then I draw another sample stack this time with only 2 parameters , and repeat this until all the correlations converge towards 0. And then finally I am left with 3 values for my input parameters which give the best SCORE that I can tell, but there is of course no way to verify this because it's completely arbitrary and manual.

It is certainly more efficient than just going through all the values which are infinite as the best value might hide between say 0.0000001 and 0.0000002 which would take 10^21 calculations which is extremely bad, so doing step by step constraints is better, but there is a risk of overshooting or undershooting and all of it is based mostly on instincts with the help of correlation and average to guide me.

Now here is my question, is there a formal and more analytical way of doing this automatically, but without errors and for me to be able to set the threshold for tolerance and the number of digits of precision I would work with.

How to find the best (highest) output for a black box model with 3 inputs in the most efficient way, where there is a relationship between the input and the output (so the output is not an independent random variable) but it can't be described by a simple formula (so no derivatives)?

  • $\begingroup$ Most likely there is no really good answer to this. A lot depends on 1) whether you have some a priori information about your "black box" and 2) what exactly you mean by an "efficient computation". Look at proceedings.mlr.press/v70/malherbe17a/malherbe17a.pdf, for instance. They try to minimize the number of function calls but do not care much about the time needed to choose the next evaluation point and the whole idea is just to ignore the regions in which a low value is observed. That's not very impressive but, in general, it is not clear if one can do any better. $\endgroup$ – fedja Jan 11 '20 at 13:23
  • $\begingroup$ @fedja 1) I don't have a priori information hence the black box. I only know the probable minimum and maximum range for each param, but within that haystack nothing, I don't know whether there is a smooth tendency of any kind of if it's just fragments of local maximums randomly scattered, most likely the second. $\endgroup$ – user421473 Jan 13 '20 at 12:42
  • $\begingroup$ @fedja 2) Unfortunately I don't like linear searchers, they all get stuck in stationary zig-zags, at which point you need to set a tolerance to jump over them, which at that point becomes no better than random search and in multidimensional cases it becomes very slow. However random search is too ignorant about local patterns so I'd rather use a genetic algorithm then, which I already have, but it needs to be optimized too, that is what I am trying to do here. It's quite annoying that you need to optimize the optimizer which then becomes an infinite regression. $\endgroup$ – user421473 Jan 13 '20 at 12:45
  • $\begingroup$ OK, I'll try a couple of ideas and see if they lead anywhere. Meanwhile, you can evaluate the efficiency of your current algorithm on some randomly generated functions, say $f(X)=\max_{i=1}^M[V_i-s|X-X_i|]$ where $X_i$ are chosen randomly in the unit cube and $V_i$ are randomly chosen in $[0,1]$, say, with varying $M$ and $s$. Since the maximum is attained at one of $X_i$, you can compute it but let the machine keep it secret from you until you finish your algorithm. Then compare the value you got with the truth and see what the success rate is. I wonder if you can do better than 80% :-) $\endgroup$ – fedja Jan 13 '20 at 13:00
  • $\begingroup$ Just curious (to compare what I'm doing to what you have): what is the total number of test points you think you can afford? (I'll use the same number and try to optimize the scheme within that constraint). $\endgroup$ – fedja Jan 13 '20 at 22:10

Edit2 One more code (this time in C++; if you get interested, I can try to convert to Python, but otherwise I would not bother: Python is a good language, but the idea to count the number of spaces to check the block structure when debugging the code drives me nuts). This code can handle maxima like in example 2 (within reasonable limits, of course) and it tries to learn the shape of the level sets on the go. Again, I bounded the number of calls by 20000 and aimed at 60% success rate. If you can afford more calls, let me know: that should increase the performance quite a bit.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

const long int N=30, TIME=20000, K=18, M=55, NN=50, D=3;
int cwrite=1000;
double powmean=1.3, powspan=0.8, sc=2*exp(-log(N)/3), speedadj=0.1, matrspeed=20;
double SETTIME=3*N/speedadj, REL=2;
double ssl=5;

double X[M][3], V[M], cf[M][D][3], pp[M][D], sl[M]; 

double ff(double t)
return t*t;

double gg(double v)
return (v>0?v*v:-v*v);  

double abs(double x)
return (x>0?x:-x);

double Norm(double x[3])
double s2=0;
for(int i=0;i<3;++i) s2+=x[i]*x[i];
return sqrt(s2);

void update(double X[][3], double eps, double unit[3])
double pr[3][3], Y[3][3];
for(int i=0;i<3;++i) for(int j=0; j<3;++j) pr[i][j]=(i==j?1:0)-(1-exp(-eps))*unit[i]*unit[j]/Norm(unit);
for(int i=0;i<3;++i) for(int j=0; j<3;++j)
Y[i][j]=0; for(int k=0;k<3;++k) Y[i][j]+=X[i][k]*pr[k][j];
for(int i=0;i<3;++i) for(int j=0; j<3;++j) X[i][j]=Y[i][j];

double HS(double X[][3])
double hs2=0;
for(int i=0;i<3;++i) for(int j=0; j<3;++j) hs2+=X[i][j]*X[i][j];
return sqrt(hs2/3);

double dist(double x[3], double y[3])
double s2=0;
for(int i=0;i<3;++i) s2+=(x[i]-y[i])*(x[i]-y[i]);
return sqrt(s2);

double unit[3];

double unitrand()
return ((double)rand())/RAND_MAX;

void brand()
for(int i=0;i<3;++i) unit[i]=2*unitrand()-1;
while(Norm(unit)>1) {for(int i=0;i<3;++i) unit[i]=2*unitrand()-1;}

void crand()
for(int i=0;i<3;++i) unit[i]=unitrand();

void proj(double x[3])
for(int i=0;i<3;++i) {if(x[i]<0) x[i]=0; if (x[i]>1) x[i]=1;}

double Log(double x)
double y=(x>0?x:-x)+1.0/100000000000000000;
return log(y);

double dot(double x[3], double y[3])
double s=0;
for(int i=0;i<3;++i) s+=x[i]*y[i];
return s;

void assign(double x[3], double y[3])
for(int i=0;i<3;++i) x[i]=y[i];

void mult(double S[3][3], double x[3])
double y[3];
for(int i=0;i<3;++i) y[i]=dot(S[i],x);

void add(double x[3],double y[3])
for(int i=0;i<3;++i) x[i]+=y[i];

void multnum(double S[3][3], double fact)
for(int i=0;i<3;++i) for(int j=0;j<3;++j) S[i][j]*=fact;

double decr(double sl, double cf[D][3], double pp[D], double x[3])
double s=0;
for(int i=0;i<D;++i) s+=sl*exp(pp[i]*Log(dot(cf[i],x)));
return s;   

struct T{double x[3]; double v; double speed; int succ; double S[3][3]; T *t;}; 
int NNN;

void del(T *P)
P->t=P->t->t; --NNN; 

double f(double x[3])
double s=0;
for(int i=0;i<M;++i)
double y[3]; for(int j=0;j<3;++j) y[j]=x[j]-X[i][j];    
double val=V[i]-decr(sl[i],cf[i],pp[i],y); 
if(i==0 || val>s) s=val;
return s;

int main()


int count=0, ccount=0;

for(int k=0;k<NN;++k)

for(int i=0;i<M;++i)
crand(); assign(X[i],unit); 
for(int j=0;j<D;++j) {brand(); assign(cf[i][j],unit); pp[i][j]=powmean+powspan*(2*unitrand()-1);} 

int ti;
double tv;

ti=0; tv=V[0];
for(int i=0;i<M;++i) {double vv=V[i]; if(vv>tv) {ti=i; tv=vv;};}

T *P,*Q,*QQ; 

double xmax[3]; double vmax;

P=new T; crand(); assign(P->x,unit); P->v=f(P->x);  
P->speed=0; P->succ=0; 
for(int i=0;i<3;++i) for(int j=0;j<3;++j) P->S[i][j]=(i==j?sc:0);
assign(xmax,P->x); vmax=P->v; 

Q=P; QQ=P;

for(int n=1;n<N;++n) 
Q->t=new T; Q=Q->t; 
crand(); assign(Q->x,unit); Q->v=f(Q->x);  
Q->speed=0; Q->succ=0; 
for(int i=0;i<3;++i) for(int j=0;j<3;++j) Q->S[i][j]=(i==j?sc:0);
if(Q->v>vmax){assign(xmax,Q->x); vmax=Q->v;} 

for(long int ttime=0;ttime<TIME;++ttime)
T *W=Q->t; 

double shift[3], x[3];
brand(); assign(x,W->x); assign(shift,unit); mult(W->S,shift); add(x,shift); proj(x); double v=f(x);

double maxx=(W->speed>abs(v-W->v)?W->speed:abs(v-W->v));

double newsc=sc*exp(-K*ff((double)ttime/TIME));
double expfact=K*(ff((double)(ttime+NNN)/TIME)-ff((double)ttime/TIME)); 

double relspeed=(W->speed>0?(W->v-v)/W->speed:0); 
if(relspeed>1) relspeed=1;
if(relspeed<-1) relspeed=-1;

double matr=matrspeed*expfact; 
double eps=matr*gg(relspeed);
if(abs(eps)>2) {eps=2*eps/abs(eps);} eps=eps*(4-abs(eps))/2; 

T *WW=W; double ddist=dist(x,W->x);
for(int qq=0;qq<NNN;++qq)
if(dist(x,W->x)<ddist && W->v<v) {WW=W; ddist=dist(x,W->x);}
assign(WW->x,x); WW->v=v; ++WW->succ;
if(v>vmax) { assign(xmax,x); vmax=v; QQ=WW;}

double oldhs=HS(W->S);
if(ttime>SETTIME) update(W->S,eps,unit); 
double fact=HS(W->S)/oldhs, corrfact=newsc/HS(W->S); 


if(ttime>SETTIME && W->v+REL*W->speed<vmax) del(Q); else Q=Q->t;

printf("true max\n");
printf("%.12f    %.12f    %.12f   %.12f\n",X[ti][0],X[ti][1],X[ti][2],tv);

printf("found with %d candidates left\n", NNN);
printf("%.12f    %.12f    %.12f      %.12f    %d\n",QQ->x[0],QQ->x[1],QQ->x[2], QQ->v, QQ->succ);
if(dist(xmax,X[ti])<0.000001 || tv-vmax<0.0000001) {++count; ++ccount; printf("success!\n");} 
double ddist=3;
for(int i=0;i<M;++i) if(ddist>dist(QQ->x,X[i])) ddist=dist(QQ->x,X[i]); 
if (ddist<0.00001) ++ccount;
printf("distance to a local max= %.12f\n",ddist);


printf("success rates %.2f  %.2f", (double)count/NN, (double)ccount/NN);    

return 0;


Let me try to explain to you the general idea about random search algorithms as I understand it myself (I'm by no means an expert, but after looking at a few papers in the subject like the one I quoted, I suspect that I'm not a total ignoramus either).

The general algorithm that you want to finish in $T$ steps with (hopefully) an error $exp(-K)$ with one running point $x$ and a scale $s$ associated with is as follows. You initialize the scale $s$ to $1$, $x$ to a random point in the unit cube, and your supposed maximum to $v=f(x)$. At step $t$ you choose a random point $y$ in the cube centered at $x$ of radius $s$ (I assume the $\ell^\infty$ metric), compute $w=f(y)$, and if $w>v$, replace $x$ by $y$ and $v$ by $w$. Otherwise you keep them as they were. After that you reduce your scale to $s=\exp(-K\varphi(t/T))$ where $\varphi$ is some increasing bijection from $[0,1]$ to itself (the most obvious regime to try is $\varphi(t)=t$ but you may want to go a bit faster in the beginning and a bit slower near the end for the reasons to be explained below. Notice that by time $t=T$ you'll be searching the area of size $e^{-K}$ no matter what. The question is why you have some reason to believe that your maximum is still there.

The story is that if $x_0$ is the true global maximum, then the level sets $\{f\ge eta\}$ have certain volume $V(eta)$ and diameter $D(\eta)$. You also have to assume that the maximum is not too weird, so when eta approaches the maximum, both $V(\eta)$ and $D(\eta)$ to $0$ and there is some nice relation between them. Let's look at two examples:

Example 1:

$f(x)=-a|x|-b|y|-c|z|$. Then the level sets are octahedra and $V(\eta)\approx |\eta|^{3}(abc)^{-1}$, while $D(\eta)\approx |\eta|\min(a,b,c)^{-1}$, so in this case $V$ is proportional to $D^3$ (which is good) but the coefficient of the proportionality is about $\frac{\min(a,b,c)^3}{abc}$, which can be pretty small if $a,b,c$ are different from each other.

Example 2:

$f(x)=-|x|^3-y^2-|z|$. Then the level sets are still octahedra-like but their shape changes when the level increases. Here $V(\eta)= |\eta|^{11/6}$, while $D(\eta)\approx |\eta|^{1/3}$, so in this case $V$ is proportional to $D^{11/2}$ (which is bad, as we will see).

If you are at some point $x$ with scale $s$ and value $f(x)=v$, your hope is that $x_0$ is within distance $s$ away, i.e., that the level set $\{f>v\}$ satisfies $D<s$. At this point you are searching in a cube of size $s$, so the probability to hit that level set is about $V/s^3$. Once you approach $s\approx D$ when you descend, you have the chance $V/D^3$ to hit it, so you'd better stay at that size of $s$ at least $D^3/V$ steps to get any chance to improve before you reduce your scale by some noticeable factor (say, $2$). The good news is that if the hit is scored, then the distribution over the level set is still uniform, so (unless the level set shape is really weird) the level set shrinks a constant number of times, so the reduction of scale after the hit is justified and you are again in a good situation to start it all over.

Notice that in example 1, the ratio $V/D^3$ stays fixed but in example $2$ it is about $D^{5/2}$, which means that in the first example it is enough to spend some constant (though, perhaps, big) time at scale $D$ while in the second example the time you need to spend at that scale is about $D^{-5/2}$, so if you aim at the precision $e^{-K}$ and time $T$ and choose the linear regime, then you'd better have $T/K>e^{\frac 52 K}$ (roughly speaking) or you will be destined to lose your grip on the maximum at some point (you'll still get something, of course, but certainly not the precision you are aiming at).

This is the one point story and the main moral of it is "Go down as slowly as you can afford to still have a chance to reach your destination and, if you can, decelerate towards the end". If you try to accelerate, you'll just end far from your destination.

One point search is not too bad especially if you run it a few times independently (at each run there is a chance that you'll end up at some local max instead) but you can try to do a bit better with running several points at once and excluding the hopeless candidates (which will have the desirable effect of deceleration of the scale reduction for the winners). The implementation of the one point algorithm is trivial, so just try it and see what it gives. Of course, you want to try things on the models you create yourself rather than on the ones for which you have no idea what the truth is to evaluate the efficiency, so think what would be a good family to consider as the testing ground. I suggested one, but you may want to choose some other one if you believe that my suggestion is too simplistic (example 2 certainly doesn't show there).

Let me know if you are interested in continuing this discussion or if you knew all I said already :-).

Since you didn't respond, I took the liberty to formalize the setup myself. The notion of an "arbitrary function" is not well defined, so I considered the family of random functions $f(x)=\max_{i=1}^M [V_i-\|A_i(x-x_i)\|_1]$ where $V_i$ are independent random variables uniformly distributed in $[0,1]$, $x_i$ are random points in the cube, and $A_i$ are random linear operators (within reasonable limits) and the norm is the $\ell^1$ norm (so the function is piecewise linear). This generates something with line cross-sections looking like this (for $M=55$):

enter image description here

I hope it looks "arbitrary" enough. The program (in Asymptote) is attached. It should be more or less self-explanatory. The success rate (finding the maximum with more or less machine precision by throwing random points in certain changing regions) is about 60% for $M=55$. It is still short of 80% but you shouldn't expect miracles :-).

I decided that once you have 10,000 samples for preliminary analysis only, I can afford 20,000 function calls with constant time of call processing. Try it on your data and let me know if it works to your satisfaction (why it works if it does is another story; we can discuss it if you find it useful). If somebody can do a better job on this model by pure random selection without heavy reliance on the piecewise linear structure, I'll be most interested.

import three;

///////////// Parameters to play with (I hope, nearly optimal for 20000 function calls) ///////

int N=64, TIME=20000, K=10; 
real speedadj=0.1, sc=1;
real SETTIME=2*N/speedadj, REL=5;

//////////// End of parameters //////////////////

int NN=100; //Number of trials to determine efficiency
int count=0; //Just the counter for successes

struct T{triple x; real v; real sc; real speed; int succ; T t;} /// The list of running points type ////
T P,Q,QQ; 

/////// Stupid nonsense for generating points in the cube and such ////////////

triple trand()
triple a=2*(unitrand(),unitrand(),unitrand())-(1,1,1);
return a;

real cut(real x)
if(x<0) return 0;
if(x>1) return 1;
return x;

triple proj(triple a)
return (cut(a.x),cut(a.y),cut(a.z));

/// End stupid nonsense ////

for(int k=0;k<NN;++k)

int sec=seconds()+k; //You can set it to a reported bad seed to see if the corresponding profile is really problematic

int NNN=N; // The number of moving test points

int time=0; // The function call counter

/////////////////////// Generating a random function /////////////////////////////////

int M=55; //The number of potential local maxima 
real sl=2; //The steepness of humps

triple[] X; real[] V,s;

for(int i=0;i<M;++i) {X[i]=(0.5,0.5,0.5)+0.5*trand(); V[i]=unitrand(); s[i]=sl*(1+2*unitrand());}

typedef real F(triple);
F[] f;

for(int i=0;i<M;++i) 
triple cfx=trand(), cfy=trand(), cfz=trand();  
f[i]=new real(triple a){return abs(dot(cfx,a))+abs(dot(cfy,a))+abs(dot(cfz,a));};

real f(triple x)
real v=V[0]-s[0]*f[0](x-X[0]);
for(int i=1;i<M;++i) {real vv=V[i]-s[i]*f[i](x-X[i]); if(vv>v) v=vv;}
return v;

int cmax=0;
int ti=0; real tv=f(X[0]);
for(int i=0;i<M;++i) {real vv=f(X[i]); if(vv>tv) {ti=i; tv=vv;}; if(vv==V[i]) ++cmax;}

write(string(cmax)+" maxima");
write("true max");

//s[ti]*=2; // Try to make the maximal peak narrower, if you want, to see how it all breaks down

/////////////////// End of generating a random function /////////////////////////

srand(k^3); // Allows to randomly treat the same function many times ///

////////////////// The algorithm body ////////////////// 

void delete(T P)
P.t=P.t.t; --NNN; 

triple xmax; real vmax;

P=new T; P.x=(0.5,0.5,0.5)+0.5*trand(); P.v=f(P.x); P.sc=sc; 
P.speed=0; P.succ=0; 
xmax=P.x; vmax=P.v; 

Q=P; QQ=P;

for(int n=1;n<N;++n) 
Q.t=new T; Q=Q.t; Q.x=(0.5,0.5,0.5)+0.5*trand(); Q.v=f(Q.x); Q.sc=sc; 
Q.speed=0; Q.succ=0; 
if(Q.v>vmax) {xmax=Q.x; vmax=Q.v; QQ=Q;}

T W=Q.t; 

triple x=proj(W.x+W.sc*trand()); real v=f(x);
if(v>vmax) {xmax=x; vmax=v; QQ=W;}
if(v>W.v){W.x=x; W.v=v; ++W.succ;}

if(time>SETTIME && W.v+REL*W.speed<vmax) delete(Q); else Q=Q.t;

///////////////// Report of success/failure for a single run //////////////////

write("found with "+string(NNN)+" candidates left");
if(abs(xmax-X[ti])<0.000001 && tv-vmax<0.000001) {++count; write("success!");} 
else {write("failure at seed "+string(sec));}

/////////////// End single run report ////////////////////


///////// End of the algorithm body /////////////////////

write("success rate is ", suffix=none); write(count/NN); //Just the overall statistics
  • $\begingroup$ The random function is good enough, hopefully you only ran it on 1 set of outputs, and then if affordable on other ones, if not then selecting an output with many local zigzags (gray line on the chart) will test it better. Now as for your algo, I am not familiar with Asymptote lang, I am a Python person, but I see your update code Q.x=(0.5,0.5,0.5)+0.5*trand() is horribly inefficient, you are just randomly jumping all over the place with a fixed threshold of 0.5 which is not optimized either. On one hand you jump over ground that you already covered with no logic in it. $\endgroup$ – user421473 Jan 21 '20 at 18:16
  • $\begingroup$ Since if the function is truly random the maximum can be anywhere (and can repeat) so focusing on 1 area is better and not go over the same place that you already covered. On the other hand this is still fixed jumping with a fixed threshold, the threshold has to change over time too. You could also keep track of the rate of change, since if the function is random ,it doesn't matter, but if it isn't then it does help zoom into interesting areas. So this algo is probably equivalent in success rate to a particle swarm algo, but more efficient, but less efficient than a genetic algo. $\endgroup$ – user421473 Jan 21 '20 at 18:21
  • $\begingroup$ @user421473 The line you quoted is not "the update code", just the initial setup of 60 running points. The update code comes later in the program. What I did was just a slight modification of the standard (as it turns out) approach slowly reducing the search area. I guess I'll have to talk a bit about what the general idea of random search is and when it is efficient and when it isn't, but that has to wait a bit :-). $\endgroup$ – fedja Jan 22 '20 at 4:17
  • $\begingroup$ @user421473 I added some common knowledge stuff to the answer. Meanwhile I'm trying to see if I can handle the more general model when the maxima are like in example 2. That requires learning the local shape in the process, which is close to what you wanted to do though my approach to it is somewhat different. $\endgroup$ – fedja Jan 23 '20 at 4:48
  • $\begingroup$ I played around with scipy that @ACheca suggested, and it seems like the shgo algorithm is the most efficient one but only for low dimensions, otherwise the differential evolution that you or him suggested, however they are still slower than the optimized optimizer I have currently, so I doubt your solution will be better (it can't outperform differential evolution), so I stick to what I have. Thanks for all the help anyway! $\endgroup$ – user421473 Feb 2 '20 at 18:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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