Quadprog 'interior-point-convex' failure

Dear Matlab Gurus,

I have a problem with Matlab quadprog function.

A few years ago I implemented a code using quadprog with the 'active-set' algorithm, and it worked pretty well. The problem is characterised by a structurally dense hessian and inequality constraint matrix (see data attached).

Now this algorithm is no more available and only 'interior-point-convex' and 'trust-region-reflective' can be used. However, the 'trust-region-reflective' cannot be used with inequality constraints. So the choice is reduced to 'interior-point-convex' that is not able to solve the problem:

clear variables, close all

% load problem
load('testData.mat');

% quadprog setup (R2018b)
opts = optimoptions(@quadprog,...
    'Algorithm','interior-point-convex');
tic
[x,fval,exitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,[],[],x0,opts);
toc

The output is

The interior-point-convex algorithm does not accept an initial point.
Ignoring X0.

No feasible solution found.

quadprog stopped because it was unable to find a point that satisfies
the constraints within the default value of the constraint tolerance.

<stopping criteria details>

The funny is that I already have a feasible solution x0, but the 'interior-point-convex' algorithm does not accept an initial point:

% feasibility check
any(A*x0-b > 0)

ans =

     0

When using the 'active-set' algorithm with a Matlab version < 2016 everything is fine and the algorithm ends correctly.

clear variables, close all

% load problem
load('testData.mat');

% quadprog setup (R2015b)
opts = optimoptions(@quadprog,...
    'Algorithm','active-set');
tic
[x,fval,exitflag,output,lambda] = quadprog(H,f,A,b,Aeq,beq,[],[],x0,opts);
toc

with output

Warning: The 'active-set' algorithm will be removed in a future release. To avoid this warning or a future error, choose a
different algorithm: 'interior-point-convex' or 'trust-region-reflective'.
 
> In quadprog (line 427)
  In test_quadprog_solvers (line 10) 
Optimization terminated.
Elapsed time is 2.239872 seconds.

I have seen that fmincon still have the 'active-set' algorithm and I am wondering if with a proper setup of options I can trick the optimizer to solve a purely quadratic programming. I tried:

clear variables, close all

% load problem
load('testData.mat');

% fmincon setup (R2018b)
opts = optimoptions(@fmincon,...
    'Algorithm','active-set',...
    'SpecifyObjectiveGradient',true);

tic
[xred,fval,exitflag,output,lambda] = fmincon(@(x)quadFcn(x,H,f),x0,A,b,Aeq,beq,[],[],[],opts);
toc

With

function [fval,gval] = quadFcn(x,H,f)

fval = .5*x'*H*x+f'*x;

if nargout > 1 % gradient required
    gval = H*x+f;
end
end

And the result is fine

Local minimum possible. Constraints satisfied.

fmincon stopped because the predicted change in the objective function
is less than the default value of the function tolerance and constraints 
are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

Elapsed time is 3.966116 seconds.

But the calculation time is almost doubled with respect to the quadprog function with 'active-set' algorithm, and things are worse for more complex problems than this simple benchmark.

Any idea to restore the original behavior or performance?

Thank you in advance

[军 卫]

THANKS！！

Answer 1

[Matt J]

Simply l2-normalizing the rows of A also seems to help,

load('testData-Matt.mat');

opts = optimoptions(@quadprog,'Algorithm','interior-point-convex');

anorm=vecnorm(A,2,2);
AAA=A./anorm;
bbb=b./anorm;

tic;
[y,fval,exitflag,output,lambda] = quadprog(H,f,AAA,bbb,[],[],[],[],[],opts);
toc%Elapsed time is 1.533334 seconds.

>> exitflag

exitflag =

     1

[Fabio Freschi]

This seems to solve the original problem. These are the results of the tests:

Algorithm 1: R2015b - quadprog with 'active-set' and inital guess

• fval = 2.341611223732766e+08
• time = 2.525043 seconds;

Algorithm 2: R2018b - quadprog with 'interior-point-convex' and orthonormalization of rows of A

• fval = 2.341611223732780e+08
• time = 1.465645 seconds

Algorithm 3: R2018b - fmincon with 'active-set' and initial guess

• fval = 2.341611223732788e+08
• time = 4.290164 seconds

I'll come back later with additional results for more complex cases.

Many thanks to Bruno and Matt for the helpful comments and quick replies.

[Matt J]

Just keep in mind that the fastest observed time was with the technique of translating and normalizing b(i) as shown here,

https://www.mathworks.com/matlabcentral/answers/467328-quadprog-interior-point-convex-failure#comment_715290

This should always be advantageous because, when given a problem with all b(i)>0, quadprog is probably smart enough to know right away that x=0 will always be an interior point.

[Fabio Freschi]

Thanks once again. I am going to make some extensive tests.

It is still puzzling me why quadprog does not make scaling by itself. It is like asking the user to provide pivoting before calling the backslash operator.

[Bruno Luong]

"It is still puzzling me why quadprog does not make scaling by itself."

Because to do it in the "right" way, user can do more than simple centering and scaling, such as linear transformation (or even non-linear in case of fmincon) of the input variables. Mathematically it's call right preconditioning of the problem.

Just keep in mind that when you do numerical minimization, precondingtioning is an important aspect that can break or make your problem resolution.

Answer 2

[Bruno Luong]

You can shift the input variable, it seems to work

[y,fval,exitflag,output,lambda] = quadprog(H,f,A,b-A*x0,[],[],[],[],[],opts);
x=x0+y

[Matt J]

The result is almost indistinguishable numerically from x0.

>>norm(x-x0,inf)./norm(x0,inf)*100

ans =

   5.8197e-16

[Bruno Luong]

Yes but at least it *find* a solution, which is obviolusly x0 feed in by OP.

[Matt J]

@Fabio,

In addition to what Bruno mentions, you should rescale your constraints so that all b(i)=1. This noticeably speeds things up. Note that I've chosen a random objective vector f here because the one in testData.mat is just zero.

f=rand(size(f));

opts = optimoptions(@quadprog,'Algorithm','interior-point-convex');

bb=b-A*x0;
AA=A./bb;
bb=ones(size(bb));

tic
[y,fval,exitflag,output,lambda] = quadprog(H,f,A,b-A*x0,[],[],[],[],[],opts);
toc%Elapsed time is 15.233883 seconds - and still non-convergent

tic;
[y,fval,exitflag,output,lambda] = quadprog(H,f,AA,bb,[],[],[],[],[],opts);
toc %Elapsed time is 0.604501 seconds.

[Fabio Freschi]

First and foremost thanks for your answers.

Some additional insight about the problem:

1) My problem is not to use the initial guess. I am inclined to believe that this could be useful because for the specific problem I always have an initial feasible guess. The problem is to find the solution for the quadratic problem, that is the solution that minimize the energy of a magnetic system (the matrix H is, in fact, the inductance matrix). I am not picky with the algorithm: any effective and efficient solution is fine, with or without the initial guess.

2) I don't think x0 is the optimal solution, with exact arithmetic. x0 is for sure feasible but it is not optimal. Maybe the conditioning of the problem is such that with the given tolerances, x0 is the optimal solutions of quadprog.

3) f is a vector of zeros because it is not used in the formulation of the minimization problem. I am keeping the vector for possible future use, but it will be unlikely used.

Answer 3

[Matt J]

I don't have a pre-R2016a Matlab release readily at hand, and so I cannot compare performance. However, if the problem is simply that you would like to be able to incorporate the initial x0 into quadprog's interior-point-convex algorithm, you could try adding one more inequality constraint,

f.'*x <= f.'*x0

With this constraint, quadprog should restrict its search to points with objective better than x0.

[Bruno Luong]

How adding even more constraint can help error "No feasible solution found." ? Obviously the interior-point-convex algo fails to find the feasible point, possibly because it's more sensitive to whatever numerical instabily or error tan the active set algo.

[Matt J]

If the problem is infeasible, it means that x0 is already optimal.

[Bruno Luong]

??? The problem is feasible, just QUADPROG fails to find it, if you read OP description.

[Bruno Luong]

Let me try it for you Matt:

>> [x,fval,exitflag,output,lambda] = quadprog(H,f,[A],[b],[],[],[],[],[],opts);

No feasible solution found.

quadprog stopped because it was unable to find a point that satisfies
the constraints within the value of the constraint tolerance.

<stopping criteria details>
>> [x,fval,exitflag,output,lambda] = quadprog(H,f,[A;f'],[b;f'*x0],[],[],[],[],[],opts);

No feasible solution found.

quadprog stopped because it was unable to find a point that satisfies
the constraints within the value of the constraint tolerance.

<stopping criteria details>

% However

>> all(A*x0<=b)

ans =

  logical

   1

>> all([A;f']*x0<=[b; f'*x0])

ans =

  logical

   1

>> 

[Matt J]

OK, my approach doesn't address the difficulty of not being able to find an initial feasible point to the posted problem, but it does address the more general issue that the interior-point quadprog implementation doesn't let you benefit from a prior x0. If you have a good intial guess, one should be able to exploit that to give the search a head start and reduce computation...

[Bruno Luong]

You might convince me if you show if there is any saving in the CPU time. Not sure adding this constraint help QUADPROG to better start or it make it even worse by dealing with more constraints.

[Matt J]

Here is a case where it helped CPU time, though perhaps it is difficult to guarantee.

load('testData-Matt.mat');

anorm=vecnorm(A,2,2);
AAA=A./anorm;
bbb=b./anorm;

tic;
[y,fval,exitflag,output,lambda] = quadprog(H,f,AAA,bbb,[],[],[],[],[],opts);
toc %Elapsed time is 1.527962 seconds.



fff=f./norm(f);
bbbf=fff.'*x0;
AAAA=[AAA;fff.'];
bbbb=[bbb;bbbf];


tic;
[y,fval,exitflag,output,lambda] = quadprog(H,f,AAAA,bbbb,[],[],[],[],[],opts);
toc%Elapsed time is 1.068352 seconds.