%% Nonlinear Least Squares Without and Including Jacobian
% This example shows how to solve a nonlinear least-squares problem in two ways. The example first solves the problem without using a Jacobian function. Then it shows how to include a Jacobian, and illustrates the resulting improved efficiency.
% The problem has 10 terms with two unknowns: find x, a two-dimensional vector, that minimizes
% \[\sum\limits_{k = 1}^{10} {{{\left( {2 + 2k - {e^{k{x_1}}} - {e^{k{x_2}}}} \right)}^2}} ,\]
% starting at the point x0 = [0.3,0.4].
% Because lsqnonlin assumes that the sum of squares is not explicitly formed in the user function, the function passed to lsqnonlin must compute the vector-valued function
% \[{F_k}(x) = 2 + 2k - {e^{k{x_1}}} - {e^{k{x_2}}},\]
% for k = 1 to 10 (that is, F must have 10 components).
%% Solve Problem Without Jacobian
% The helper function myfun defined at the end of this example implements the vector-valued objective function with no derivative information. Solve the minimization starting from the point x0.

x0 = [0.3,0.4]; % Starting guess
[x,resnorm,res,eflag,output] = lsqnonlin(@myfun,x0); % Invoke optimizer

% Examine the solution and number of function evaluations.

disp(x)
disp(resnorm)
disp(output.funcCount)

%
%% Solve Problem Including Jacobian
% The objective function is simple enough that you can calculate its Jacobian. Following the definition in Jacobians of Vector Functions, a Jacobian function represents the matrix
% \[{J_{kj}}(x) = \frac{{\partial {F_k}(x)}}{{\partial {x_j}}}.\]
% Here, F_k(x) is the kth component of the objective function. This example has
% \[{F_k}(x) = 2 + 2k - {e^{k{x_1}}} - {e^{k{x_2}}},\]
% so
% \[\begin{array}{l}
% {J_{k1}}(x) =  - k{e^{k{x_1}}}\\
% {J_{k2}}(x) =  - k{e^{k{x_2}}}.
% \end{array}\]
% The helper function myfun2 defined at the end of this example implements the objective function with the Jacobian. Set options so the solver uses the Jacobian.

opts = optimoptions(@lsqnonlin,'SpecifyObjectiveGradient',true);

% Run the solver.

lb = []; % No bounds
ub = [];
[x2,resnorm2,res2,eflag2,output2] = lsqnonlin(@myfun2,x0,lb,ub,opts);

% The solution is the same as the previous solution.

disp(x2)
disp(resnorm2)

% The advantage of using a Jacobian is that the solver takes many fewer function evaluations.

disp(output2.funcCount)

%
%% Helper Functions
% This code creates the myfun helper function.

function F = myfun(x)
k = 1:10;
F = 2 + 2*k-exp(k*x(1))-exp(k*x(2));
end

% This code creates the myfun2 helper function.

function [F,J] = myfun2(x)
k = 1:10;
F = 2 + 2*k-exp(k*x(1))-exp(k*x(2));
if nargout > 1
    J = zeros(10,2);
    J(k,1) = -k.*exp(k*x(1));
    J(k,2) = -k.*exp(k*x(2));
end
end

% Copyright 2019 The MathWorks, Inc.