Issue with (nonconvex) quadratic constraints

[tkoolen]

Consider the following:

using JuMP
using SCIP
using LinearAlgebra
m = Model(with_optimizer(SCIP.Optimizer))
@variable m x[1 : 3]
@variable m y[1 : 3]
@variable m z[1 : 3]
w = x × y - x × z
@constraint m w .== 1
optimize!(m)

SCIP reports that the problem is solved, but value.(w) results in

3-element Array{Float64,1}:
 1.9999999962747097
 1.9999999925494194
 1.9999999962747097

i.e., twice what it should be.

[tkoolen]

I think the issue is

https://github.com/SCIP-Interfaces/SCIP.jl/blob/5507c001e7988d95972692f07d2f511c94917224/src/MOI_wrapper/quadratic_constraints.jl#L20

Note the difference between JuMP.QuadExpr and MOI.ScalarQuadraticFunction, https://github.com/JuliaOpt/JuMP.jl/blob/c01802926c34aa30e5d272ac1063f1595084aa4a/src/quad_expr.jl#L336-L349. MOI.ScalarQuadraticFunction already has the 1/2 factor for terms involving distinct variables, as shown in the (confusing) documentation.

[rschwarz]

Thanks for reporting. It's quite possible that there is a problem with the transformation.

However, I read the documentation differently. Actually, I first took the coefficients as they are, and got failing MOI tests, to then look up the documentation and add the scaling factor. Indeed, MOI's qcp tests are passing, in particular qcp2, where we have

#  x^2 <= 2 (c)
cf = MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(0.0, x)], [MOI.ScalarQuadraticTerm(2.0, x, x)], 0.0)
c = MOI.add_constraint(model, cf, MOI.LessThan(2.0))

with the the factor 2.0 for x*x which I'm correcting for SCIP.

Unfortunately, your example model is not completely obvious to me (and I get an ERROR: UndefVarError: × not defined in the Julia REPL), but I will try again with a simpler model (from JuMP), such as x^2 - 4 == 0.

[rschwarz]

OK, so I find that the JuMP to MOI conversion makes a distinction between quadratic (x^2) and bilinear (x*y) terms. For

using JuMP
using SCIP

model = Model(with_optimizer(SCIP.Optimizer))
@variable(model, x >= 0)
@variable(model, y >= 0)
@constraint(model, x * x == 4)
@constraint(model, x * y == 2)

optimize!(model)
@show value(x) value(y)

I get

value(x) = 2.0
value(y) = 1.999999999

So, the term x^2 is handled properly, while x*y is scaled incorrectly.

I guess the section of the MOI doc that says:

Duplicate indices in a or Q are accepted, and the corresponding coefficients are summed together. "Mirrored" indices (q,r) and (r,q) (where r and q are VariableIndexes) are considered duplicates; only one need be specified.

means that mirrored indices are expected, and otherwise should have a factor of 2 in the coefficient.

[tkoolen]

Exactly, I was typing up just that.

Just because I was almost done with it, here's my version:

(and I get an ERROR: UndefVarError: × not defined in the Julia REPL)

× (cross, the cross product) is in LinearAlgebra. Maybe something's going wrong with copy-pasting unicode. Anyway, this simpler version already demonstrates the problem:

using JuMP
using SCIP

m = Model(with_optimizer(SCIP.Optimizer))
@variable m x
@variable m y
w = x * y
@constraint m w == 1
optimize!(m)
@show value(w)

The issue is that the case of distinct variables (off-diagonal terms in the Q matrix of the MOI.ScalarAffineFunction) needs to be treated differently than the case of non-distinct variables (diagonal terms), as witnessed by https://github.com/JuliaOpt/JuMP.jl/blob/c01802926c34aa30e5d272ac1063f1595084aa4a/src/quad_expr.jl#L343-L345. This is because MOI.ScalarAffineTerm. This is because MOI.ScalarQuadraticFunction uses its list ofScalarQuadraticTerms to define a symmetric matrix in sparse form:

• a term with distinct variables, like x * y, would have two entries in Q, corresponding to 1/2 * x * y and 1/2 * y * x, and since 1) ""Mirrored" indices (q,r) and (r,q) [...] are considered duplicates; only one needs to be specified." and 2) "Duplicate indices in a or Q are accepted, and the corresponding coefficients are summed together, you only need the single ScalarQuadraticTerm(1/2, x, y) in the ScalarQuadraticFunction
• a term like x * x will have only a single entry in Q on the diagonal, and so the "mirrored indices" rule doesn't apply.

[rschwarz]

× (cross, the cross product) is in LinearAlgebra.

Yes, I saw that using later. I'll produce a patch for this issue shortly...

[rschwarz]

I tagged a new release with the fix: https://github.com/JuliaLang/METADATA.jl/pull/22146

[tkoolen]

Great, and thanks very much for your hard work on the wrapper!