The discussion revolves around the implementation of a function minimization algorithm in MATLAB and C, specifically using the Nelder-Mead method. Participants are comparing the outputs of their respective codes and attempting to identify discrepancies in results when minimizing a function.
Participants generally agree that the C implementation is not producing the same results as MATLAB, but there is no consensus on the cause of the discrepancy. Multiple competing views regarding the correctness of the implementations and the nature of the solutions remain unresolved.
Participants mention potential issues with parameter passing, boundary conditions, and the nature of the minimization algorithm, but these aspects remain unresolved and are subject to further investigation.
Readers interested in numerical optimization, function minimization algorithms, and comparisons between MATLAB and C implementations may find this discussion relevant.
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 % 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;
endmathworks said:
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]));
}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;
}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.
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;double p[5], // 5 is max(2, 3, 5)
x[16]; // 16 is max(2, 3, 5, 6, 16) for(i=0; i<LM_INFO_SZ; ++i)
printf("%g ", info[i]);
printf("\n");