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)

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'); [sol3,fval3,exitflag3] = solve(nonneglsq,'Solver','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);

Local minimum possible. lsqlin stopped because the relative change in function value is less than the function tolerance.

disp(sol.x(8))

2.2832e-15

options = optimoptions('quadprog','Algorithm','trust-region-reflective'); [sol3,fval3,exitflag3] = solve(nonneglsq,'Solver','quadprog','Options',options);

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.