fmincon; produces different answers against theoretically the same question ...

Hi all, I have a curiosity, I'm optimizing a non-linear function. There are some equality conditions (i.e the sum of the answers must equal 1), and then some boundaries 0 < x_i < 0.2.

There are eight variables. I can imagine two ways to solve the problem;

          weights = fmincon(@absexp,x0,A,b,Aeq,beq,[],[],[],options);
          weights2 = fmincon(@absexp,x0,[],[],Aeq,beq,lb,ub,[],options);

In the second case, I impose the constraints as bounds which are passed directly to fmincon.

In the first case, A and b represent inequalities which setup the same constraints... by saying that;

- x_i >= 0 and x_i <= 0.2

What's interesting is that the two solutions yield different results from the same starting conditions. Nearly 1% on x_1 ... so it's not negligible, and I believe it's way outside of the default tolerance for a solution that fmincon advertises as it's default (I think?).

Is this reasonable? Is there something here that would cause matlab to wheel in fundamentally different algorithms for the solution?

Any opinion greatly appreciated! Simon

Answer 1

[Alan Weiss]

Bound constraints are indeed handled differently than linear inequality constraints. You are likely to get a faster, more accurate solution using bounds than linear inequalities, but not always. That is one reason why the documentation recommends using bounds when possible instead of linear or nonlinear constraints.

As the documentation describes, several fmincon algorithms can satisfy bounds at every iteration. There is no such guarantee for linear constraints.

Alan Weiss

MATLAB mathematical toolbox documentation

[Matt J]

several fmincon algorithms can satisfy bounds at every iteration. There is no such guarantee for linear constraints.

Further to Alan's point, if your objective function is differentiable only in the interior of the region where the bounds are satisfied, then confining the algorithm's iterations to that interior region is the only way to ensure that the algorithm will converge properly.

I don't know what your @absexp function looks like, but if it looks something like this,

    f(x)=exp(abs(x))

then it is an example of what I'm talking about. If you specify explicit positivity bounds lb=0, then fmincon's default interior-point algorithm should ensure that x>0 at all iterations and avoid the non-differentiability (except asymptotically) at x=0. You might also need to use GradObj='on' to avoid finite differencing samples being chosen from the bad region x<=0.

[Nick]

Does this still hold up in R2021b (specifically for fmincon with interior point algorithm)? I notice that the link documenting that it is recommneded to use bounds when possible no longer exists.

[John D'Errico]

It is surely still a good idea.

[Matt J]

I could have sworn seeing an indication in the code that the toolbox now pre-processes A,b to see if any of the general linear inequality constraints are bound constraints. The following test indicates, though, that that is not the case,

 A=[-eye(2);eye(2)]; b=[0; 0; 1 ;1];

 opts=optimoptions('fmincon','Display','iter');
 
 fmincon(@(x) norm(x)^2,[10,5],A,b,[],[],[0,0],[1,1],[],opts)

Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to strictly satisfy the bounds.

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       3    6.050000e-01    0.000e+00    6.841e-01
    1       6    8.717252e-02    0.000e+00    4.341e-02    4.826e-01
    2       9    1.295851e-01    0.000e+00    9.322e-02    6.473e-02
    3      12    5.159471e-02    0.000e+00    2.753e-02    1.328e-01
    4      15    1.458103e-02    0.000e+00    7.387e-03    1.064e-01
    5      18    4.660130e-03    0.000e+00    2.381e-03    5.249e-02
    6      21    2.329639e-03    0.000e+00    1.200e-03    2.000e-02
    7      24    8.089295e-04    0.000e+00    4.083e-04    1.982e-02
    8      27    4.470690e-04    0.000e+00    2.266e-04    7.298e-03
    9      30    1.152813e-04    0.000e+00    5.779e-05    1.041e-02
   10      33    3.091217e-05    0.000e+00    1.546e-05    5.177e-03
   11      36    9.852413e-06    0.000e+00    4.931e-06    2.421e-03
   12      39    4.864689e-06    0.000e+00    2.435e-06    9.333e-04
   13      42    1.649849e-06    0.000e+00    8.253e-07    9.211e-04

Local 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.

ans = 1×2

1.0e+-3 * 0.9083 0.9083

 
 fmincon(@(x) norm(x)^2,[10,5],A,b,[],[],[],[],[],opts)

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       3    1.250000e+02    9.000e+00    2.000e+01
    1       6    1.514284e-01    0.000e+00    7.679e-01    1.097e+01
    2      10    1.123232e-01    0.000e+00    3.846e-01    2.945e-01
    3      13    7.168794e-02    0.000e+00    1.040e-01    8.443e-02
    4      16    3.241175e-02    0.000e+00    3.524e-02    8.772e-02
    5      19    8.872412e-03    0.000e+00    9.449e-03    8.585e-02
    6      22    2.432764e-03    0.000e+00    2.462e-03    4.487e-02
    7      25    7.102976e-04    0.000e+00    7.195e-04    2.267e-02
    8      28    2.900992e-04    0.000e+00    2.938e-04    9.619e-03
    9      31    9.423019e-05    0.000e+00    9.480e-05    7.325e-03
   10      34    4.766173e-05    0.000e+00    4.792e-05    2.803e-03
   11      37    1.216450e-05    0.000e+00    1.218e-05    3.416e-03
   12      40    3.244082e-06    0.000e+00    3.248e-06    1.687e-03
   13      43    1.023184e-06    0.000e+00    1.024e-06    7.896e-04
   14      46    4.947461e-07    0.000e+00    4.952e-07    3.081e-04

Local 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.

ans = 1×2

1.0e+-3 * 0.4973 0.4975

Answer 2

[Matt J]

I thought the Optimization Toolbox solvers now preprocess the linear constraints in A,b, Aeq,beq to see if any can be re-expressed as pure bounds, but apparently not.

In any case, the separateBounds() function from,

https://www.mathworks.com/matlabcentral/fileexchange/30892-analyze-n-dimensional-convex-polyhedra

will do so, e.g.,

A=[-eye(2);
    eye(2);
    1 1]; 

b=[0 0  1 1 1]';

Aeq=[0 1]; 
beq=[0.5];

[A,b,Aeq,beq,lb,ub]=separateBounds(A,b,Aeq,beq)

A = 1×2

1 1

b = 1

Aeq =

  0×2 empty double matrix


beq =

  1×0 empty double row vector

lb = 2×1

0 0.5000

ub = 2×1

1.0000 0.5000

[Nick]

Oh ok I see. Thank you so much!