Numerical issues solving a simple QCLP model
AnsweredI 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
:

Hi Jack,
Suboptimal 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:
 set BarHomogeneous parameter to 1
 set NumericFocus to 2 or 3
 experiment with the value of BarCorrectors, e.g., set it to 10 or 1000
 turn off Presolve
 decrease the tolerances BarConvTol, BarQCPConvTol, FeasibilityTol
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 thirdparty 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ł 
Hi Jaromił,
Thanks for the help. I still get some suboptimal termination warnings when playing with the various tuning parameters in both the original and reformulated 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.9999999999994527e01
C0 C1 5.0000000000000033e01
C1 C0 5.0000000000000033e01
C1 C1 5.0000000000000033e01
ENDATAMany thanks!
Jack

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.5These 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 suboptimal 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 suboptimal solution. This confirms the influence of the roundoff 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ł
Please sign in to leave a comment.
Comments
3 comments