%% Nonnegative Linear Least Squares, Problem-Based
% This example shows how to use several algorithms to solve a linear least squares problem with the bound constraint that the solution is nonnegative. A linear least squares problem has the form
% \min_x \|Cx - d\|^2.
% In this case, constrain the solution to be nonnegative, x \ge 0.
% To begin, load the arrays C and d into your workspace.

load particle

% View the size of each array.

sizec = size(C)
sized = size(d)

% Create an optimization variable x of the appropriate size for multiplication by C. Impose a lower bound of 0 on the elements of x.

x = optimvar('x',sizec(2),'LowerBound',0);

% Create the objective function expression.

residual = C*x - d;
obj = sum(residual.^2);

% Create an optimization problem called nonneglsq and include the objective function in the problem.

nonneglsq = optimproblem('Objective',obj);

% Find the default solver for the problem.

opts = optimoptions(nonneglsq)

% Solve the problem using the default solver.

[sol,fval,exitflag,output] = solve(nonneglsq);

% To see details of the optimization process, examine the output structure.

disp(output)

% The output structure shows that the lsqlin solver uses a sparse internal linear solver for the interior-point algorithm and takes 9 iterations to arrive at a first-order optimality measure of about 1e-6.
%
%% Change Algorithm
% The trust-region-reflective algorithm handles bound-constrained problems. See how well it performs on this problem.

opts.Algorithm = 'trust-region-reflective';
[sol2,fval2,exitflag2,output2] = solve(nonneglsq,'Options',opts);
disp(output2)

% This time, the solver takes more iterations and arrives at a solution with a lower (better) first-order optimality measure.
% To get an even better first-order optimality measure, try setting the SubproblemAlgorithm option to 'factorization'.

opts.SubproblemAlgorithm = 'factorization';
[sol3,fval3,exitflag3,output3] = solve(nonneglsq,'Options',opts);
disp(output3)

% Using this option brings the first-order optimality measure nearly to zero, which is the best possible.
%
%% Change Solver
% The lsqnonneg solver is specially designed to handle nonnegative linear least squares. Try that solver.

[sol4,fval4,exitflag4,output4] = solve(nonneglsq,'Solver','lsqnonneg');
disp(output4)

% lsqnonneg does not report a first-order optimality measure. Instead, investigate the residual norms, which are returned in the fval outputs. To see the lower-significance digits, subtract 22.5794 from each residual norm.

t = table(fval - 22.5794, fval2 - 22.5794, fval3 - 22.5794, fval4 - 22.5794,...
    'VariableNames',{'default','trust-region-reflective','factorization','lsqnonneg'})

% The default solver has a slightly higher (worse) residual norm than the other three, whose residual norms are indistinguishable at this level of display precision. To see which is lowest, subtract the lsqnonneg result from the two results.

disp([fval2 - fval4,fval3 - fval4])

% The lsqnonneg residual norm is the smallest by a nearly negligible amount. However, lsqnonneg takes the most iterations to converge.
%
% Copyright 2018-2019 The MathWorks, Inc.