Generate Code for lsqcurvefit or lsqnonlin

Generate Code for `lsqcurvefit` or `lsqnonlin`

This example shows how to generate C code for nonlinear least squares.

Data and Model for Least Squares

In this example, the vector xdata represents 100 data points, and the vector ydata represents the associated measurements. The modeled relationship between xdata and ydata is

${ydata}_{i} = a_{1} + a_{2} \exp (- a_{3} {xdata}_{i}) + ε_{i} .$

Generate the data for the problem.

rng(5489,'twister') % For reproducibility
xdata = -2*log(rand(100,1));
ydata = (ones(100,1) + .1*randn(100,1)) + (3*ones(100,1)+...
    0.5*randn(100,1)).*exp((-(2*ones(100,1)+...
    0.5*randn(100,1))).*xdata);

The code generates xdata from 100 independent samples of an exponential distribution with mean 2. The code generates ydata from its defining equation using a = [1;3;2], perturbed by adding normal deviates with standard deviations [0.1;0.5;0.5].

Solve Generating Code for `lsqcurvefit`

Solver Approach

The goal is to find parameters for the model ${\hat{a}}_{i}$ , i = 1, 2, 3 that best fit the data.

To fit the parameters to the data using lsqcurvefit, you need to define a fitting function. For lsqcurvefit, the fitting function takes a parameter vector a and the data xdata and returns a prediction of the response, which should be equal to ydata with no noise and a perfect model. Define the fitting function predicted as an anonymous function.

predicted = @(a,xdata) a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata);

To fit the model to the data, lsqcurvefit needs an initial estimate a0 of the parameters.

a0 = [2;2;2];

Call lsqcurvefit to find the best-fitting parameters ${\hat{a}}_{i}$ .

[ahat,resnorm,residual,exitflag,output,lambda,jacobian] =...
   lsqcurvefit(predicted,a0,xdata,ydata);

Code Generation Approach

To solve the same problem using code generation, complete the following steps.

Write a function that incorporates all of the previous steps: generate data, create a fitting function, create an initial point, and call lsqcurvefit.

function [x,res] = solvelsqcurve
rng(5489,'twister') % For reproducibility
xdata = -2*log(rand(100,1));
ydata = (ones(100,1) + .1*randn(100,1)) + (3*ones(100,1)+...
    0.5*randn(100,1)).*exp((-(2*ones(100,1)+...
    0.5*randn(100,1))).*xdata);
predicted = @(a,xdata) a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata);
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt','Display','off');
a0 = [2;2;2];
lb = [];
ub = [];
[x,res] = lsqcurvefit(predicted,a0,xdata,ydata,lb,ub,options);
end

Create a configuration for code generation. In this case, use 'mex'.
```
cfg = coder.config('mex');
```
Generate code for the solvelsqcurve function.
```
codegen -config cfg solvelsqcurve
```
Test the generated code by running the generated file, which is named solvelsqcurve_mex.mexw64 or similar.
```
[x,res] = solvelsqcurve_mex
```
```
x =

    1.0169
    3.1444
    2.1596


res =

    7.4101
```

Solve Generating Code for `lsqnonlin`

Solver Approach

The goal is to find parameters for the model ${\hat{a}}_{i}$ , i = 1, 2, 3 that best fit the data.

To fit the parameters to the data using lsqnonlin, you need to define a fitting function. For lsqnonlin, the fitting function takes a parameter vector a, the data xdata, and the data ydata. The fitting function returns the difference between the prediction of a response and the data ydata, which should equal 0 with no noise and a perfect model. Define the fitting function predicted as an anonymous function.

predicted = @(a)(a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata) - ydata)

To fit the model to the data, lsqnonlin needs an initial estimate a0 of the parameters.

a0 = [2;2;2];

Call lsqnonlin to find the best-fitting parameters ${\hat{a}}_{i}$ .

[ahat,resnorm,residual,exitflag,output,lambda,jacobian] =...
   lsqnonlin(predicted,a0);

Code Generation Approach

To solve the same problem using code generation, complete the following steps.

Write a function that incorporates all of the previous steps: generate data, create a fitting function, create an initial point, and call lsqnonlin.

function [x,res] = solvelsqnon
rng(5489,'twister') % For reproducibility
xdata = -2*log(rand(100,1));
ydata = (ones(100,1) + .1*randn(100,1)) + (3*ones(100,1)+...
    0.5*randn(100,1)).*exp((-(2*ones(100,1)+...
    0.5*randn(100,1))).*xdata);
predicted = @(a) (a(1)*ones(100,1)+a(2)*exp(-a(3)*xdata) - ydata);
options = optimoptions('lsqnonlin','Algorithm','levenberg-marquardt','Display','off');
a0 = [2;2;2];
lb = [];
ub = [];
[x,res] = lsqnonlin(predicted,a0,lb,ub,options);
end

Create a configuration for code generation. In this case, use 'mex'.
```
cfg = coder.config('mex');
```
Generate code for the solvelsqnon function.
```
codegen -config cfg solvelsqnon
```
Test the generated code by running the generated file, which is named solvelsqnon_mex.mexw64 or similar.
```
[x,res] = solvelsqnon_mex
```
```
x =

    1.0169
    3.1444
    2.1596


res =

    7.4101
```
The solution is identical to the one generated by solvelsqcurve_mex because the solvers have identical underlying algorithms. So, you can use the solver you find most convenient.