Incorrect reading of MPS files with QCMATRIX field

[sebastiendesignolle]

When reading MPS files with quadratic constraints, off-diagonal elements seem to get a factor two. Here is a minimal working example:

NAME          mwe
ROWS
 N  obj
 E  c1
 G  q1
COLUMNS
    x1        obj       1.0
    x1        c1        1.0
    x2        c1        -1.0
RHS
    rhs       c1        0.0
    rhs       q1        1.0
QCMATRIX      q1
    x1        x1        1.0
    x1        x2        1.0
    x2        x1        1.0
ENDATA

Solving this problem with Gurobi within Julia (reading the file with MOI) leads to a solution of 1/√3, while solving it directly with Gurobi gives the expected solution of 1/√2.

[odow]

Ugh. I just cannot get this correct. I'll fix.

[odow]

What version of MOI are you using?

This looks correct to me:

julia> import MathOptInterface as MOI

julia> begin
           data = """
           NAME          mwe
           ROWS
           N  obj
           E  c1
           G  q1
           COLUMNS
               x1        obj       1.0
               x1        c1        1.0
               x2        c1        -1.0
           RHS
               rhs       c1        0.0
               rhs       q1        1.0
           QCMATRIX      q1
               x1        x1        1.0
               x1        x2        1.0
               x2        x1        1.0
           ENDATA
           """
           io = IOBuffer(data)
           model = MOI.FileFormats.MPS.Model()
           seekstart(io)
           read!(io, model)
           m = MOI.Utilities.Model{Float64}()
           MOI.copy_to(m, model)
           print(m)
       end
Minimize ScalarAffineFunction{Float64}:
 -0.0 + 1.0 x1

Subject to:

ScalarAffineFunction{Float64}-in-EqualTo{Float64}
 0.0 + 1.0 x1 - 1.0 x2 == 0.0  (c1)

ScalarQuadraticFunction{Float64}-in-GreaterThan{Float64}
 0.0 + 1.0 x1² + 2.0 x1*x2 >= 1.0  (q1)

VariableIndex-in-GreaterThan{Float64}
 x1 >= 0.0
 x2 >= 0.0

[odow]

I'm going to guess that this was fixed by #2628

[odow]

How did you get the sqrt(2) solution?

$ cat model.mps 
NAME          mwe
ROWS
 N  obj
 E  c1
 G  q1
COLUMNS
    x1        obj       1.0
    x1        c1        1.0
    x2        c1        -1.0
RHS
    rhs       c1        0.0
    rhs       q1        1.0
QCMATRIX      q1
    x1        x1        1.0
    x1        x2        1.0
    x2        x1        1.0
ENDATA

$ gurobi_cl model.mps
Set parameter LicenseID to value 890341
Set parameter LogFile to value "gurobi.log"
Using license file /Users/oscar/gurobi.lic

Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (mac64[x86] - Darwin 24.1.0 24B83)
Copyright (c) 2024, Gurobi Optimization, LLC

Read MPS format model from file model.mps
Reading time = 0.00 seconds
mwe: 1 rows, 2 columns, 2 nonzeros

CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 1 rows, 2 columns and 2 nonzeros
Model fingerprint: 0x6e9e8535
Model has 1 quadratic constraint
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [1e+00, 2e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [0e+00, 0e+00]
  QRHS range       [1e+00, 1e+00]
Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored
Presolve time: 0.00s
Presolved: 4 rows, 5 columns, 8 nonzeros
Presolved model has 1 second-order cone constraint
Ordering time: 0.00s

Barrier statistics:
 AA' NZ     : 5.000e+00
 Factor NZ  : 1.000e+01
 Factor Ops : 3.000e+01 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   2.01798596e+00  0.00000000e+00  2.02e+00 0.00e+00  3.25e-01     0s
   1   8.38756216e-01  1.92176195e-01  6.88e-01 2.78e-17  8.74e-02     0s
   2   6.59824808e-01  3.63072890e-01  7.57e-07 6.59e-17  5.94e-02     0s
   3   6.23053714e-01  5.67578847e-01  3.20e-10 1.11e-16  1.11e-02     0s
   4   5.79265651e-01  5.76947619e-01  8.36e-13 1.11e-16  4.64e-04     0s
   5   5.77413011e-01  5.77323204e-01  1.55e-12 2.22e-16  1.80e-05     0s
   6   5.77350495e-01  5.77349299e-01  5.24e-12 3.55e-15  2.39e-07     0s

Barrier solved model in 6 iterations and 0.00 seconds (0.00 work units)
Optimal objective 5.77350495e-01

[odow]

And getting Gurobi to write back out its view of the model shows the same thing:

(base) oscar@MacBookPro /tmp % cat model.mps 
NAME          mwe
ROWS
 N  obj
 E  c1
 G  q1
COLUMNS
    x1        obj       1.0
    x1        c1        1.0
    x2        c1        -1.0
RHS
    rhs       c1        0.0
    rhs       q1        1.0
QCMATRIX      q1
    x1        x1        1.0
    x1        x2        1.0
    x2        x1        1.0
ENDATA
(base) oscar@MacBookPro /tmp % gurobi_cl ResultFile=model.lp model.mps 
Set parameter LicenseID to value 890341
Set parameter LogFile to value "gurobi.log"
Using license file /Users/oscar/gurobi.lic

Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (mac64[x86] - Darwin 24.1.0 24B83)
Copyright (c) 2024, Gurobi Optimization, LLC

Read MPS format model from file model.mps
Reading time = 0.00 seconds
mwe: 1 rows, 2 columns, 2 nonzeros

CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 1 rows, 2 columns and 2 nonzeros
Model fingerprint: 0x6e9e8535
Model has 1 quadratic constraint
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [1e+00, 2e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [0e+00, 0e+00]
  QRHS range       [1e+00, 1e+00]
Warning for adding constraints: zero or small (< 1e-13) coefficients, ignored
Presolve time: 0.00s
Presolved: 4 rows, 5 columns, 8 nonzeros
Presolved model has 1 second-order cone constraint
Ordering time: 0.00s

Barrier statistics:
 AA' NZ     : 5.000e+00
 Factor NZ  : 1.000e+01
 Factor Ops : 3.000e+01 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   2.01798596e+00  0.00000000e+00  2.02e+00 0.00e+00  3.25e-01     0s
   1   8.38756216e-01  1.92176195e-01  6.88e-01 2.78e-17  8.74e-02     0s
   2   6.59824808e-01  3.63072890e-01  7.57e-07 6.59e-17  5.94e-02     0s
   3   6.23053714e-01  5.67578847e-01  3.20e-10 1.11e-16  1.11e-02     0s
   4   5.79265651e-01  5.76947619e-01  8.36e-13 1.11e-16  4.64e-04     0s
   5   5.77413011e-01  5.77323204e-01  1.55e-12 2.22e-16  1.80e-05     0s
   6   5.77350495e-01  5.77349299e-01  5.24e-12 3.55e-15  2.39e-07     0s

Barrier solved model in 6 iterations and 0.00 seconds (0.00 work units)
Optimal objective 5.77350495e-01


Wrote result file 'model.lp'

(base) oscar@MacBookPro /tmp % cat model.lp
\ Model mwe
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
  x1
Subject To
 c1: x1 - x2 = 0
 q1: [ x1 ^2 + 2 x1 * x2 ] >= 1
Bounds
End

[odow]

Part of the confusion is that Gurobi's QCMATRIX does NOT have an associated /2. I don't know why they chose that convention...

[sebastiendesignolle]

Thanks for the quick answer! I'm using MOI v1.37.0, but it could be that the problem comes from Gurobi, as the 1/√2 was obtained with Gurobi 10, but your log for Gurobi 11 and my test with Gurobi 12 both give 1/√3. Here's the Gurobi 10 log:

Gurobi Optimizer version 10.0.1 build v10.0.1rc0 (linux64)
Copyright (c) 2023, Gurobi Optimization, LLC

Read MPS format model from file mwe.mps
Reading time = 0.00 seconds
mwe: 1 rows, 2 columns, 2 nonzeros

Optimize a model with 1 rows, 2 columns and 2 nonzeros
Model fingerprint: 0x6e9e8535
Model has 1 quadratic constraint
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  QMatrix range    [1e+00, 2e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [0e+00, 0e+00]
  RHS range        [0e+00, 0e+00]
  QRHS range       [1e+00, 1e+00]
Presolve time: 0.01s
Presolved: 4 rows, 5 columns, 9 nonzeros
Presolved model has 1 second-order cone constraint
Ordering time: 0.00s

Barrier statistics:
 AA' NZ     : 5.000e+00
 Factor NZ  : 1.000e+01
 Factor Ops : 3.000e+01 (less than 1 second per iteration)
 Threads    : 1

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   2.04290743e+00  0.00000000e+00  6.15e-01 0.00e+00  4.73e-01     0s
   1   8.74170935e-01  2.45674909e-01  1.86e-01 2.78e-17  1.13e-01     0s
   2   7.36994006e-01  6.67146028e-01  2.05e-07 1.11e-16  1.40e-02     0s
   3   7.08894130e-01  7.07019104e-01  3.42e-13 1.11e-16  3.75e-04     0s
   4   7.07108725e-01  7.07106684e-01  1.08e-10 4.88e-15  4.08e-07     0s
   5   7.07106783e-01  7.07106781e-01  2.82e-08 1.84e-12  4.68e-10     0s

Barrier solved model in 5 iterations and 0.01 seconds (0.00 work units)
Optimal objective 7.07106783e-01

I'm checking with the actual problem the issue comes from and keep you tuned. If it was a false alarm, I'll apologise and close the issue.

[odow]

The solution from Gurobi 10 is incorrect.

If the model is

Minimize
  x1
Subject To
 c1: x1 - x2 = 0
 q1: [ x1 ^2 + 2 x1 * x2 ] >= 1

we can project out c1 to get

Minimize
  x1
Subject To
 q1: [ x1 ^2 + 2 x1 * x1 ] >= 1

which simplifies to

Minimize
  x1
Subject To
 q1: [ 3 x1 ^2 ] >= 1

and then to

Minimize
  x1
Subject To
 q1: x1 >= 1 / sqrt(3)

Closing this issue because everything seems to be working as expected on our end.