Minimizing function in Matlab and C

[I like Serena]

If you upload neldermead.c, I'll take a look at it.

 

[a.mlw.walker]

neldermead.c, i have had to attach it as a textfile

 

[I like Serena]

Well, here are my first observations.

There does not appear to be a method you can select.

However, one of the parameters of the function NelderMeadSimplexMethod(), which I assume is the one you are using, is "eps". Another one is "timeout",
Apparently you have passed eps=0.1 to it.
You might try using e.g. eps=0.01.
"timeout" is the maximum number of iterations.

 

[a.mlw.walker]

my eps is:
double eps = 1.0e-12;//accuracy of Nelder Mead

?

 

[I like Serena]

All right, that would mean that you have indeed a local minimum that won't approach closer than an sse=0.99.
I find the number a little weird since it is so close to a nice round number.

What did you use for the "timeout" and for the initial guess "xinit"?

 

[a.mlw.walker]

initial guess is the same as matlan, 1, 1, -1.
timeout is
int timeout = 1000000;

it is strange, i asked some guys on a c forum, they have told me i need to get rid of my global variables before i can know whether it is me or the algorithm. I think I am at a local minimum

 

[I like Serena]

Yes, I also think you're a another local minimum, and as yet I can't explain why the C version would find another minimum than the MATLAB version.

It appears as if MATLAB is taking smaller steps, so perhaps they are using different (smaller?) parameters (al, bt, gm) to Nelder-Mead as you already tried.

Did you already print intermediate values of a,b,c?
For the C version as well as the MATLAB version?
This would show which stepsize has been taken, and also whether the same algorithm is used.

Btw, it won't matter (in this case) that you have global variables.

 

[a.mlw.walker]

Yeah that's what I reckon about the globals, they are not the solution to the problem. I'll get the full first 3 iterations and see what's doing what. I read about and nelder-mead is apparently known to fail.
EDIT:
The below values are the MATLAB values for the first three iterations - however sse is the total sum of the error, where as in the C (image) it shows the steps adding up sum, so the last value of sum in each iteration, should match sse in matlab.

a =
     1
b =
     1
c =
    -1
Fitted_Curve =
    6.1378   12.4210   18.7041   24.9873   31.2705
Error_Vector =
   -4.9378   -9.7280  -14.3791  -18.8563  -23.1455
sse =
  1.2171e+003
 
 Iteration   Func-count     min f(x)         Procedure
     0            1          1217.05         
a =
    1.0500
b =
     1
c =
    -1
Fitted_Curve =
    6.1447   12.4279   18.7111   24.9943   31.2774
Error_Vector =
   -4.9447   -9.7349  -14.3861  -18.8633  -23.1524
sse =
  1.2180e+003
a =
     1
b =
    1.0500
c =
    -1
Fitted_Curve =
    5.8455   11.8295   17.8135   23.7975   29.7814
Error_Vector =
   -4.6455   -9.1365  -13.4885  -17.6665  -21.6564
sse =
  1.0681e+003

Interesting, I just noticed Mtlab has a fourth parameter deciding 'next move' called shrink, where the C only has expand contract and reflect, MATLAB has sigma for shrinking, and each expand/contract/reflect then 'shrinks' it:

 % Perform an inside contraction
                xcc = (1-psi)*xbar + psi*v(:,end);
                x(:) = xcc; fxcc = funfcn(x,varargin{:});
                func_evals = func_evals+1;
                
                if fxcc < fv(:,end)
                    v(:,end) = xcc;
                    fv(:,end) = fxcc;
                    how = 'contract inside';
                else
                    % perform a shrink
                    how = 'shrink';
                end
            end
            if strcmp(how,'shrink')
                for j=two2np1
                    v(:,j)=v(:,1)+sigma*(v(:,j) - v(:,1));
                    x(:) = v(:,j); fv(:,j) = funfcn(x,varargin{:});
                end
                func_evals = func_evals + n;
            end

 

[a.mlw.walker]

The attached is the latest information. Top left is MATLAB producing the required minimum
a = 0.025
b = 1.899
c = -0.369
top right is the c code when the parameters
al = 1, bt = 0.5, gm = 2;
a = 3.52
b = 3.63
c = -0.000073
bottom left is when parameters
al = 2, bt = 0.5, gm = 1;
a = 0.076
b = 3.45
c = -0.672
Very odd.

 

[I like Serena]

As you can see MATLAB uses an entirely different stepsize and also takes entirely different steps.

I think MATLAB does not use Nelder-Mead, but something like the Conjugate gradient method.
As you can see it takes a small step in each of the parameters a,b,c.
This is the method to deduce the gradient.
Then the gradient is used to make a step.

I see that the MATLAB fminsearch() function returns a set of values [x,fval,exitflag,output].
The 'output' value contains the algorithm being used.
Can you print these return values?

 

[a.mlw.walker]

No probs, the ouput is:
Estimates =
0.0250 1.8944 -0.3690
FVAL =
0.0037
EXITFLAG =
1
OUTPUT =
iterations: 85
funcCount: 151
algorithm: 'Nelder-Mead simplex direct search'
message: [1x196 char]

 

[I like Serena]

That looks like MATLAB does use Nelder-Mead. :)

So the question becomes: what should the first iteration of Nelder-Mead do?
You can find the algorithm with Google.

I'm not sure yet what's going on with the shrink step.
It may be somewhat "hidden" in the C code, but I haven't looked yet.
And there should be a 4th parameter for the shrink step.

 

[a.mlw.walker]

I think the scalar could be this line in the C
scalarvector(nvar, simp, 0.50, simp);
the 0.5 is the same as the shrink value in the matlab, but it doesn't solve the prob...

 

[I like Serena]

From: http://www.mathworks.com/help/techdoc/math/bsotu2d.html#bsgpq6p-11

The algorithm first makes a simplex around the initial guess x0 by adding 5% of each component x0(i) to x0, and using these n vectors as elements of the simplex in addition to x0.

This is what MATLAB does, but not what the C code does.
Looking at the C code, you'll see that the function initial_simplex() has a different way to initialize the initial simplex.
This means that MATLAB and the C code both implement different variants of the Nelder-Mead algorithm.
This is not wrong, but will lead to different results.

If you want more optimal results (more robust against degeneracy) you should for instance use Levenberg-Marquardt.
The catch is that this algorithm requires you to supply an extra function that calculates the derivatives of your curve fitting function with respect to each parameter.

 

[a.mlw.walker]

Yeah, I thought that. I don't want to use leveneberg marquardt tho for that reason. I was wondering about a ant colony search - that seems pretty robust, but I don't think my roblem should require that. Do you think it is just the initial simplex that is different, and I could adjust the initial simplex code to match the 5% thingy?

 

[I like Serena]

Ah, well, it's not that difficult to calculate the derivatives for your curve-fitting function and you'll get way better results!

But yes, I suspect you can adjust the initial simplex to match the 5% thingy.
Beyond the initial simplex, the algorithm appears to be deterministic.

 

[a.mlw.walker]

Levenberg Marquardt. Do you know of any code? It seems harder to implement than the downhill simplex

 

[I like Serena]

I have the "Numerical Recipes in C" code (copyrighted I'm afraid though not expensive) that implements it and I have used that in the past.
I was very happy with it and swore to myself that would be the only algorithm I'd ever use to fit models (assuming I'd be able to calculate the derivatives).

It's just a function with a bunch of parameters just like your Nelder-Mead function.
The only real difference is that you need to supply a function as an argument that also calculates the derivatives for the set of parameters.

A quick Google search yielded this link:
http://www.ics.forth.gr/~lourakis/levmar/
which appears to be a freely downloadable implementation.

 

[a.mlw.walker]

levmar hey. if i get it, would you be able to help me implement it if i struggle?

 

[I like Serena]

Sure! :)

 

[a.mlw.walker]

Cool! Ok, I'll get it, take a look, and see how far I can get.

 

[a.mlw.walker]

So I got levmar running on the rosenbrock function. Quite neat actually. Interesting that if you give it the jacobian it solves the problem with 10 times fewer iterations. Side note.

I can't find the function that is actually solving the problem - called dlevmar_dif. I want to add my arrays into the function calls so that I can pass the function information. the rosenbrock function doesn't need to be passed anything. anyway just to let you know where i;m at.

 

[I like Serena]

Nice!
I take it you're getting a good solution?
Do you get the same solution or yet another, better or worse?

 

[a.mlw.walker]

Still looking for that function - can only find the pointer at the moment. I am going to do it globally because i can't find it - and just get an asnwer. I'll get back to you

 

[a.mlw.walker]

Hey, So I am trying to get my function now to match the style of the function being solved in levmar.c
The supplied function as the example is:

void ros(double *p, double *x, int m, int n, void *data)
{
register int i;

  for(i=0; i<n; ++i)
    x[i]=((1.0-p[0])*(1.0-p[0]) + ROSD*(p[1]-p[0]*p[0])*(p[1]-p[0]*p[0]));
}

My function is:

static double function(double *p, double *sum, int m, int n, void *data)
{

    double Fitted_Curve[5] = {0};
    double Error_Vector[5] = {0};
    int i;
    double a, b, c;
    sum = 0;
//    double  v1[5], v2[5], geom_inv[5], norm = 0, norm_2 = 0, v2sum = 0, x_coeff = 0;
//    Actual_Output[5] = {1.2, 2.693, 4.325, 6.131, 8.125};
   // printf("Actual_Output[0[ = %f", Actual_Output[0]);
a = p[0];//parameters
b=p[1];
c=p[2];

  for (i = 0; i <= 4; i++)
  {

    Fitted_Curve[i] = (1/(a*b))*(c-asinh(sinh(c)*exp(a*Input[i]*2*pi)));
    Error_Vector[i] = Actual_Output[i]-Fitted_Curve[i];
    printf("a(x[0]) = %f, b(x[1]) = %f, c(x[2]) = %f\n",p[0], p[1], p[2]);
    printf("Fitted_Curve= %f\n", Fitted_Curve[i]);
    printf("Error_Vector= %f\n", Error_Vector[i]);
    }

    for (i = 0; i <= 4; i++)
    {
        sum = sum + pow(Error_Vector[i],2);
        printf("sum= %f\n", sum);
    }

    a_global = a;
    b_global = b;
//    x_coeff_global = x_coeff;
    return sum;

}

p is the vector of my parameters, n and m at the moment I can't tell why there are two of them because I can't find the function that seems to solve the problem - although it runs, but between n and m, I think they state the number of parameters being solved.
I can't find any definition of data, but my problem is that the demo code doesn't have a return - it returns the array x[]. However I am trying to minimize the sum of the errors squared, so I want to return sum (I think), however I am not sure how to adjust the function so that I can do this. Did that make sense?
I have attached the output - It seems to only somplete one iteration, but I am not sure it knows I want to minimize sum?

EDIT: I found this on the website:

int dlevmar_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata), /* functional relation describing measurements.
                                                                 * A p \in R^m yields a \hat{x} \in  R^n
                                                                 */
double *p,         /* I/O: initial parameter estimates. On output contains the estimated solution */
double *x,         /* I: measurement vector. NULL implies a zero vector */
int m,             /* I: parameter vector dimension (i.e. #unknowns) */
int n,             /* I: measurement vector dimension */
int itmax,         /* I: maximum number of iterations */
double opts[5],    /* I: opts[0-4] = minim. options [\tau, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively the
                    * scale factor for initial \mu, stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2 and the
                    * step used in difference approximation to the Jacobian. If \delta<0, the Jacobian is approximated
                    * with central differences which are more accurate (but slower!) compared to the forward differences
                    * employed by default. Set to NULL for defaults to be used.

 

[I like Serena]

How did you call the dlevmar_dif() function? Which parameters?

In your case m=3 (3 parameters).
And n=1 (your result vector only contains your sse/sum).

For now, I'd call dlevmar_dif() with opts=NULL.

The void* adata parameter looks like the common way to pass custom-data to your function. This should also be a parameter to dlevmar_dif() that does not do anything except pass the parameter to your function.
I guess you don't need this (yet).

What is this "minimization info" that you're displaying?
I can tell you that "-1#.IND" represents "not-a-number", meaning that an calculation was done that is invalid (such as taking the power of a negative number).

And what does the "reason 6" mean?

 

[I like Serena]

I just realized! You've changed the code of your function!

You replaced:
sum = sum + Error_Vector*Error_Vector;
by
sum = sum + pow(Error_Vector, 2);

That would cause the error you are seeing, since pow() is not defined for negative numbers, and your numbers do get negative!

 

[a.mlw.walker]

Hi thanks for the replies.
I changed the equation for sum however that didnt get rid of the #IND. It is only doing one iteration from the console screen.
The code where the function is called in main is

  m=3; n=5;
    p[0]=1.0; p[1]=1.0; p[2]=-1.0;
    for(i=0; i<n; i++) x[i]=0.0;
    //ret=dlevmar_der(ros, jacros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic Jacobian
    ret=dlevmar_dif(function, p, x, m, n, 10000, NULL, info, NULL, NULL, NULL);  // no Jacobian
  break;

However at the beginning of main, is this:

double p[5], // 5 is max(2, 3, 5)
	   x[16]; // 16 is max(2, 3, 5, 6, 16)

x is what the originial example was passing back - second parameter in the function ros.
I am just passing a variable back, is that ok?

From changing n to 1, I get a message out saying that the minimum number of measurement can't be less than the number of unknowns, so i changed it to 5 - my number of measurements.
The minimization info comes from:

  for(i=0; i<LM_INFO_SZ; ++i)
    printf("%g ", info[i]);
  printf("\n");

I just don't know what info holds, and where it is set at the moment

 

[I like Serena]

I just looked the web page over.
According to my interpretation you should:

At the beginning of main, you should have:
double p[3] = {0,0,0};
double x[5] = {0,0,0,0,0};
int n = 5;

In your function, set x to Error_Vector for each i (i=0..4).
That is, you should not calculate the sse/sum yourself.

The contents of 'info' are documented on the web page, but let's first call the function correctly, before interpreting it.

 

[a.mlw.walker]

Cool thanks for that. It fixed all the so far known issues. Check image - says completed no iterations - due to reason 7 (I need to find that). I have set m to 3 as that is what it was before.

Edit: found reason 7 on website:
7 - stopped by invalid (i.e. NaN or Inf) "func" values; a user error

The minimizing error:

O: information regarding the minimization. Set to NULL if don't care
* info[0]= ||e||_2 at initial p.
* info[1-4]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, \mu/max[J^T J]_ii ], all computed at estimated p.
* info[5]= # iterations,
* info[6]=reason for terminating: 1 - stopped by small gradient J^T e
* 2 - stopped by small Dp
* 3 - stopped by itmax
* 4 - singular matrix. Restart from current p with increased \mu
* 5 - no further error reduction is possible. Restart with increased mu
* 6 - stopped by small ||e||_2
* 7 - stopped by invalid (i.e. NaN or Inf) "func" values; a user error
* info[7]= # function evaluations
* info[8]= # Jacobian evaluations
* info[9]= # linear systems solved, i.e. # attempts for reducing error