Skip to content

Support QP interface via MOI #30

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mtanneau wants to merge 10 commits into
base: main
Choose a base branch
from
Open

Support QP interface via MOI #30

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
a926678
Support QP interface via MOI
mtanneau Jan 10, 2026
bbcb4f4
Move objective data extraction to separate method
mtanneau Jan 10, 2026
ace38fd
Format
mtanneau Jan 10, 2026
a7ff4a2
Fix typos
mtanneau Jan 10, 2026
7c6a48a
In-place offset + remove comments
mtanneau Jan 11, 2026
067d973
Move import
mtanneau Jan 11, 2026
1841980
Bump cuopt compat version to 25.12
mtanneau Jan 12, 2026
c108151
Explain MOI vs cuOpt convention re:quadratic functions
mtanneau Jan 12, 2026
b27e7d0
Add QP unit tests
mtanneau Jan 13, 2026
2ee820d
Fix typos + format
mtanneau Jan 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,12 @@ version = "0.1.2"
MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[compat]
MathOptInterface = "1.34"
PrecompileTools = "1"
SparseArrays = "1"
Test = "1.6"
julia = "1.6"

Expand Down
174 changes: 148 additions & 26 deletions src/MOI_wrapper.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@
# The HiGHS wrapper is released under an MIT license, a copy of which can be
# found in `/thirdparty/THIRD_PARTY_LICENSES` or at https://opensource.org/licenses/MIT.

import SparseArrays

import MathOptInterface as MOI
const CleverDicts = MOI.Utilities.CleverDicts

Expand Down Expand Up @@ -621,8 +623,13 @@ end

function MOI.supports(
::Optimizer,
::Union{MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}},
)
::MOI.ObjectiveFunction{F},
) where {
F<:Union{
MOI.ScalarAffineFunction{Float64},
MOI.ScalarQuadraticFunction{Float64},
},
}
return true
end

Expand Down Expand Up @@ -651,6 +658,7 @@ function _check_input_data(dest::Optimizer, src::MOI.ModelLike)
if attr in (
MOI.Name(),
MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}(),
MOI.ObjectiveSense(),
)
continue
Expand Down Expand Up @@ -903,17 +911,95 @@ function _get_objective_data(

objective_sense = sense == MOI.MIN_SENSE ? CUOPT_MINIMIZE : CUOPT_MAXIMIZE

objective_coefficients = zeros(Float64, numcol)
F = MOI.get(src, MOI.ObjectiveFunctionType())
f_obj = MOI.get(src, MOI.ObjectiveFunction{F}())

objective_offset,
objective_coefficients_linear,
qobj_row_offsets,
qobj_col_indices,
qobj_matrix_values = _get_objective_data(f_obj, mapping, numcol)

return objective_sense,
objective_offset,
objective_coefficients_linear,
qobj_row_offsets,
qobj_col_indices,
qobj_matrix_values
end

function _get_objective_data(
f_obj::MOI.ScalarAffineFunction,
mapping,
numcol::Int32,
)
objective_offset = f_obj.constant

objective_coefficients_linear = zeros(Float64, numcol)
for term in f_obj.terms
objective_coefficients[mapping[term.variable].value] += term.coefficient
i = mapping[term.variable].value
objective_coefficients_linear[i] += term.coefficient
end

return objective_offset,
objective_coefficients_linear,
Int32[0],
Int32[],
Float64[]
end

function _get_objective_data(
f_obj::MOI.ScalarQuadraticFunction,
mapping,
numcol::Int32,
)
objective_offset = f_obj.constant

return objective_sense, objective_offset, objective_coefficients
objective_coefficients_linear = zeros(Float64, numcol)
for term in f_obj.affine_terms
i = mapping[term.variable].value
objective_coefficients_linear[i] += term.coefficient
end

# Extract quadratic objective
Qtrows = Int32[]
Qtcols = Int32[]
Qtvals = Float64[]
sizehint!(Qtrows, length(f_obj.quadratic_terms))
sizehint!(Qtcols, length(f_obj.quadratic_terms))
sizehint!(Qtvals, length(f_obj.quadratic_terms))
for qterm in f_obj.quadratic_terms
i = mapping[qterm.variable_1].value
j = mapping[qterm.variable_2].value
v = qterm.coefficient
# MOI stores quadratic functions as `¹/₂ xᵀQx + aᵀx + b`, with `Q` symmetric...
# (https://jump.dev/MathOptInterface.jl/stable/reference/standard_form/#MathOptInterface.ScalarQuadraticFunction)
# ... whereas cuOpt expects a QP objective of the form `¹xᵀQx + aᵀx + b`,
# where `Q` need not be symmetric
# --> we need to scale diagonal coeffs. by ¹/₂ to match cuOpt convention
if i == j
v /= 2
Copy link
Collaborator

@rg20 rg20 Jan 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure if this is accurate.

cuOpt expects users to provide the true objective function.

For example: if you are minimizing x1^2 + x2^2, the matrix should be [1 0; 0 1], cuOpt internally minimizes for (1/2) [x1 x2] [2 0; 0 2] [x1; x2] in this case.

Copy link
Member

@odow odow Jan 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cuOpt expects users to provide the true objective function.

https://www.youtube.com/watch?v=M31xoZGyj9w&t=2698s

Copy link
Member

@odow odow Jan 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happens to the off-diagonal terms?

Copy link
Collaborator

@rg20 rg20 Jan 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

cuOpt is not necessarily expecting the matrix to be symmetric.

If the objective is xT Q x + cT x, we internally symmetrize and solve for (1/2) xT (Q + QT) x + cT x

So if the objective is x1^2 + x2^2 + 2 x1 x2,

Q can be: [1 2; 0 1], or [1 0; 2 1] or [1 1; 1; 1]

Copy link
Author

@mtanneau mtanneau Jan 12, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(including the links for completeness)

I believe the /2 factor for diagonal terms is required for correctness, but agreed it's not super clear why at first.

  1. The current MOI wrapper uses a JuMP-level data structure to store the optimization model as the user builds it. That internal representation abides by MOI conventions.
  2. when optimize! gets called, the problem data is flushed from the JuMP-level data store into a cuOpt-level data structure --> this is why this PR mostly modifies the copy_to function.
  3. Hence, when we access the quadratic objective in copy_to, we get as input an MOI.ScalarQuadraticFunction, which abides by the MOI convention described in the docs above.

Here is a small example to illustrate this

using JuMP

model = Model()
@variable(model, x)
@variable(model, y)
@objective(model, Min, x*x + 2*x*y + 3*y*y)
objective_function(model)  # x² + 2 x*y + 3 y²

# Now, access the MOI representation
F = MOI.get(model, MOI.ObjectiveFunctionType())
f = MOI.get(model, MOI.ObjectiveFunction{F}())
f.quadratic terms

the last line outputs

3-element Vector{MathOptInterface.ScalarQuadraticTerm{Float64}}:
 MathOptInterface.ScalarQuadraticTerm{Float64}(2.0, MOI.VariableIndex(1), MOI.VariableIndex(1))
 MathOptInterface.ScalarQuadraticTerm{Float64}(2.0, MOI.VariableIndex(1), MOI.VariableIndex(2))
 MathOptInterface.ScalarQuadraticTerm{Float64}(6.0, MOI.VariableIndex(2), MOI.VariableIndex(2))

--> you can see that the diagonal coefficients were multiplied by 2. This was done by MOI in accordance with its convention.

Copy link
Author

@mtanneau mtanneau Jan 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@rg20 I can add the following:

  • add link to the MOI docs on quadratic function in that comment, to provide additional context
  • add a couple of unit tests to validate the cuOpt solution and objective value when solving QP from MOI

Copy link
Member

@odow odow Jan 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My conclusion with this all is that the only valid way is to test, test, test, test, and test. There are too many subtleties to try and logically reason about the transformations, and every time I do, I end up making a mistake.

There are tests in MOI for various cases, but these are precisely the ones that you're skipping because of the upstream bug... 😢 ("test_objective_qp_ObjectiveFunction_zero_ofdiag" and "test_objective_qp_ObjectiveFunction_edge_cases").

Copy link
Collaborator

@rg20 rg20 Jan 13, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mtanneau thanks for explaining the subtleties in the MOI wrapper.

Copy link
Author

@mtanneau mtanneau Jan 13, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added more unit tests that specifically trigger the QP objective.
Notes:

  • the upstream issue regarding QP with no constraints might get fixed in cuopt v26.02 (tracking Use augmented system when there are no constraints NVIDIA/cuopt#765), which would allow to run more MOI-level tests.
  • if the team is OK with adding JuMP as a test dependency, I'm happy to add similar tests where the QP is built from JuMP, not from MOI.

Copy link
Member

@mlubin mlubin Jan 14, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

JuMP shouldn't be a test dependency; MOI-style tests should be sufficient

end

# We are building a COO of Qᵀ --> swap i and j
push!(Qtrows, j)
push!(Qtcols, i)
push!(Qtvals, v)
end
# CSC of Qᵀ is CSR of Q
Qt = SparseArrays.sparse(Qtrows, Qtcols, Qtvals, numcol, numcol)
qobj_matrix_values = Qt.nzval
qobj_row_offsets = Qt.colptr
qobj_col_indices = Qt.rowval
# Revert to 0-based indexing
qobj_row_offsets .-= Int32(1)
qobj_col_indices .-= Int32(1)

return objective_offset,
objective_coefficients_linear,
qobj_row_offsets,
qobj_col_indices,
qobj_matrix_values
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
Expand Down Expand Up @@ -970,27 +1056,63 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
has_integrality = true
end

objective_sense, objective_offset, objective_coefficients =
_get_objective_data(dest, src, mapping, numcol)

ref_problem = Ref{cuOptOptimizationProblem}()
ret = cuOptCreateRangedProblem(
numrow,
numcol,
objective_sense,
objective_offset,
objective_coefficients,
constraint_matrix_row_offsets,
constraint_matrix_column_indices,
constraint_matrix_coefficients,
rowlower,
rowupper,
collower,
colupper,
var_type,
ref_problem,
)
_check_ret(ret, "cuOptCreateRangedProblem")
objective_sense,
objective_offset,
objective_coefficients,
qobj_row_offsets,
qobj_col_indices,
qobj_matrix_values = _get_objective_data(dest, src, mapping, numcol)

# Is this a QP or an LP?
has_quadratic_objective = length(qobj_matrix_values) > 0
if has_quadratic_objective && has_integrality
error(
"cuOpt does not support models with quadratic objectives _and_ integer variables",
)
end

if has_quadratic_objective
ref_problem = Ref{cuOptOptimizationProblem}()
ret = cuOptCreateQuadraticRangedProblem(
numrow,
numcol,
objective_sense,
objective_offset,
objective_coefficients,
qobj_row_offsets,
qobj_col_indices,
qobj_matrix_values,
constraint_matrix_row_offsets,
constraint_matrix_column_indices,
constraint_matrix_coefficients,
rowlower,
rowupper,
collower,
colupper,
ref_problem,
)
_check_ret(ret, "cuOptCreateQuadraticRangedProblem")
else
ref_problem = Ref{cuOptOptimizationProblem}()
ret = cuOptCreateRangedProblem(
numrow,
numcol,
objective_sense,
objective_offset,
objective_coefficients,
constraint_matrix_row_offsets,
constraint_matrix_column_indices,
constraint_matrix_coefficients,
rowlower,
rowupper,
collower,
colupper,
var_type,
ref_problem,
)
_check_ret(ret, "cuOptCreateRangedProblem")
end

dest.cuopt_problem = ref_problem[]

ref_settings = Ref{cuOptSolverSettings}()
Expand Down
2 changes: 1 addition & 1 deletion src/cuOpt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ function __init__()
error("Failed to get cuOpt library version (status code: $status)")
end
version = VersionNumber(major[], minor[], patch[])
min, max = v"25.08", v"25.13"
min, max = v"25.12", v"25.13"
if !(min <= version < max)
error(
"Incompatible cuOpt library version. Got $version, but supported versions are [$min, $max)",
Expand Down
132 changes: 132 additions & 0 deletions test/MOI_wrapper.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,10 @@ function test_runtests_cache_optimizer()
"test_constraint_ZeroOne_bounds_3",
# Upstream bug: https://github.com/NVIDIA/cuopt/issues/112
"test_solve_TerminationStatus_DUAL_INFEASIBLE",
# Upstream bug: https://github.com/NVIDIA/cuopt/issues/759
# (cuOpt crashes when given a QP with no linear constraints)
"test_objective_qp_ObjectiveFunction_zero_ofdiag",
"test_objective_qp_ObjectiveFunction_edge_cases",
],
)
return
Expand All @@ -78,6 +82,134 @@ function test_air05()
return
end

function test_qp_objective()
model = MOI.instantiate(cuOpt.Optimizer; with_cache_type = Float64)
MOI.set(model, MOI.Silent(), true)
x, _ = MOI.add_constrained_variable(
model,
(MOI.GreaterThan(0.0), MOI.LessThan(1.0)),
)
y, _ = MOI.add_constrained_variable(
model,
(MOI.GreaterThan(0.0), MOI.LessThan(1.0)),
)

# x + y == 1
MOI.add_constraint(
model,
MOI.ScalarAffineFunction(
[MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(1.0, y)],
0.0,
),
MOI.EqualTo(1.0),
)

F = MOI.ScalarQuadraticFunction{Float64}
MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

# 1. Homogeneous QP with diagonal objective
# Min ¹/₂ * (x² + y²) s.t. x+y == 1
fobj = MOI.ScalarQuadraticFunction(
[
MOI.ScalarQuadraticTerm(1.0, x, x),
MOI.ScalarQuadraticTerm(1.0, y, y),
],
MOI.ScalarAffineTerm{Float64}[],
0.0,
)
MOI.set(model, MOI.ObjectiveFunction{F}(), fobj)
MOI.optimize!(model)

@test isapprox(
MOI.get(model, MOI.VariablePrimal(), x),
0.5;
atol = 1e-4,
rtol = 1e-4,
)
@test isapprox(
MOI.get(model, MOI.VariablePrimal(), y),
0.5;
atol = 1e-4,
rtol = 1e-4,
)
# ¹/₂ * (2 * 0.5² + 2*0.5²) == 0.5
@test isapprox(
MOI.get(model, MOI.ObjectiveValue()),
0.25;
atol = 1e-4,
rtol = 1e-4,
)

# Change diagonal coefficients
# Min ¹/₂ * (2x² + 2y²) s.t. x+y == 1
fobj = MOI.ScalarQuadraticFunction(
[
MOI.ScalarQuadraticTerm(2.0, x, x),
MOI.ScalarQuadraticTerm(2.0, y, y),
],
MOI.ScalarAffineTerm{Float64}[],
0.0,
)
MOI.set(model, MOI.ObjectiveFunction{F}(), fobj)
MOI.optimize!(model)
# Same solution, but different objective value
# ¹/₂ * (2 * 0.5² + 2*0.5²) == 0.5
@test isapprox(
MOI.get(model, MOI.ObjectiveValue()),
0.5;
atol = 1e-4,
rtol = 1e-4,
)
@test isapprox(
MOI.get(model, MOI.VariablePrimal(), x),
0.5;
atol = 1e-4,
rtol = 1e-4,
)
@test isapprox(
MOI.get(model, MOI.VariablePrimal(), y),
0.5;
atol = 1e-4,
rtol = 1e-4,
)

# QP with off-diagonal term
# (x - y - 2)² = x² - 2xy + y² - 4x + 4y + 4
# = ¹/₂ (2x² - xy - yx + 2y²) + (-4x + 4y) + 4
# Solution is (x, y) = (1, 0)
fobj = MOI.ScalarQuadraticFunction(
[
MOI.ScalarQuadraticTerm(2.0, x, x),
MOI.ScalarQuadraticTerm(2.0, y, y),
MOI.ScalarQuadraticTerm(-1.0, x, y),
],
[MOI.ScalarAffineTerm(-4.0, x), MOI.ScalarAffineTerm(+4.0, y)],
4.0,
)
MOI.set(model, MOI.ObjectiveFunction{F}(), fobj)
MOI.optimize!(model)
@test isapprox(
MOI.get(model, MOI.ObjectiveValue()),
1;
atol = 1e-4,
rtol = 1e-4,
)
@test isapprox(
MOI.get(model, MOI.VariablePrimal(), x),
1.0;
atol = 1e-4,
rtol = 1e-4,
)
@test isapprox(
MOI.get(model, MOI.VariablePrimal(), y),
0.0;
atol = 1e-4,
rtol = 1e-4,
)

return nothing
end

end # TestMOIWrapper

TestMOIWrapper.runtests()