Quadratically constrained linear maximisation problem: issues with fmincon

I would like to solve the following quadratically constrained linear programming problem.

I have written a Matlab code (R2091b) that solve the problem using Gurobi. Now, I would like to rewrite the code using fmincon instead of Gurobi. This is because the optimisation problem will have to be solve thousands of times and Gurobi's academic license does not allow to parallelise via array jobs in a cluster. However, I'm encountering a huge problem: Gurobi takes 0.23 second to give a solution, fmincon takes 13 sec. I suspect this should be due to my mistakes/inefficiences in providing gradient, hessian, etc. Could you kindly help me to improve below?

Also, Gurobi gives me 0.2 as solution, fmincon gives 0.089. Can the accuracy of fmincon be improved without trying other starting points?

This is my code with Gurobi

    clear 
    rng default
    load matrices
    
    %1) GUROBI
    model.A=[Aineq; Aeq];
    model.obj=f;
    model.modelsense='max';
    model.sense=[repmat('<', size(Aineq,1),1); repmat('=', size(Aeq,1),1)];
    model.rhs=[bineq; beq]; 
    model.ub=ub;
    model.lb=lb;
    model.quadcon(1).Qc=Q; 
    model.quadcon(1).q=q;
    model.quadcon(1).rhs=d;
    params.outputflag = 0; 
    result=gurobi(model, params);
    max_problem_Gurobi=result.objval;

This is my code with fmincon

    %2) FMINCON
    options = optimoptions(@fmincon,'Algorithm','interior-point',...
                                             'SpecifyObjectiveGradient',true,...
                                             'SpecifyConstraintGradient',true,...
                                             'HessianFcn',@(z,lambda)quadhess(z,lambda,Q));
    fun = @(z)quadobj(z,f.');
    nonlconstr = @(z)quadconstr(z,Q,d);
    [~,fval] = fmincon(fun,z0,Aineq,bineq,Aeq,beq,lb,ub,nonlconstr,options);
    max_problem_fmincon=-fval; 

    function [y,yeq,grady,gradyeq] = quadconstr(z,Q,d)
    y= z'*Q*z -d;
    yeq = []; %no quadratic inequalities 
    if nargout > 2
        grady = 2*Q*z;
    end
    gradyeq = [];  %no quadratic inequalities
    end

    function hess = quadhess(z,lambda,Q) %#ok<INUSL>
    hess = 2*lambda.ineqnonlin*Q;
    end
    
    function [y,grady] = quadobj(z,f)
    y = -f'*z;
    if nargout > 1
        grady = -f;
    end

Further, if I run the code with fmincon with as starting point the optimal point given by Gurobi, I still get the solution 0.089 (instead of 0.2 as in Gurobi). Why?

    %3) FMINCON
    z0=result_Gurobi.x;
    options = optimoptions(@fmincon,'Algorithm','interior-point',...
                                             'SpecifyObjectiveGradient',true,...
                                             'SpecifyConstraintGradient',true,...
                                             'HessianFcn',@(z,lambda)quadhess(z,lambda,Q));
    fun = @(z)quadobj(z,f.');
    nonlconstr = @(z)quadconstr(z,Q,d);
    [~,fval] = fmincon(fun,z0,Aineq,bineq,Aeq,beq,lb,ub,nonlconstr,options);
    max_problem_fmincon=-fval; 

    function [y,yeq,grady,gradyeq] = quadconstr(z,Q,d)
    y= z'*Q*z -d;
    yeq = []; %no quadratic inequalities 
    if nargout > 2
        grady = 2*Q*z;
    end
    gradyeq = [];  %no quadratic inequalities
    end

    function hess = quadhess(z,lambda,Q) %#ok<INUSL>
    hess = 2*lambda.ineqnonlin*Q;
    end
    
    function [y,grady] = quadobj(z,f)
    y = -f'*z;
    if nargout > 1
        grady = -f;
    end
    end

[Torsten]

grady = 2*Q*z

Same for hess.

[CT]

Thanks, question updated. But it still runs forever.

[Torsten]

In quadobj, y=-f'*x and grady = -f since you want to maximize.

[CT]

Thanks!!! Now it runs in 13 sec (still too much). Also it gives as solution 0.089 (while Gurobi gives 0.2). I understand I should try several starting points but it is not feasible to do that in my real exercise.

[Torsten]

What if you don't specify derivatives ?

[CT]

I tried this (without calling options)

fun = @(z)quadobj(z,f.');
nonlconstr = @(z)quadconstr(z,Q,d);
[~,fval] = fmincon(fun,z0,Aineq,bineq,Aeq,beq,lb,ub,nonlconstr);
max_problem_fmincon=-fval; 

but it tells me

Solver stopped prematurely.
fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+03.

When I increase MaxFunctionEvaluations to e.g. 10^6 it runs forever.

[Torsten]

Strange.

[CT]

Thanks: is it working normally in your computer?

[CT]

Any further help?

[Torsten]

I don't have access to Matlab at the moment.

Did you check whether both solutions (Gurobi and Matlab) satisfy the constraints ?

Did you try a different solver in Matlab (e.g. sqp) ?

[CT]

Thanks.

1) Did you check whether both solutions (Gurobi and Matlab) satisfy the constraints?

For the equality constraints, fmincon satisfies the constraints with more precision. For example, suppose a constraint requires , then fmincon is such that while Gurobi is such that

The inequality constraints are satisfied by bothe solvers.

2) Did you try a different solver in Matlab (e.g. sqp) ?

sqp runs forever.

[Torsten]

What about the quadratic constraint and the bound constraints ?

[CT]

Inequality constraints are satisfied by both solvers.

Quadratic constraint satisfied by both solvers: with Gurobi , with fmincon

[Matt J]

I suggest you include the results from Gurobi in the matrices.mat file as well, so we can study it.

[CT]

Thanks: they are now in the matrices under the name result_Gurobi

[Matt J]

I suspect this should be due to my mistakes/inefficiences in providing gradient, hessian, etc.

I don't think that is the reason. I think the main reason is that Gurobi apparently has specific mechanisms for handling quadratic constraints, as suggested by this segment of code.

    model.quadcon(1).Qc=Q; 
    model.quadcon(1).q=q;
    model.quadcon(1).rhs=d;

Conversely, fmincon has no way of distinguishing quadratic constraints from other more general nonlinear constraints and giving them special handling. It handles all nonlinear constraints in the same way.

That said, the following implementation of your functions would be slightly more efficient.

    function [y,yeq,grady,gradyeq] = quadconstr(z,Q,d)
    
        Qz=Q*z;
        y= z'*Qz - d;
        yeq = []; %no quadratic inequalities 
        if nargout > 2
            grady = 2*Qz;
            gradyeq = [];  %no quadratic inequalities
        end
        
    end

    function [y,grady] = quadobj(z,f)
     grady = -f;
     y = grady.'*z;
    end

[CT]

Thanks, but it is not more efficient, still 36 sec.

[Matt J]

There was no expectation of great gains, but I think it has to be marginally more efficient, maybe a reduction from 36.05 sec to 36 sec. You can plainly see that fewer vector arithmetic operations are done in this version.

Answer 1

[Matt J]

Well, it would be interesting to know what algorithm Gurobi uses, but the issue of the objective function difference appears to be a matter of the tolerances

    options = optimoptions(@fmincon,'Algorithm','interior-point',...
                                             'SpecifyObjectiveGradient',true,...
                                             'SpecifyConstraintGradient',true,...
                                             'HessianFcn',@(z,lambda)quadhess(z,lambda,Q),...
                                             'StepTolerance',1e-30,'OptimalityTolerance',1e-10);
    fun = @(z)quadobj(z,f);
    nonlconstr = @(z)quadconstr(z,Q,d);
    tic;
    [~,fval] = fmincon(fun,z0(:),Aineq,bineq,Aeq,beq,lb,ub,nonlconstr,options);
    toc
    max_problem_fmincon=-fval

    max_problem_fmincon =

    0.2000

[CT]

Thanks. Time difference is huge though: Gurobi 0.23 sec, fmincon 34 sec. Is there a way to improve on fmincon's speed?

[Matt J]

Maybe transform the problem into a space where Q is diagonal. Also, maybe impose a quadratic equality constraint instead of inequality. It is easy to verify in advance whether the solution satisfies the quadratic constraint with strict inequality. The optimality conditions in that case are the same as for when the quadratic constraint is absent, leaving only the linear constraints. You can therefore use linprog, and see if the solution it gives happens to satisfy the quadratic constraints.

[CT]

Thanks. Can you elaborate on that? My question was indeed about improving the speed of fmincon in my specific example.

1) Transform the problem into a space where Q is diagonal: how? can you link some examples?

2) Your second suggestion is to run the following LP

Let be the solution of the LP. Then, if , stop here. If , run

Correct?

[Matt J]

1) Your Q is symmetric with eigendecomposition Q=R*D*R.', so if we make the change of variables x=R.'*z, the problem will have the same form in terms of x (linear objective and constraints), and the quadaratic constraint will become

Since D is diagonal, this may be simpler for an iterative solver to try to satisfy.

2) Correct.

[Matt J]

Maybe transform the problem into a space where Q is diagonal.

On 2nd thought, that will probably ruin the sparsity of your linear constraints.

[CT]

"ruin"= reduce speed?

[Matt J]

Yes, I think the sparsity is important for speed.

[Matt J]

How fast does Gurobi solve the problem when the quadratic inequality is replaced with equality? And how fast when the quadratic constraint is removed altogether?

[CT]

1) With quadratic equality: I need to investigate because it is non-convex

2) Without quadratic constraint: 0.16 sec.

[Matt J]

And does the problem data from the thousands of problem instances that you are trying to solve change in a continuous incremental way? If you had the optimal solution for one instance of the problem, would it serve as a good initial estimate for the next instance?