MATLAB: Quadprog and fmincon provide very different results on the same problem

fminconMATLABquadprog

Dear all,
I am trying to solve a quadratic minimization problem with quadprog, but I obtain results that are very different from fmincon. I do not understand why. My problem is the following:
where || || is the Euclidean norm. I have written the following code (with quantities similar to my problem):
% Setup values (just as example)
n = 4;
AA = [12, -1, 0, 0; ...
0, 0, 111, -1];
bb = [0; 0];
x = [1; 10; 1; 61];
% options for fmincon
optionsFmincon = optimset('LargeScale', 'off', ...
'Algorithm', 'interior-point', ...
'UseParallel', 'never', ...
'MaxFunEvals', 20000, ...
'GradObj', 'on', ...
'TolX', 1e-9, ...
'TolFun', 1e-9, ...
'TolCon', 1e-9, ...
'MaxIter', 10000, ...
'Display', 'none' ...
);
% options for quadprog
optionsQuadprog = optimset('Algorithm', 'interior-point-convex', ...
'TolX', 1e-12, ...
'TolFun', 1e-12, ...
'TolCon', 1e-12, ...
'MaxIter', 10000, ...
'Display', 'none' ...
);
% Fmincon minimization
lowerBound = -1000000000000000*ones(params.n, 1); %fictitious lower and upper bounds
upperBound = +1000000000000000*ones(params.n, 1);
[out1, Jopt1, exitflag1, output1] = fmincon(@(y) cost(x, y, n), x, [], [], [], [], lowerBound, upperBound, @(y) constraints(y, AA, bb), optionsFmincon);
% Solution with quadprog
H = eye(n);
f = -2*eye(n)*x;
A = [];
b = [];
Aeq = AA;
beq = bb;
lb = [];
ub = [];
x0 = [];
[out, Jopt, exitflag, output] = quadprog(H, f, A, b, Aeq, beq, lb, ub, x0, optionsQuadprog);
% Comparison of results
[out1, out]
[out1'*H*out1-2*x'*eye(n)*out1+x'*eye(n)*x, out'*H*out-2*x'*eye(n)*out+x'*eye(n)*x]
[norm(x-out1)^2, norm(x-out)^2]
[norm(AA*out1-bb), norm(AA*out-bb)]
% cost function for fmincon
function [J, gradient] = cost(x, y, n)
J = (x-y)'*eye(n)*(x-y); % equivalent to J = y'*eye(n)*y - (2*eye(n)*x)'*y; since I can omit x*eye(n)*x
gradient = -2*(x-y);
end
% constraints function for fmincon
function [c, ceq] = constraints(x, AA, bb)
ceq = AA*x-bb;
c = [];
end
If you execute the code, you get very different solutions. The optimal cost of the fmincon is equal to 0.23048, while the optimal cost of quadprog is 3823 (therefore, also the values of the optimal y are very different). In both cases the constraints are satisfied.
Probably I am making something wrong, but I do not understand what.
Any help from you is really appreciated. Thank you in advance

Best Answer

  • Why are you doing this? Really, why are you trying to solve it using that kluge of approaches? Go simple. Don't over-complicate things.
    n = 4;
    AA = [12, -1, 0, 0; ...
    0, 0, 111, -1];
    bb = [0; 0];
    x = [1; 10; 1; 61];
    Lets see. Minimize norm(x-y), subject to AA*y = bb. Trivial, using the proper tool. In the optimization TB, that is lsqlin. As I show below, my own LSE also gives it directly, found on the file exchange.
    Ylsqlin = lsqlin(eye(4),x,[],[],AA,bb)
    Minimum found that satisfies the constraints.
    Optimization completed because the objective function is non-decreasing in
    feasible directions, to within the default value of the optimality tolerance,
    and constraints are satisfied to within the default value of the constraint tolerance.
    <stopping criteria details>
    Ylsqlin =
    0.834482758578957
    10.0137931029475
    0.549586106124114
    61.0040577797767
    Verification:
    AA*Ylsqlin - bb
    ans =
    -1.77635683940025e-15
    0
    The sum of squares of the difference was:
    norm(x - Ylsqlin)^2
    ans =
    0.230475348269706
    Using my own LSE for comparison, found on the file exchange.
    lse(eye(4),x,AA,bb)
    ans =
    0.83448275862069
    10.0137931034483
    0.549586106151599
    61.0040577828275
    Now, how would it be solved using quadprog? The sum of squares of differences is:
    (x-y)'*eye(4)*(x-y)
    Expanding, we get
    x'*x + y'*y - 2*x'*y
    x'*x is a constant, so irrelevant to the optimization, BUT it explains why your objectives are so different.
    x'*x
    ans =
    3823
    To convert this objective to a quadprog usable form, we want it as:
    FVAL = 0.5*X'*H*X + f'*X.
    So
    H = 2*eye(4);
    F = -2*x';
    Now solve.
    Yqp = quadprog(H,F,[],[],AA,bb)
    Minimum found that satisfies the constraints.
    Optimization completed because the objective function is non-decreasing in
    feasible directions, to within the default value of the optimality tolerance,
    and constraints are satisfied to within the default value of the constraint tolerance.
    <stopping criteria details>
    Yqp =
    0.834482758599826
    10.0137931031979
    0.549586106137858
    61.0040577813022
    It should be no surprise this is the same as lsqlin, but it took more thought to write. Use the proper tools. And now, you want to use fmincon.
    fun = @(y) norm(x-y);
    x0 = rand(4,1);
    Yfmincon = fmincon(fun,x0,[],[],AA,bb)
    Local minimum found that satisfies the constraints.
    Optimization completed because the objective function is non-decreasing in
    feasible directions, to within the default value of the optimality tolerance,
    and constraints are satisfied to within the default value of the constraint tolerance.
    <stopping criteria details>
    Yfmincon =
    0.834482752432331
    10.013793029188
    0.549586102080044
    61.0040573308849
    So esentially identical. And very easy to write in all cases, though lsqlin was by far the simplest.
    Anyway, there was absolutely NO need to provide a gradient function to fmincon. (In fact, there is almost never a need to provide a gradient.) A waste of time certainly in this case. Nor even an m-file cost function. The other differences in your results I would put down to simple algebra mistakes that are hardly worth debugging, because you massively over-complicated the problem.