%% Nonnegative Linear Least Squares, Solver-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)

% The C matrix has 2000 rows and 400 columns. Therefore, to have the correct size for matrix multiplication, the x vector has 400 rows. To represent the nonnegativity constraint, set lower bounds of zero on all variables.

lb = zeros(size(C,2),1);

% Solve the problem using lsqlin.

[x,resnorm,residual,exitflag,output] = ...
    lsqlin(C,d,[],[],[],[],lb);

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

disp(output)

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

options = optimoptions('lsqlin','Algorithm','trust-region-reflective');
[x2,resnorm2,residual2,exitflag2,output2] = ...
    lsqlin(C,d,[],[],[],[],lb,[],[],options);
disp(output2)

% This time, the solver takes more iterations and reaches a solution with a higher (worse) first-order optimality measure.
% To improve the first-order optimality measure, try setting the SubproblemAlgorithm option to 'factorization'.

options.SubproblemAlgorithm = 'factorization';
[x3,resnorm3,residual3,exitflag3,output3] = ...
    lsqlin(C,d,[],[],[],[],lb,[],[],options);
disp(output3)

% Using this option brings the first-order optimality measure nearly to zero, which is the best possible result.
%
%% Change Solver
% Try solving t problem using the lsqnonneg solver, which is designed to handle nonnegative linear least squares.

[x4,resnorm4,residual4,exitflag4,output4] = lsqnonneg(C,d);
disp(output4)

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

t = table(resnorm - 22.5794, resnorm2 - 22.5794, resnorm3 - 22.5794, resnorm4 - 22.5794,...
    'VariableNames',{'default','trust-region-reflective','factorization','lsqnonneg'})

% The default lsqlin algorithm has a higher residual norm than the trust-region-reflective algorithm. The factorization and lsqnonneg residual norms are even lower, and are the same at this level of display precision. See which one is lower.

disp(resnorm3 - resnorm4)

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