MATLAB Examples

# Solver Behavior with a Nonsmooth Problem

This example shows the importance of choosing an appropriate solver for optimization problems. It also shows that a single point of non-smoothness can cause problems for Optimization Toolbox™ solvers.

In general, the solver decision tables provide guidance on which solver is likely to work best for your problem. For smooth problems, see docid:optim_ug.brhkghv-19. For nonsmooth problems, see docid:gads.bsa_e9p first, and for more information consult docid:gads.bu9r6el-1.

## A Function with a Single Nonsmooth Point

The function is nonsmooth at the point 0, which is the minimizing point. Here is a 2-D plot using the matrix norm for the 4-D point .

figure
x = linspace(-5,5,51);
[xx,yy] = meshgrid(x);
zz = zeros(size(xx));
for ii = 1:length(x)
for jj = 1:length(x)
zz(ii,jj) = sqrt(norm([xx(ii,jj),yy(ii,jj);0,0]));
end
end

surf(xx,yy,zz)
xlabel('x(1)')
ylabel('x(2)')
title('Norm([x(1),x(2);0,0])^{1/2}')

This example uses matrix norm for a 2-by-6 matrix x. The matrix norm relates to the singular value decomposition, which is not as smooth as the Euclidean norm. See docid:matlab_ref.bt00e3k-1.

## Minimize Using patternsearch

patternsearch is the recommended first solver to try for nonsmooth problems. See docid:gads.bsa_e9p. Start patternsearch from a nonzero 2-by-6 matrix x0, and attempt to locate the minimum of . For this attempt, and all others, use the default solver options.

Return the solution, which should be near zero, the objective function value, which should likewise be near zero, and the number of function evaluations taken.

fun = @(x)norm([x(1:6);x(7:12)])^(1/2);
x0 = [1:6;7:12];
rng default
x0 = x0 + rand(size(x0))
x0 =

1.8147    2.1270    3.6324    4.2785    5.9575    6.1576
7.9058    8.9134    9.0975   10.5469   11.9649   12.9706

[xps,fvalps,eflagps,outputps] = patternsearch(fun,x0);
xps,fvalps,eflagps,outputps.funccount
Optimization terminated: mesh size less than options.MeshTolerance.

xps =

1.0e-04 *

0.1116   -0.1209    0.3503   -0.0520   -0.1270    0.2031
-0.3082   -0.1526    0.0623    0.0652    0.4479    0.1173

fvalps =

0.0073

eflagps =

1

ans =

10780

patternsearch reaches a good solution, as evinced by exit flag 1. However, it takes over 10,000 function evaluations to converge.

## Minimize Using fminsearch

The documentation states that fminsearch sometimes can handle discontinuities, so this is a reasonable option.

[xfms,fvalfms,eflagfms,outputfms] = fminsearch(fun,x0);
xfms,fvalfms,eflagfms,outputfms.funcCount

Exiting: Maximum number of function evaluations has been exceeded
- increase MaxFunEvals option.
Current function value: 3.197063

xfms =

2.2640    1.1747    9.0693    8.1652    1.7367   -1.2958
3.7456    1.2694    0.2714   -3.7942    3.8714    1.9290

fvalfms =

3.1971

eflagfms =

0

ans =

2401

Using default options, fminsearch runs out of function evaluations before it converges to a solution. Exit flag 0 indicates this lack of convergence. The reported solution is poor.

## Use particleswarm

particleswarm is recommended as the next solver to try. See docid:gads.bu9r6hl-6.

[xpsw,fvalpsw,eflagpsw,outputpsw] = particleswarm(fun,12);
xpsw,fvalpsw,eflagpsw,outputpsw.funccount
Optimization ended: relative change in the objective value
over the last OPTIONS.MaxStallIterations iterations is less than OPTIONS.FunctionTolerance.

xpsw =

1.0e-12 *

Columns 1 through 7

-0.0386   -0.1282   -0.0560    0.0904    0.0771   -0.0541   -0.1189

Columns 8 through 12

0.1290   -0.0032    0.0165    0.0728   -0.0026

fvalpsw =

4.5222e-07

eflagpsw =

1

ans =

37200

particleswarm finds an even more accurate solution than patternsearch, but takes over 35,000 function evaluations. Exit flag 1 indicates that the solution is good.

## Use ga

ga is a popular solver, but is not recommended as the first solver to try. See how well it works on this problem.

[xga,fvalga,eflagga,outputga] = ga(fun,12);
xga,fvalga,eflagga,outputga.funccount
Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

xga =

Columns 1 through 7

-0.0061   -0.0904    0.0816   -0.0484    0.0799   -0.1925    0.0048

Columns 8 through 12

0.3581    0.0848    0.0232    0.0237   -0.1294

fvalga =

0.6257

eflagga =

1

ans =

68600

ga does not find as good a solution as patternsearch or particleswarm, and takes about twice as many function evaluations as particleswarm. Exit flag 1 is misleading in this case.

## Use fminunc from Optimization Toolbox

fminunc is not recommended for nonsmooth functions. See how it performs on this one.

[xfmu,fvalfmu,eflagfmu,outputfmu] = fminunc(fun,x0);
xfmu,fvalfmu,eflagfmu,outputfmu.funcCount
Local minimum possible.

fminunc stopped because the size of the current step is less than
the default value of the step size tolerance.

xfmu =

-0.5844   -0.9726   -0.4356    0.1467    0.3263   -0.1002
-0.0769   -0.1092   -0.3429   -0.6856   -0.7609   -0.6524

fvalfmu =

1.1269

eflagfmu =

2

ans =

429

The fminunc solution is not as good as the ga solution. However, fminunc reaches the rather poor solution in relatively few function evaluations. Exit flag 2 means you should take care, the first-order optimality conditions are not met at the reported solution.

## Use fmincon from Optimization Toolbox

fmincon can sometimes minimize nonsmooth functions. See how it performs on this one.

[xfmc,fvalfmc,eflagfmc,outputfmc] = fmincon(fun,x0);
xfmc,fvalfmc,eflagfmc,outputfmc.funcCount
Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints are
satisfied to within the default value of the constraint tolerance.

xfmc =

1.0e-10 *

-0.5368   -0.3821    0.1441   -0.6352    0.4227    0.0358
0.4836    0.1967    0.3806   -0.5431   -0.2634    0.2564

fvalfmc =

1.0051e-05

eflagfmc =

2

ans =

743

fmincon with default options produces an accurate solution after fewer than 1000 function evaluations. Exit flag 2 does not mean that the solution is inaccurate, but that the first-order optimality conditions are not met. This is because the gradient of the objective function is not zero at the solution.

## Summary of Results

Choosing the appropriate solver leads to better, faster results. This summary shows how disparate the results can be. The solution quality is 'Poor' if the objective function value is greater than 0.1, 'Good' if the value is smaller than 0.01, and 'Mediocre' otherwise.

Solver = {'patternsearch';'fminsearch';'particleswarm';'ga';'fminunc';'fmincon'};
SolutionQuality = {'Good';'Poor';'Good';'Poor';'Poor';'Good'};
FVal = [fvalps,fvalfms,fvalpsw,fvalga,fvalfmu,fvalfmc]';
NumEval = [outputps.funccount,outputfms.funcCount,outputpsw.funccount,...
outputga.funccount,outputfmu.funcCount,outputfmc.funcCount]';
results = table(Solver,SolutionQuality,FVal,NumEval)
results =

6x4 table

Solver         SolutionQuality       FVal       NumEval
_______________    _______________    __________    _______

'patternsearch'    'Good'              0.0072656    10780
'fminsearch'       'Poor'                 3.1971     2401
'particleswarm'    'Good'             4.5222e-07    37200
'ga'               'Poor'                0.62572    68600
'fminunc'          'Poor'                 1.1269      429
'fmincon'          'Good'             1.0051e-05      743

Another view of the results.

figure
hold on
for ii = 1:length(FVal)
clr = rand(1,3);
plot(NumEval(ii),FVal(ii),'o','MarkerSize',10,'MarkerEdgeColor',clr,'MarkerFaceColor',clr)
text(NumEval(ii),FVal(ii)+0.2,Solver{ii},'Color',clr);
end
ylabel('FVal')
xlabel('NumEval')
title('Reported Minimum and Evaluations By Solver')
hold off

While particleswarm achieves the lowest objective function value, it does so by taking over three times as many function evaluations as patternsearch, and over 30 times as many as fmincon.

fmincon is not generally recommended for nonsmooth problems. It is effective in this case, but this case has just one nonsmooth point.