Solve minimization problem where constraint is the system of linear inequation with MathOptInterface efficiently

I want to solve

min cx
s.t. Ax + b < 0

with MathOptInterface, where A is a matrix and b and c are vectors.
However, I couldn’t understand the documentation of MathOptInterface.

Currently, my code is as follows:
(Note that A is m * n matrix, b and c are vector with n and m elements respectively.)

using MathOptInterface
const MOI = MathOptInterface
using GLPK

optimizer = GLPK.Optimizer()
x = MOI.add_variables(optimizer, m)
MOI.set(optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0))
MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
for i in 1:n
        MOI.add_constraint(optimizer,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A[i,:],x), b[i]),
                                 MOI.LessThan(0.0))
end
MOI.optimize!(optimizer)

That is, I read the matrix with each rows and create scalar affine constraint.
Actually, this code works, but I think there is more efficient way to solve this problem, in particular using VectorAffineFunction and VectorAffineTerm. I couldn’t get it when I read the docs.

Does anyone knows the solution of this?

This code looks fine to me. Why do you think that it’s inefficient?

Sorry, I said “inefficient” because it takes long time when the matrix A is large. (Is it true?)
Also, I think the code looks sophisticated when I use VectorAffineFunction.
I just want to know how to use VectorAffineFunction.

I just want to know how to use VectorAffineFunction .

You can do it in a way similar to how you build ScalarAffineFunction with broadcasting in your first comment.
Let’s take these A and b for example (any value would work)

julia> A = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> b = collect(1:3)
3-element Array{Int64,1}:
 1
 2
 3

You can do as follows

julia> using MathOptInterface

julia> const MOI = MathOptInterface
MathOptInterface

julia> x = MOI.VariableIndex.(1:3) # should normally comes from `MOI.add_variables`
3-element Array{MathOptInterface.VariableIndex,1}:
 MathOptInterface.VariableIndex(1)
 MathOptInterface.VariableIndex(2)
 MathOptInterface.VariableIndex(3)

julia> terms = MOI.VectorAffineTerm.(1:2, MOI.ScalarAffineTerm.(A, reshape(x, 1, 3)))
2×3 Array{MathOptInterface.VectorAffineTerm{Int64},2}:
 VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(1, VariableIndex(1)))  …  VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(5, VariableIndex(3)))
 VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(2, VariableIndex(1)))     VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(6, VariableIndex(3)))

julia> f = MOI.VectorAffineFunction(vec(terms), b)
MathOptInterface.VectorAffineFunction{Int64}(MathOptInterface.VectorAffineTerm{Int64}[VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(1, VariableIndex(1))), VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(2, VariableIndex(1))), VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(3, VariableIndex(2))), VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(4, VariableIndex(2))), VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(5, VariableIndex(3))), VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(6, VariableIndex(3)))], [1, 2, 3])
1 Like

I got the following error

julia> using MathOptInterface

julia> const MOI = MathOptInterface
MathOptInterface

julia> using GLPK

julia> A = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> b = collect(1:2)
2-element Array{Int64,1}:
 1
 2

julia> c = [1.0, 2.0, 3.0]
3-element Array{Float64,1}:
 1.0
 2.0
 3.0

julia> optimizer = GLPK.Optimizer()
A LinQuadOptInterface model with backend:
Prob(Ptr{Nothing} @0x00007fa352a88810)

julia> x = MOI.add_variables(optimizer, 3)
3-element Array{MathOptInterface.VariableIndex,1}:
 MathOptInterface.VariableIndex(1)
 MathOptInterface.VariableIndex(2)
 MathOptInterface.VariableIndex(3)

julia> MOI.set(optimizer,
               MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
               MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(c, x), 0.0))

julia> MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
MIN_SENSE::OptimizationSense = 0

julia> terms = MOI.VectorAffineTerm.(1:2, MOI.ScalarAffineTerm.(A, reshape(x, 1, 3)))
2×3 Array{MathOptInterface.VectorAffineTerm{Int64},2}:
 VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(1, VariableIndex(1)))  …  VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(5, VariableIndex(3)))
 VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(2, VariableIndex(1)))     VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(6, VariableIndex(3)))

julia> f = MOI.VectorAffineFunction(vec(terms), b)
MathOptInterface.VectorAffineFunction{Int64}(MathOptInterface.VectorAffineTerm{Int64}[VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(1, VariableIndex(1))), VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(2, VariableIndex(1))), VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(3, VariableIndex(2))), VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(4, VariableIndex(2))), VectorAffineTerm{Int64}(1, ScalarAffineTerm{Int64}(5, VariableIndex(3))), VectorAffineTerm{Int64}(2, ScalarAffineTerm{Int64}(6, VariableIndex(3)))], [1, 2])

julia> MOI.add_constraint(optimizer, f, MOI.Nonpositives(2))
ERROR: MathOptInterface.UnsupportedConstraint{MathOptInterface.VectorAffineFunction{Int64},MathOptInterface.Nonpositives}: `MathOptInterface.VectorAffineFunction{Int64}`-in-`MathOptInterface.Nonpositives` constraints is not supported by the model.
Stacktrace:
 [1] #correct_throw_add_constraint_error_fallback#30(::MathOptInterface.AddConstraintNotAllowed{MathOptInterface.VectorAffineFunction{Int64},MathOptInterface.Nonpositives}, ::Function, ::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Int64}, ::MathOptInterface.Nonpositives) at /Users/myname/.julia/packages/MathOptInterface/C3lip/src/constraints.jl:122
 [2] correct_throw_add_constraint_error_fallback(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Int64}, ::MathOptInterface.Nonpositives) at /Users/myname/.julia/packages/MathOptInterface/C3lip/src/constraints.jl:119
 [3] #throw_add_constraint_error_fallback#27(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Int64}, ::MathOptInterface.Nonpositives) at /Users/myname/.julia/packages/MathOptInterface/C3lip/src/constraints.jl:99
 [4] throw_add_constraint_error_fallback(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Int64}, ::MathOptInterface.Nonpositives) at /Users/myname/.julia/packages/MathOptInterface/C3lip/src/constraints.jl:99
 [5] add_constraint(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Int64}, ::MathOptInterface.Nonpositives) at /Users/myname/.julia/packages/MathOptInterface/C3lip/src/constraints.jl:83
 [6] top-level scope at none:0

How can I fix this error? Why VectorAffineFunction doesn’t support Nonpositives? Linear constraints says it supports Nonpositives.

Use Float64 instead of Int coefficients in A and b. This is slight usability bug from the MathOptInterface perspective, but isn’t visible from JuMP since it automatically converts values to Float64s.

Thanks, my code works well with using Float64 !

The example script is really helpful to understand how to use VectorAffineTerm and Function.

One thing that I cannot understand is: why do we need to apply reshape(x, 1, n) in ScalarAffineTerm. ? I can hardly couple the codes and their corresponding constraints.

Also, when matrices and vectors are given for an optimization problem (e.g. A, c, lb, ub, etc.), is there a simple way to solve it with MathOptInterface as we do so with MathProgBase?

One thing that I cannot understand is: why do we need to apply reshape(x, 1, n) in ScalarAffineTerm. ? I can hardly couple the codes and their corresponding constraints.

You can to create a matrix of terms. The term (i, j) corresponds to the entry (i, j) of the matrix so it should have row i and variable x[j].
In the broadcast, 1:2 and x are one dimensional and A is 2-dimensional so they will be broadcasted constant along the missing dimension.
As 1:2 is constant in the second dimension, it is broadcasted correctly but x should be broadcasted constant in the first dimension so we reshape it into a row vector for that.

Also, when matrices and vectors are given for an optimization problem (e.g. A , c , lb , ub , etc.), is there a simple way to solve it with MathOptInterface as we do so with MathProgBase ?

The simple way is to do it in JuMP with

@variable(model, lb[i] <= x[i=1:n] <= ub[i])
@constraint(model, l .<= A * x .<= u)

For a way similar to MathProgBase, see a relevant discussion here:

Hi, This is still giving me errors. I’m running the following code:

using MathOptInterface
const MOI = MathOptInterface
using GLPK

optimizer = GLPK.Optimizer()

A = reshape(1:6, 2, 3)
A = convert(AbstractArray{Float64}, A)
b = collect(1:3)
b = convert(AbstractArray{Float64}, b)

x = MOI.VariableIndex.(1:3) # should normally comes from `MOI.add_variables`

terms = MOI.VectorAffineTerm.(1:2, MOI.ScalarAffineTerm.(A, reshape(x, 1, 3)))

f = MOI.VectorAffineFunction(vec(terms), b)

MOI.add_constraint(optimizer, f, MOI.Nonpositives(2))

and getting the following error:

ERROR: MathOptInterface.UnsupportedConstraint{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Nonpositives}: `MathOptInterface.VectorAffineFunction{Float64}`-in-`MathOptInterface.Nonpositives` constraint is not supported by the model.
Stacktrace:
 [1] correct_throw_add_constraint_error_fallback(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.Nonpositives; error_if_supported::MathOptInterface.AddConstraintNotAllowed{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Nonpositives}) at C:\Users\pnvol\.julia\packages\MathOptInterface\VjkSQ\src\constraints.jl:158
 [2] correct_throw_add_constraint_error_fallback(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.Nonpositives) at C:\Users\pnvol\.julia\packages\MathOptInterface\VjkSQ\src\constraints.jl:155
 [3] throw_add_constraint_error_fallback(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.Nonpositives; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\pnvol\.julia\packages\MathOptInterface\VjkSQ\src\constraints.jl:135
 [4] throw_add_constraint_error_fallback(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.Nonpositives) at C:\Users\pnvol\.julia\packages\MathOptInterface\VjkSQ\src\constraints.jl:135
 [5] add_constraint(::GLPK.Optimizer, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.Nonpositives) at C:\Users\pnvol\.julia\packages\MathOptInterface\VjkSQ\src\constraints.jl:119
 [6] top-level scope at C:\Users\pnvol\Dropbox\JYC\Code\J15_2017_chart\MWE_MOI.jl:90

Any ideas? I just really miss linprog.

GLPK doesn’t support VectorAffineFunction

Either use

optimizer = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)

or pass ScalarAffineFunctions.

I just really miss linprog.

Use JuMP:

model = Model(GLPK.Optimizer)
@variable(model, l[i] <= x[i = 1:n] <= u[i])
@objective(model, c' * x)
@constraint(model, A * x .>= b)
optimize!(model)

Thanks. Both those solutions work. Much appreciated.