Nonnegative Least-Squares, Problem-Based

This example shows problem-based approaches for solving the problem

C*x = d

in the least-squares sense. Here C and d are given arrays and x is the unknown array, and x is restricted to be nonnegative. In other words,

minimize sum((C*x - d).^2)

subject to x >= 0.

To begin, load arrays C and d into your workspace.

load particle

View the sizes of the arrays.

sizec = size(C)
sizec = 1×2

        2000         400

sized = size(d)
sized = 1×2

        2000           1

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)
opts = 
  lsqlin options:

   Options used by current Algorithm ('interior-point'):
   (Other available algorithms: 'trust-region-reflective')

   Set properties:
     No options set.

   Default properties:
              Algorithm: 'interior-point'
    ConstraintTolerance: 1.0000e-08
                Display: 'final'
           LinearSolver: 'auto'
          MaxIterations: 200
    OptimalityTolerance: 1.0000e-08
          StepTolerance: 1.0000e-12

   Show options not used by current Algorithm ('interior-point')

Solve the problem using the default solver.

[sol,fval,exitflag] = solve(nonneglsq)
Solving problem using lsqlin.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
sol = struct with fields:
    x: [400x1 double]

fval = 22.5795
exitflag = 
    OptimalSolution

The default solver, lsqlin, worked quickly. But there are other solvers that handle problems of this nature. lsqnonneg solves problems of this type, as does quadprog. Try each.

[sol2,fval2,exitflag2] = solve(nonneglsq,'Solver','lsqnonneg');
Solving problem using lsqnonneg.
[sol3,fval3,exitflag3] = solve(nonneglsq,'Solver','quadprog');
Solving problem using quadprog.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

Examine the first ten entries of each solution.

solutions = table(sol.x(1:10),sol2.x(1:10),sol3.x(1:10),'VariableNames',{'lsqlin','lsqnonneg','quadprog'})
solutions=10×3 table
      lsqlin      lsqnonneg     quadprog 
    __________    _________    __________

      0.049416    0.049416       0.049416
      0.082985    0.082985       0.082985
      0.078469    0.078469       0.078469
       0.11682     0.11682        0.11682
       0.07165     0.07165        0.07165
      0.060525    0.060525       0.060525
    1.2196e-09           0     1.2196e-09
     2.703e-10           0      2.703e-10
    1.5017e-10           0     1.5017e-10
    1.2331e-10           0     1.2331e-10

Several lsqnonneg entries appear to be exactly zero, while the other solutions are not exactly zero. For example, the eighth entry in each solution is:

fprintf('The entries are %d, %d, and %d.\n',sol.x(8),sol2.x(8),sol3.x(8))
The entries are 2.702960e-10, 0, and 2.702960e-10.

Perhaps the other quadprog and lsqlin algorithm would be a bit more accurate.

options = optimoptions('lsqlin','Algorithm','trust-region-reflective');
[sol,fval,exitflag] = solve(nonneglsq,'Solver','lsqlin','Options',options);
Solving problem using lsqlin.

Local minimum possible.

lsqlin stopped because the relative change in function value is less than the function tolerance.
disp(sol.x(8))
   2.2800e-15
options = optimoptions('quadprog','Algorithm','trust-region-reflective');
[sol3,fval3,exitflag3] = solve(nonneglsq,'Solver','quadprog','Options',options);
Solving problem using quadprog.

Local minimum possible.

quadprog stopped because the relative change in function value is less than the sqrt of the function tolerance, the rate of change in the function value is slow, and no negative curvature was detected.
disp(sol3.x(8))
   6.7615e-13

The solutions are slightly closer to zero using the 'trust-region-reflective' algorithm than the default interior-point algorithms.

Related Topics