CONEPROG: Precision Issue (Exitflag -7)

These days I tried to port a symmetric positive semi-definite quadratic programming problem to a conic optimization problem. So-far I was able to transcribe the problem into a cone program, however it seems that still some options tuning is necessary, since all the iterations end up with an exitflag -7. Before I shall present the problem specification as a conic optimization problem and a snippet of the evaluation process, I shortly present a result of an example that is close to the solution. Moreover, it is worth noting that some examples finished up with the correct solution despite the consecutive appearance of the exitflag -7. In addition, I followed the piece of advice of the documentation to change the solver, but this had no effect. The same happened while changing the tolerance value from up to .

Due to the following returned outcome of an example, I suppose that I am on the right track, though the result below is not correct

x =
    8.4966   11.5004    5.5087   23.0781   50.0829

The returned result of the evaluation is close to the expected solution below

pk_bv =

    8.5000   11.5000    5.5000   23.0833   50.0833    

which let me suppose that I encounter some numerical issues that I cannot localize. In the hope that a different angle may solve the puzzle, I present now the parameters of the problem. Moreover, we have , it should hold with equality.

The subsequent options setting has been chosen:

 opts = optimoptions("coneprog","OptimalityTolerance",1e-12,"ConstraintTolerance",1e-12,"LinearSolver","auto")

And the cone problem is specified as:

    Q=E'*E; % Matrix Q is symmetric positive semi-definite. Matrix E is made up by the numbers -1,0,1, it is not square.
    b=-E'*a;
    Aineq=[];
    bineq=[];
    c=[b(:);1];  % objective function
    Asc = sqrtm(Q);
    Asc((end+1),(end+1)) = 1;
    d = [zeros(size(b(:)));1];
    bsc = zeros(size(d));
    gamma = -1;
    socConstraints = secondordercone(Asc,bsc,d,gamma);
    Aeq=[E(m,:),0];
    lb=[vi(:);-inf];
    ub=[ra(:);inf];
    [u,fval,exitflag,output,lambda]=coneprog(c,socConstraints,Aineq,bineq,Aeq,a(m),lb,ub,opts)

In this respect, it might be useful to compare the above problem with the symmetric positive semi-definite quadratic programming problem:

 [x,fval,exitflag,output,lambda] = quadprog(Q,b,[],[],E(m,:),a(m),vi,ra,x,opts);

Furthermore, I want now to present some snippets of the evaluation process to observe the options and output structure settings.

 >> x=cp_kernel(bv)

opts = 

  coneprog options:

   Set properties:
    ConstraintTolerance: 1.0000e-12
           LinearSolver: 'auto'
    OptimalityTolerance: 1.0000e-12

   Default properties:
                Display: 'final'
          MaxIterations: 200
                MaxTime: Inf

Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.

u =

   1.0e+04 *

    0.0008
    0.0012
    0.0005
    0.0023
    0.0050
    5.2886


fval =

  -5.2785e+04


exitflag =

    -7


output = 

  struct with fields:

           iterations: 30
    primalfeasibility: 1.0206e-11
      dualfeasibility: 6.6721e-07
           dualitygap: 9.6139e-12
            algorithm: 'interior-point'
         linearsolver: 'augmented'
              message: 'Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.'


lambda = 

  struct with fields:

      eqlin: 94.9844
    ineqlin: [0x1 double]
        soc: 9.9244e-09
      lower: [6x1 double]
      upper: [6x1 double]


ans =

    8.4783   11.8906    5.2411   23.1990   49.8576


         .
	 .
	 .
	 .
	 .
  exitflag =

    -7


output = 

  struct with fields:

           iterations: 28
    primalfeasibility: 1.0126e-11
      dualfeasibility: 2.3589e-09
           dualitygap: 1.0380e-11
            algorithm: 'interior-point'
         linearsolver: 'augmented'
              message: 'Search direction is small and current iterate is not within the specified constraint and/or optimality tolerance.'


lambda = 

  struct with fields:

      eqlin: -6.4447
    ineqlin: [0x1 double]
        soc: 2.7285e-08
      lower: [6x1 double]
      upper: [6x1 double]

Warning: Probably no kernel point found! 
> In cp_kernel>computePrk (line 215)
In cp_kernel (line 102) 

x =

    8.4966   11.5004    5.5087   23.0781   50.0829   

My question is now: which options must be tuned to get the correct result? Or has the cone problem specification a major flaw, which I have overlooked?

Thank you very much in advance for the kind support and any suggestion of improvements.

EDIT 15 Sep 2021

Thanks to the kind and excellent support of Aykut Bulut I could solve the above precision issue by reformulating the cone problem. For those who are looking for some inspiration for their own work, I present the reformulated cone problem by

    Aineq=[];
    bineq=[];
    c=[zeros(size(b(:)));1];  % objective function
    Asc = E;
    Asc((end+1),(end+1)) = 1;
    d = [zeros(size(b(:)));1];
    bsc=a';
    bsc(end+1)=0;
    gamma = 0;
    socConstraints = secondordercone(Asc,bsc,d,gamma);
    Aeq=[E(m,:),0];
    lb=[-inf(n+1,1)];
    ub=[inf(n+1,1)];
    [u,fval,exitflag,output,lambda]=coneprog(c,socConstraints,Aineq,bineq,Aeq,a(m),lb,ub,opts);

Answer 1

[Aykut Bulut]

Hi Holger,

The method you use for conversion is usefull when you do not have factors of Q. In your case, you can skip computing Q and using sqrtm() function to factor it. I think some precision might be lost in these operations and skipping them might improve the numerics.

Objective function is

First observe that we can drop the constant term in the objective function. It can be added back once the optimal solution is computed. Second, instead of minimizing , we can minimize (note that we dropped the exponent 2 and \frac{1}{2} coefficient, these problems have the same set of optimal solutions and effectively same). Next, instead of minimizing , we define a new variable t and force it to be greater than equal to in the constraints, and minimize t. After this, objective function becomes t, and we have the following additional constraint.

I think solving this modified problem (together with other constraints) will yield more precise solutions. This new problem will have a different optimal objective value. However, once you have optimal x for the new problem, you can compute the optimal objective function value of the original problem by using optimal x values computed.

[Holger I. Meinhardt]

Thank you very much for your kind and excellent support, it works now. Your proposal was accompanied with my expectation. However, I supposed that the function secondordercone() can only handle square matrices, this information I had overlooked in the documentation. All examples I studied from the documentation consider square matrices, so I assumed this was the only possible way to handle my problem. Anyway, I get now the correct solutions, though the exitflag -7 is still persistent, but it has no repercussion on the accuracy of the final results, and in most cases I get on the final stage an exitflag of 1.

[Aykut Bulut]

Thanks for your feedback. I will inform the team about the examples in the documentation.