MATLAB Minimizing function in Matlab and C

[a.mlw.walker]

The foloowing code is the MATLAB code to minimize a function using the nelder mead method.

Starting=[1 1 -1]; options=optimset('Display','iter'); Estimates=fminsearch(@myfit,Starting,options,t_d,Data) Calling: function sse=myfit(params,Input,Actual_Output) a=params(1); b=params(2); c=params(3);

Fitted_Curve = (1/(a*b))*(c-asinh(sinh(c)*exp(a*Input*2*pi)));
Error_Vector=Actual_Output-Fitted_Curve ;
sse=sum(Error_Vector.^2);

This is supposed to be the equivelent C code, however it doesn't produce the same results. Can anyone see what I may have done wrong?

static double function(int n, double x[])
{
	double c;
    double Fitted_Curve[5];
    double Error_Vector[5];
    int i;
    double a, b, 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};
a = x[0];
b=x[1];
c=x[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];
      }
      for (i = 0; i <= 4; i++)
      {
          sum = sum + Error_Vector[i]*Error_Vector[i];
      }
      printf("sum = %f\n", sum);
      a_global = a;
      b_global = b;
     // x_coeff_global = x_coeff;
      return sum;
}

Thanks

 

[I like Serena]

Hi a.mlw.walker!

The C function that you're showing is functionally a match for the myFit() function that you specified as MATLAB code.

Do you also have an implementation of the fminsearch() function?
If so, what does it expect as input?

 

[a.mlw.walker]

the fminsearch function is a MATLAB function implementing the NelderMead algorithm. I know the MATLAB is correct because I have the solution. However the C does not produce the same final result but does produce the same first result!
Matlab produces:
a = 0.0250 b= 1.8944 c= -0.3690
C code produces
a = 3.52, b = 3.634 c= -0.000073
Check the pics:

 

[I like Serena]

You pictures seem to basically show that the C function is indeed a match for the MATLAB function.

So either your fminsearch is not doing its work like it does in matlab.
Or, what seems more likely, the parameters are syntactically not passed properly to the function.
Or the boundary conditions are not set properly.

Did you already print the values of a, b, and c during the iterations?

What are the values in your Input[] array?

Did you already check if the given solution by the C code is actually a solution for the problem, since it may have found another minimum (if your boundary conditions are set up differently)?

 

[a.mlw.walker]

I'll get it to print earlier values of a, b and c and post them. Yeah i plotted the points found by C and it doesn't find a curve that fits. input is just [1 2 3 4 5], that's what it is in both MATLAB and C.
Its a strange error considering the values for the first iteration are the same, but I know the Matlab is correct, I have the solution to that from previous work - and it is a solution. The problem is in the C.
I'll get a b and c for the first iteration

 

[I like Serena]

While you're at it, also print the values of 'sum' (the sse) in your function, since we'd like to see if it converges or not.

 

[a.mlw.walker]

in the C code right, not the matlab?

 

[I like Serena]

Both! :)

 

[a.mlw.walker]

Ok here we go. three iterations of both, MATLAB in white, c in black...

 

[I like Serena]

I'm just noticing that in all your pictures the sum/sse does not match the Error_Vector.
The sum/sse should be the sum of the squared entries of the Error_vector.
In all cases (except perhaps the last) it isn't.

In your original 2 pictures the sum and sse appear to be the square of only the first entry in the Error_Vector.

So what did you print for the sum/sse?

In your last 2 pictures, there are 5 sums displayed in the console dump.
They appear to be the squares of the individual Error_Vector entries, but they should still be added.
In the MATLAB output, there is only one sse, that may well be the sum needed.

This means the code you have shown in the opening post apparently does not match the code you have just executed.

 

[a.mlw.walker]

The original images were just to see whether it was doing the sum my calculator was expecting - i.e just the first Error_Vector^2. The latest images (my last post) is the actual output after three iterations from NelderMead. look at the latest images.

 

[a.mlw.walker]

Sorry your right, the C code for sum is not producing the correct values...

 

[I like Serena]

The original images were just to see whether it was doing the sum my calculator was expecting - i.e just the first Error_Vector^2. The latest images (my last post) is the actual output after three iterations from NelderMead. look at the latest images.

Yes, I believe I have.
The displayed sums for 1 iteration need to be added, and that has to be used as the result.
That should match the MATLAB output.

 

[a.mlw.walker]

Hmm, sorry I had removed that in a recent debugging session. I have put it back in, and attached again the first three iterations, but its still not right
I have changed the sum code to

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

 

[I like Serena]

All right.
Now it shows the same result for the first iteration as matlab.

As you can see, it also converges.
But it converges much faster than the MATLAB version.

This suggests that you're using a different algorithm.
Since you've presumably used a function with the same name (fminsearch), do you have the option of giving method-parameters?
What you have in C might for instance be the Conjugate gradient method.

(You did use a function by the name of fminsearch didn't you?)

 

[a.mlw.walker]

The function for the C is called NelderMead.c. Fminsearch is a built in MATLAB function that according to Matlab is the Neldermead method. I'm not bother how fast it converges if the answer isn't correct though...

 

[I like Serena]

There could be more than one solution.
Could you show the final a, b, c, and sse to judge that?

And perhaps you could also show the a, b, and c in each iteration to see how the algorithm progresses?

 

[I like Serena]

I just checked.
Your problem does have multiple solutions.
So perhaps you want the solution closest to your initial guess?
Or doesn't it matter to you?

 

[a.mlw.walker]

really, mulitple solutions? i tried the values i earlier posted and it did not produce a solution - I am looking for a curve to fit data points, and the curve those values gave did not fit the data like the matlab

 

[I like Serena]

Ah, I missed your solutions for a, b, c before.

They are both solutions.
The MATLAB version has an sse = 0.0037.
The C version has an sse = 0.099

This suggests that the C version stopped at a tolerance of 0.1.
You can probably tweak your parameters somewhere.

Btw, you should know that your data is almost on a straight line, and that whatever parameters you feed your formula, it's also almost a straight line.
So I think you shouldn't expect very interesting solutions.

What are you using this for?

 

[a.mlw.walker]

My boss wanted me to research finding parameters that will fit the curve to an equation. I agree the data is basically on a straight line, however the Matlab code gets the solution I am expecting the C goes off in a completely different direction, even though the methods are supposed to be the same. If it is as simple as the c is finding a different route then what can I do about it - the MATLAB is finding the right one?
By the way I have attached two graphs, the left is the atlab the right is the C, MATLAB is a much better ft...

 

[a.mlw.walker]

Guys what do you reckon. My graph above shows that my code may be finding a minimum of sorts, but that is not a good curve. I think it is finding a local minimum as supposed to the global minimum. Interesting that the MATLAB implementation finds the global minimum. The following website states that nelder mead method can fail in simple situations. i wonder if MATLAB is doing "more" than just a nelder mead. i have heard the rosenbrock algorithm can find global minimum, anyone used it? Or an ant search or something? anyone used that?

 

[a.mlw.walker]

Ok i have noticed something. In the file neldermead.c on line 29 it says

static dbl	al = 2, bt = 0.5, gm = 1;

these are parameters that control the adjustments made by the algorithm.
Matlab has a similar line, however the variable names are different. I fiddled with the values (just changing which one was equal to 1,2 and 0.5) and I get the correct answer for the three parameters I have been looking for, however using an alternative function to solve for, and it gives incorrect results again.

 

[a.mlw.walker]

My last post tho - about the parameters in NelderMead, swapping this around got me a solution to the function - the correct one, however for another function I have written it does not solve it correctly once I change these values.
When I change them back (static dbl al = 1, bt = 0.5, gm = 2, the second function works!
I think the method may be finding local minimums and when i swap the parameters it adjusts to a global minimum! Wotz goin On!

 

[I like Serena]

I can think of 2 things.

First your C version uses a threshold of 0.1.
You should be able to find this threshold somewhere and adjust it.
Then you can probably get a better fit.

Second the MATLAB algorithm has an "options" parameter.
Among others it specifies an algorithm to use (I think).
With the right option, you should be able to get a fit that is as close to the initial guess as possible (I think that would be something like Levenberg-Marquardt).

Finally, you just have a formula for a curve that just won't fit very well when the curve is almost a straight line and you data points are almost a straight line.

 

[a.mlw.walker]

did you look at post #21? the image of the MATLAB fits is great, the C one not so good. I am working on it, I'll see if i get anywhere...

 

[I like Serena]

Yes, the C one has a tolerance of 0.1, which is bigger than the MATLAB one...
And apparently the algorithm doesn't try to stay close to the initial gues...

 

[a.mlw.walker]

Have you seen the 0.1 tolerance anywhere, because I can't find it.

 

[I like Serena]

I only saw it in the result (as posted earlier).
It has to be somewhere in your neldermead.c file, or another one that is used by it.
Since I do not have access to this file I can't tell you.

 

[a.mlw.walker]

would you like access to the files? i only haven't put them because its quite a lot of code. right now i am trying to get rid of all my global variables so it doesn't compile right now, because i have an error I am not sure about. but i can give it to you if your knowledge of c is any better than mine?

 

[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! :)