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,D ata) 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 doesnt 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! :smile:

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 doesnt find a curve that fits. input is just [1 2 3 4 5], thats 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 isnt 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. Im 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 cant 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 havent put them because its quite alot of code. right now i am trying to get rid of all my global variables so it doesnt compile right now, because i have an error im 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[i], 0.50, simp[i]);
the 0.5 is the same as the shrink value in the matlab, but it doesnt 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 dont want to use leveneberg marquardt tho for that reason. I was wondering about a ant colony search - that seems pretty robust, but I dont 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 cant 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 doesnt need to be passed anything. anyway just to let you know where i;m at.

[I like Serena]

Nice! :smile:
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. Im gonna do it globally because i cant 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 cant tell why there are two of them because I cant 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 cant find any definition of data, but my problem is that the demo code doesnt 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! :surprised

You replaced:
sum = sum + Error_Vector[i]*Error_Vector[i];
by
sum = sum + pow(Error_Vector[i], 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 begining 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 cant 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 dont 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[i] to Error_Vector[i] 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

[I like Serena]

According to the web page:

reason 7 = stopped by invalid (i.e. NaN or Inf) "func" values; a user error
return value -1 = no iterations performed.

As you can see you have a Fitted_Curve of 1#.QNAN0, meaning your calculation failed.
You're somehow feeding your function data that it can't calculate.
I see no reason for that, so I think you made a programming error.

Can you show your function as it is now, and also the values you're feeding it, and the Fitted_Curve/Error_Vector results, when it gives the 1#.QNAN0?

Note that the only possible reason that you might not be able to calculate the result, is if you have either parameter a or b given as zero.

[a.mlw.walker]

definitions out side of functions:

#define pi 3.141592653589793238462643

double Input[5] = {0, 1, 2, 3, 4};
double Actual_Output[5] = {1.2, 2.693, 4.325, 6.131, 8.125};

Function

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

double Fitted_Curve[5] = {0};
// double Error_Vector[5] = {0};
int i;
double a, b, c;

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)));
x[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", x[i]);
}
}

main:

int main()
{
register int i;
int problem, ret;
double p[3] = {0,0,0};
double x[5] = {0,0,0,0,0};
int n = 5;
int m;
double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
char *probname[]={
"function"
};

opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
//opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute!

problem=
0; // Can have more functions if user requires

#ifdef HAVE_LAPACK

#else // no LAPACK
#ifdef _MSC_VER
#pragma message("LAPACK not available, some test problems cannot be used")
#else
#warning LAPACK not available, some test problems cannot be used
#endif // _MSC_VER

#endif /* HAVE_LAPACK */

switch(problem){
default: fprintf(stderr, "unknown problem specified (#%d)! Note that some minimization problems require LAPACK.\n", problem);
exit(1);
break;

case 0:
/* function */
m=3;
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, 100000, NULL, NULL, NULL, NULL, NULL); // no Jacobian
break;


} /* switch */

printf("Results for %s:\n", probname[problem]);
printf("Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);
for(i=0; i<m; ++i)
printf("%.7g ", p[i]);
printf("\n\nMinimization info:\n");
for(i=0; i<LM_INFO_SZ; ++i)
printf("%g ", info[i]);
printf("\n");

return 0;
}

[I like Serena]

For some reason I'm not able to look at your attached screenshot.
Can you view it yourself?
Can you perhaps upload it again?

First things that attract my attention, is that the return value of your function is 'double', while it should be 'void'.
Then, I recommend initializing the info-array to zeroes.
Furthermore, you should pass the info-array as a parameter, if you want to inspect the output results.

Finally, your compiler can warn you of mistakes like these.
You should increase the warning-level of your compiler and fix mismatches such as the prototype of your function.
Which compiler do you use?

[a.mlw.walker]

I just changed the prototype to void, and initialised infor to zero's by putting it equal to {0}.
It runs much more succesfully now - 174 iterations, but stops due to reason 2 - stopped by small Dp.
The solution isnt correct though.
I have two warnings from GCC Compiler now, but both are generated due to warnings within the code such as
#warning LAPACK not available, some test problems cannot be used.
p[0] should = 0.025
p[1] should = 1.8944
p[2] should = -0.369
The solution it has found I think is the solution from when I was using neldermead. Now im questioning the function

[I like Serena]

I believe you used a different Input[] array before.
Shouldn't you use Input[] = {1,2,3,4,5}?
I believe an Input[0] value of 0, yields a Fitted_Curve of 0, but that was never in your output before.

Btw, it has to be a different solution than the one from Nelder-Mead.
It reports an SSE=0.0015893 which is even more accurate than what you got from Matlab!

What do you get with the Input[] array of {1,2,3,4,5}?

And if you want even more accurate results, you should use dlevmar_der() with the Jacobian specified.
Do you know how to calculate the Jacobian?

[a.mlw.walker]

Hey, sorry I did use the correct input, i just showed you the wrong one.
So now I am confused, the matlab produces a different result, but not a hugely different result. Those values produce a really good fit (matlab plotter). But the original author did not get those numbers.
Now to get rid of my globals...

[I like Serena]

Well, you can expect an even better fit from your new solution. ;)
You should try it.

[a.mlw.walker]

I have the green line is the new fit, the red is the matlab. I trust our new levmar one now more!
So this *data pointer can pass the variables required by the function. Do you know how its used?
Also I have another minimization I would like to use levmar for, however it is the sum of the modulus of the function as supposed to the sum of the squares, I know that if you square and squareroot a function it is like the modulus, but will it make any difference?

[I like Serena]

You seem to have left out the picture...

To minimize the sum of absolute values, I guess you can make it work by returning the square root of the absolute value of each Error_Vector[i].


The *data pointer should be used to pass your Actual_Data[5] array.
It's like this:


void function(..., void *adata)
{
double *Actual_Data[5] = (double *)adata;
...
}

main()
{
double Actual_Data[5] = {...};
int ret = dlevmar_dif(..., Actual_Data);
}

[a.mlw.walker]

forgot graph...

[I like Serena]

Yep, the green line fits better! :smile:

[a.mlw.walker]

Im getting an error on this line:
double *Actual_Data[5] = (double *)adata;
Invalid Initializer.


This function (my other function) runs but produces very different values from the matlab

void function2(double *p, double *x, int m, int n, void *data)//equations for parameters eta, phi and omega^2 that hold wheel properties
{
double c1[23], x_coeff_2[23], v1[23], v2[23];
double eta = p[0]; //parameters
double phi=p[1]; //parameters
double omega_2=p[2]; //parameters

int i;
a_global = 0.025;
b_global = 1.8944;
//c_global = -0.3690;
for (i = 0; i<23; i++)
{
v1[i] = (exp(a_global*2*pi)-cosh(a_global*b_global*init_T[i][0]));
v2[i] = sinh(a_global*b_global*init_T[i][0]);
x_coeff_2[i] = v1[i]/v2[i];
c1[i]= b_global*b_global*(x_coeff_2[i]*x_coeff_2[i]-1);
x[i] = sqrt(fabs(c1[i]*exp(-2*a_global*init_T[i][1])+eta*((1+0.5*(4*a_global*a_global+1))*cos(init_T[i][1]+phi)-2*a_global*sin(init_T[i][1]+phi))+b_global*b_global-omega_2));
}

}

I had to adjust it from the nelder mead (where it ran) because of the way levmar wants x[] passed back to it however the nelder mead matched the matlab. This issue is reason 2 (small Dp) again

[I like Serena]

Im getting an error on this line:
double *Actual_Data[5] = (double *)adata;
Invalid Initializer.

Sorry, I intended:
double *Actual_Data = (double *)adata;


This function (my other function) runs but produces very different values from the matlab

I had to adjust it from the nelder mead (where it ran) because of the way levmar wants x[] passed back to it however the nelder mead matched the matlab. This issue is reason 2 (small Dp) again

Well, is it a good solution that it found?
What is info[1], which is the found SSE?

To be honest, I'm not sure Levenberg-Marquardt will work as intended with your function, since it is not intended for a sum of absolute errors, but for a sum of squared errors.
I believe it should still work, but it may not find an optimal solution.
Nelder-Mead may be a better choice for this problem (minimizing the sum of absolute errors).

[a.mlw.walker]

Hi so I have integrated the levmar and neldermead into one program so that I can use levmar to solve function1 and neldermead to solve function 2.What about prototypes tho? When do you need them and when do you not. For instance my functions function() and function2() dont have prototypes, however some other functions do have prototypes.

My third function:

static double rotor_function(int n, double z[],float ti[6])//equations for parameters eta, phi and omega^2 that hold wheel properties
{

double kappa = z[0], sum = 0,Rotor_Curve[6], Error_Vector[6], v0; //parameters
int i;
v0 = (1 + (kappa/2)*ti[0]*ti[0])/ti[0];
for (i = 0; i<6; i++)
{
Rotor_Curve[i] = (v0*ti[i]-(kappa/2)*pow(ti[i],2));
Error_Vector[i] = (i+1) - Rotor_Curve[i];
}
for (i = 0; i<6; i++)
{
sum = sum + Error_Vector[i]*Error_Vector[i];
}
printf("\n z = %f v0 = %f, sum = %f\n", z[0], v0, sum);
return sum;

}
doesnt have a prototype - however it also doesnt receive the initital guess at the parameter.
It is called like:
case 2:
/* function 2*/
z[0] = 0.004;

n=1;
fopt = 0;
//This is still not running correctly, check matlab for correct parameter values:
if (NelderMeadSimplexMethod(n, rotor_function, z, length, &fopt, timeout, eps, ti) == success) {
printf("reaching to minimum ");
} else {
printf("timeout ");
}
kappa = z[0];

i need to pass the parameter as an array though I think not a variable

[I like Serena]

Hi so I have integrated the levmar and neldermead into one program so that I can use levmar to solve function1 and neldermead to solve function 2.What about prototypes tho? When do you need them and when do you not. For instance my functions function() and function2() dont have prototypes, however some other functions do have prototypes.

As long as you define your function (including the body) before using it, you don't need a function prototype.
If you need a prototype, you'll know, because the compiler will complain he doesn't know the function.

I need to pass the parameter as an array though I think not a variable

I can't recall the list of parameters to your Nelder-Mead function.
But if it doesn't have a "void *data" pointer as a parameter, you need to use a different mechanism, which is as follows.

The way to do it, is to define a global variable to pass your Actual_Data[] array.
You use this global variable in your function.
And just before you call Nelder-Mead, you copy your Actual_Data[] array into the global variable.

[a.mlw.walker]

Neldermead function is like this:

status NelderMeadSimplexMethod(n, f, xinit, length, fopt, timeout, eps, array1)

my function that is failing is like this:

static double rotor_function(int n, double z[],float ti[6])//equations for parameters eta, phi and omega^2 that hold wheel properties
{

double kappa = z[0], sum = 0,Rotor_Curve[6], Error_Vector[6], v0; //parameters
int i;
v0 = (1 + (kappa/2)*ti[0]*ti[0])/ti[0];
for (i = 0; i<6; i++)
{
Rotor_Curve[i] = (v0*ti[i]-(kappa/2)*pow(ti[i],2));
Error_Vector[i] = (i+1) - Rotor_Curve[i];
}
for (i = 0; i<6; i++)
{
sum = sum + Error_Vector[i]*Error_Vector[i];
}
printf("\n z = %f v0 = %f, sum = %f\n", z[0], v0, sum);
return sum;

}

I am calling both like this:

z[0] = 0.004;

n=1;
fopt = 0;
//This is still not running correctly, check matlab for correct parameter values:
if (NelderMeadSimplexMethod(n, rotor_function, z, length, &fopt, timeout, eps, ti) == success) {
printf("reaching to minimum ");
} else {
printf("timeout ");
}
kappa = z[0];

It seems to me like function2 however this one doesnt ever have a value for z[0]

[a.mlw.walker]

I found the problem. z is getting there howere the maths is going haywire! v0 is not being calculated correctly - is this because of a problem with multiplying floats and doubles?

[I like Serena]

It's weird to see a parameter like "float ti[6]".
I suspect it should be "double ti[6]".
The mismatch could cause your problems.

Don't you get a compiler warning?
What is the prototype of your Nelder-Mead function, and what function prototype does it expect exactly as a parameter?

[a.mlw.walker]

Hey how are you, been a while.
I am working more with the levenberg marquardt code now.
I am trying to get it to use just single precision - this is an option in levmar.h
Out of interest, what happens if a variable has been declared as double but in levmar.h i state i want single precision?
when I choose single precision I get an error - a strange one at that.
Axb_core.c|1103|error: 'FLT_EPSILON' undeclared (first use in this function)

It works fine with double but not with single. I find this odd considering the author offeres single precision accuracy.
Any ideas what the problem could be?