quadprog different output for R2020a and R2017a

Hello,

If I run quadprog minimization function I get completely different results on R2020a and R2017a. It seems that the output on the lattest release is wrong. Did something changed? Is it a bug you are aware of?

Kind regards

[Ameer Hamza]

Release notes do not mention any major change from R2017a to R2020a: https://www.mathworks.com/help/matlab/release-notes.html?rntext=quadprog&startrelease=R2017a&endrelease=R2020a&groupby=release&sortby=descending&searchHighlight=quadprog

[Jason Nicholson]

Which algorithm are you using? The Inter-Point Algorithm has problems with ill-conditioned medium size problems. The Active-set algorithm in older version of MATLAB could handle this problem but suffered from other issues such that it was deprecated. The Trust region algorithm is what you should try if you can.

[Kyril Kaufmann]

Hello,

Interor-point-convex in both releases.

If I run this code

% Load files
structure = load('quadprog_input.mat')
H = structure.quadprog_input.H;         % 404x404 matrix
f = structure.quadprog_input.f;         % 1x404 vector
Aineq = structure.quadprog_input.Aineq; % 808x404 matrix
bineq = structure.quadprog_input.bineq; % 808x1 vector
opt   = structure.quadprog_input.opt;
% with opt.TolCon = 100*eps, opt.TollFun = 100*eps, opt.Display = 'none', 
% opt.Algorithm = 'interior-point-convex
Aeq = [];
beq = [];

% solve quadprog problem
[dp,fval,~] = quadprog(2*H,f,...
  Aineq, bineq, Aeq, beq,[],[],[],opt)

I get completely different output. I inserted the '.mat' files.

[Jason Nicholson]

Okay. I will take a look. My hunch is your H has a high condition number. I'll explain more if that is the issue.

[Kyril Kaufmann]

That's really nice of you.

For the moment I run the code with the 2017 release and it gives me nice results so it's not urgent. Thanks so much

Answer 1

[Jason Nicholson]

There is an option called 'LinearSolver'. It can be dense or sparse. I set it to sparse and it converged quickly.

This problem is solvable but be careful with the condition number of the H matrix. i.e. cond(H). The higher the condition number, the more ill-conditioned the problem. ill-condition problems are harder to solve. With the generic cost: J = 1/2*x'*H*x+f'*x. Small pertubations to f will cause large changes to the solution, x.

% Load files
structure = load('quadprog_input.mat');
H = structure.quadprog_input.H;         % 404x404 matrix
f = structure.quadprog_input.f;         % 1x404 vector
Aineq = structure.quadprog_input.Aineq; % 808x404 matrix
bineq = structure.quadprog_input.bineq; % 808x1 vector
% opt   = structure.quadprog_input.opt;
% with opt.TolCon = 100*eps, opt.TollFun = 100*eps, opt.Display = 'none',
% opt.Algorithm = 'interior-point-convex
Aeq = [];
beq = [];


f = f'; % f should be 404x1

% Objective is 1/2*dp'*2*H*dp+f'*dp
%                  
% dp' is 1x404
% dp is 404x1
% 2*H is 404x404
% f is 404x1
% f' is 1x404
%
%   1/2*dp' * 2*H   *dp   + f'   *dp
%      1x404 404x404 404x1  1x404 404x1
%      1x1                + 1x1
%      1x1
% Thus, f should be 404x1
%

opt = optimoptions('quadprog', 'TolCon', 100*eps, 'TolFun', 100*eps, ... 
     'Display', 'iter-detailed', 'Algorithm', 'interior-point-convex', ...
     'MaxIter', 1500,'LinearSolver','sparse');


% Best case cost ignoring constraints
[~,fval,~] = quadprog(2*H,f,[], [], [], [],[],[],[],opt)

% solve quadprog problem. Note this doesn't converge
[dp,fval,~] = quadprog(2*H,f,Aineq, bineq, Aeq, beq,[],[],[],opt); fval

[Jason Nicholson]

FYI, if this is a linear fitting problem that you have converted to quadratic program. I have another option for you. However, I will need your C matrix in H=C'*C.

[Kyril Kaufmann]

Hello, thanks for your support. No it's not a linear fitting.

With your set of options it works! At least my power imput (dp) is within meaningfull physical range. I'm not sure though what happend here...