Why does the QUADPROG function in the Optimization Toolbox sometimes return an incorrect result when solving a degenerate problem?

I want to solve the following problem using the QUADPROG function:

% Minimize 2*x(2)^2

H = zeros(4, 4); H(2,2) = 4;

f = zeros(4, 1);

% Subject to -x(3)+x(4)<=-1 and -2*x(2)+x(3)-x(4)<=-1

A = [0 0 -1 1; 0 -2 1 -1];

b = [-1;-1];

% Create default options structure

options = optimset('quadprog'); 

% Call quadprog

[x, fval] = quadprog(H, f, A, b, [], [], [], [], [], options);

From this code, I receive a result that does not satisfy the constraints, and the first component of "x" has orders of magnitude larger than the other components.

Answer 1

[MathWorks Support Team]

This enhancement has been incorporated in Release 2011a (R2011a). For previous product releases, read below for any possible workarounds:

The QUADPROG function within the Optimization Toolbox may return improper constraints for certain problems that result in an indefinite or singular Hessian matrix. The example below illustrates such a problem:

% Minimize 2*x(2)^2

H = zeros(4, 4); H(2,2) = 4;

f = zeros(4, 1);

% Subject to -x(3)+x(4)<=-1 and -2*x(2)+x(3)-x(4)<=-1

A = [0 0 -1 1; 0 -2 1 -1];

b = [-1;-1];

% Create default options structure

options = optimset('quadprog'); 

% Call quadprog

[x, fval] = quadprog(H, f, A, b, [], [], [], [], [], options);

% Evaluate results

disp(sprintf('The function value is %f at x = [%f;%f;%f;%f]', fval, x))

Constraint = A*x;

disp(sprintf('Constraint 1 is %f <= %f and constraint 2 is %f <= %f',...

     Constraint(1), b(1), Constraint(2), b(2)))

In this example, the Hessian matrix that QUADPROG computes is indeed singular. The singular Hessian and zero "f" vector do not provide QUADPROG with enough information to determine a good direction in which to search for the solution of the problem. There are two ways to modify the problem so that QUADPROG can compute a nonsingular Hessian and minimizer that satisfy the constraints:

1) Remove the degenerate variable or variables from the problem statement. In the example above, the variable "x(1)" does not appear in either the expression to be minimized or the constraints. Therefore, if you remove the variable "x(1)" from the problem entirely, you remove enough of the degeneracy that QUADPROG is able to determine a good enough search direction to solve the problem. Modifying the example above to do this yields:

% Minimize 2*x(2)^2

H = zeros(4, 4); H(2,2) = 4;

f = zeros(4, 1);

% Subject to -x(3)+x(4)<=-1 and -2*x(2)+x(3)-x(4)<=-1

A = [0 0 -1 1; 0 -2 1 -1];

b = [-1;-1];

% Now remove variable 1 from the problem

H(1,:) = [];

H(:,1) = [];

f(1) = [];

A(:,1) = [];

% Create default options structure

options = optimset('quadprog'); 

% Call quadprog

[x, fval] = quadprog(H, f, A, b, [], [], [], [], [], options);

% Evaluate results

disp(sprintf('The function value is %f at x = [%f;%f;%f]', fval, x))

Constraint = A*x;

disp(sprintf('Constraint 1 is %f <= %f and constraint 2 is %f <= %f',...

     Constraint(1), b(1), Constraint(2), b(2)))

2) If removing the degenerate variable or variables from the problem is not possible, an alternate workaround is to add a small perturbation to the problem. This perturbation will change the Hessian matrix slightly, potentially enough to cause the Hessian to be nonsingular. Adding the perturbation to the original problem yields:

% Minimize (2+eps)*x(2)^2 + eps*x(1)^2 + eps*x(3)^2 + eps*x(4)^2

H = zeros(4, 4); H(2,2) = 4;

f = zeros(4, 1);

% Subject to -x(3)+x(4)<=-1 and -2*x(2)+x(3)-x(4)<=-1

A = [0 0 -1 1; 0 -2 1 -1];

b = [-1;-1];

% Modify original problem to include perturbation

Perturbation = eps*eye(4);

H = H + Perturbation;

% Create default options structure

options = optimset('quadprog'); 

% Call quadprog

[x, fval] = quadprog(H, f, A, b, [], [], [], [], [], options);

% Evaluate results

disp(sprintf('The function value is %f at x = [%f;%f;%f;%f]', fval, x))

Constraint = A*x;

disp(sprintf('Constraint 1 is %f <= %f and constraint 2 is %f <= %f',...

     Constraint(1), b(1), Constraint(2), b(2)))