Main Content

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

$$\underset{x}{\mathrm{min}}\Vert Cx-d{\Vert}^{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)

`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: 'active-set', '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,output] = 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.

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

disp(output)

message: '...' algorithm: 'interior-point' firstorderopt: 9.9673e-07 constrviolation: 0 iterations: 9 linearsolver: 'sparse' cgiterations: [] solver: 'lsqlin'

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.

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);

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

disp(output2)

iterations: 14 algorithm: 'trust-region-reflective' firstorderopt: 5.2187e-08 cgiterations: 54 constrviolation: [] linearsolver: [] message: 'Local minimum possible....' solver: 'lsqlin'

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);

Solving problem using lsqlin. Optimal solution found.

disp(output3)

iterations: 11 algorithm: 'trust-region-reflective' firstorderopt: 1.3973e-14 cgiterations: 0 constrviolation: [] linearsolver: [] message: 'Optimal solution found.' solver: 'lsqlin'

Using this option brings the first-order optimality measure nearly to zero, which is the best possible.

The `lsqnonneg`

solver is specially designed to handle nonnegative linear least squares. Try that solver.

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

Solving problem using lsqnonneg.

disp(output4)

iterations: 184 algorithm: 'active-set' message: 'Optimization terminated.' solver: 'lsqnonneg'

`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'})

`t=`*1×4 table*
default trust-region-reflective factorization lsqnonneg
__________ _______________________ _____________ __________
5.0804e-05 4.9179e-05 4.9179e-05 4.9179e-05

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])

1.0e-12 * 0.7070 0.6999

The `lsqnonneg`

residual norm is the smallest by a nearly negligible amount. However, `lsqnonneg`

takes the most iterations to converge.