Skip to main content

Numerical issues solving a simple QCLP model

Answered
Jack Light
Gurobi-versary
First Comment
First Question

I am trying to solve a simple QCLP model with Gurobi but keep running into numerical issues. Any help much appreciated!

A MWE is as follows:

- Consider the regression of Y on a constant and a variable X:

Y = b_1 + b_2*X + residuals

 - I use Gurobi to find lower and upper bounds of b_2 subject to the constraint that the corresponding sum of squared residuals be weakly less than the minimized sum of squared residuals.

- The bounds in this example should collapse to a point whose value is equivalent to the coefficient on X obtained from a standard regression of Y on X and a constant.

- I am using AMPL with Gurobi as the solver, but am seeing several numerical issues with the barrier algorithm used. In particular the bounds fail to converge to a point, but to a small bound either side of the point they should converge too. The issue persists even in very small models (e.g. with  N <10) despite playing with various Gurobi numerical parameters.

- E.g.: If OLS recovers a value of b_2 = 1.587, Gurobi returns the bound [1.585 1.588]

- MWE code as follows (written for AMPL):

param N, default 100;
param X1 {j in 1..N};
param Y {j in 1..N};
var beta {k in 1..3} >=-10, <=10 ;

minimize minPar:
beta[2];

maximize maxPar:
beta[2];

subject to SS_Obj: (1/N)*(sum {j in 1..N} ( Y[j] - beta[1] - X1[j]*beta[2] )^2) <= Obj_value;

Where Obj_value is obtained as the minimized value of:

(1/N)*(sum {j in 1..N} ( Y[j] - beta[1] - X1[j]*beta[2])^2);

The log file output is as follows

:

 

 

0

Comments

3 comments

Date Votes
  • Jaromił Najman
    Gurobi Staff Gurobi Staff
    • Edited

    Hi Jack,

    Sub-optimal termination is printed when the barrier algorithm does not make any (or only very small) progress in the solution value and cannot decrease the dual residual, meaning that although the final convergence criterion is met, the dual infeasibility cannot be reduced. This is often the case for shaky numerics. 

    You can try to overcome this issue by either reformulating your problem and/or experimenting with the following parameters; this may however result in a performance loss:

    You could also try reformulating your problem. It currently reads
    \[\min_{b_1,b_2} \sum_j (y_j - b_1 - x_j\cdot b_2)^2 \]

    where \(y_j, x_j\) are constants and \(b_{1},b_2\) are variables.
    You could  formulate it as

    \[\min_{b_1,b_2} \sum_j (y_j^2 - 2\cdot y_j\cdot b_1 - 2\cdot y_j\cdot x_j \cdot b_2) + \sum_j (b_1 + x_j\cdot b_2)^2\]

    If none of the above helps, could you then provide an MPS or an LP file of the model, see How do I write an MPS model file from a third-party API and How do I write an "MPS file" for my problem? To make this file available to the community, please have a look at Posting to the Community Forum.

    Best regards,
    Jaromił

    0
  • Jack Light
    • Edited
    Gurobi-versary
    First Comment
    First Question

    Hi Jaromił,

    Thanks for the help. I still get some sub-optimal termination warnings when playing with the various tuning parameters in both the original and re-formulated problem. Below is the MPS output for the following problem (the upper bound on b_2 is found analogously):

     

    \begin{align}
    \underset{b_1,b_2}{\min} \ b_2 &\quad s.t. \\
    \sum_{j=1}^N (y_j - b_1 -b_2x_j)^2 &\leq \underset{b_1,b_2}{\min} \sum_{j=1}^N (y_j - b_1 -b_2x_j)^2 \\
    b_1,b_2 &\in [-5,5]
    \end{align}

     

    MPS:

    NAME foo
    OBJSENSE MAX
    ROWS
     N  OBJ
     L  qc0     
    COLUMNS
        C0        OBJ       0
        C0        qc0       -5.0347716616104652e+00
        C1        OBJ       1
        C1        qc0       -4.0309243407684274e+00
    RHS
        RHS1      qc0       -8.6280302422805661e+00
    BOUNDS
     LO BND1      C0        -5
     UP BND1      C0        5
     LO BND1      C1        -5
     UP BND1      C1        5
    QCMATRIX   qc0     
        C0        C0        9.9999999999994527e-01
        C0        C1        5.0000000000000033e-01
        C1        C0        5.0000000000000033e-01
        C1        C1        5.0000000000000033e-01
    ENDATA

    Many thanks!

    Jack

    0
  • Jaromił Najman
    Gurobi Staff Gurobi Staff
    • Edited

    Hi Jack,

    The coefficients of your QCMatrix suffer from round off errors. The values should most likely read

    QCMATRIX   qc0     
        C0        C0        1
        C0        C1        0.5
        C1        C0        0.5
        C1        C1        0.5

    These small round off errors may cause a lot of trouble during the solution process as they are far below the feasibility tolerance of the solver but still may have influence on the solution process. Moreover, it seems like the digits way below feasibility tolerance of the linear coefficients have significant meaning to the solution of the problem. This explains the sub-optimal solution when minimizing the modified problem. You can tackle this with setting NumericFocus=3 and BarCorrectors=100 to achieve additional solution stability.

    If I write the MPS file to an LP file, which holds only limited precision of the coefficients, the barrier algorithm does not run in a sub-optimal solution. This confirms the influence of the round-off errors and digits below the feasibility tolerance in the coefficients.

    For more information on numerics, I recommend the excellent article by Ed Klotz on Numerical Instability in Linear and Integer Programs. The article focuses on linear programs, but the main ideas also apply to quadratic programs.

    Best regards,
    Jaromił

    0

Please sign in to leave a comment.

Didn't find what you were looking for?

New post