Note
JAX is not on the default solve path. It remains an optional dependency of the
differentiable-solve and learned-relaxation subsystems, but a plain solve()
imports no jax module: model derivatives come from the POUNCE AD tape in the Rust
core.
Solver Backend Selection Guide#
discopt ships with two NLP solver backends:
POUNCE (
nlp_solver="pounce", the default) – a pure-Rust port of the Ipopt interior-point method [Mehrotra, 1992, Wächter and Biegler, 2006] exposed via PyO3. A core dependency, so nothing extra is needed.cyipopt / Ipopt (
nlp_solver="ipopt") – the reference C++ Ipopt [Wächter and Biegler, 2006] accessed through cyipopt. Sparse linear algebra (MA27/MA57/MUMPS) for large problems; needs libipopt and BLAS/LAPACK.
nlp_solver="ipm" and nlp_solver="sparse_ipm" are back-compat aliases that
resolve to POUNCE, not a third engine: the pure-JAX Mehrotra IPM they once named
is retired and removed.
For LP and QP sub-classes the solver auto-dispatches to specialized routines that bypass the general NLP machinery entirely – a pure LP or MILP goes to HiGHS, and a QP to POUNCE’s QP path.
This tutorial walks through the trade-offs with head-to-head benchmarks so you can make an informed choice. The theoretical foundations of interior-point methods are covered in [Nocedal and Wright, 2006] and [Wright, 1997].
import time
import discopt.modeling as dm
import numpy as np
Backend Overview#
Feature |
POUNCE ( |
cyipopt / Ipopt ( |
|---|---|---|
Language |
Rust (PyO3) |
C++ / Fortran |
Method |
Filter interior-point (an Ipopt port) |
Filter interior-point (the reference) |
Sparse linear algebra |
No |
Yes (MA27 / MA57 / MUMPS) |
External dependencies |
None (a core dependency) |
libipopt, BLAS/LAPACK |
Best for |
The default everywhere |
Large, sparse NLP; cross-checking |
"ipm" and "sparse_ipm" are not rows in this table: they are aliases for the
POUNCE column. A cell below asserts that, rather than leaving it as prose.
Automatic Problem Classification#
Before dispatching to any NLP backend, discopt classifies the problem using Rust structure-detection on the expression DAG. LP and QP problems are routed to specialized solvers automatically, regardless of the nlp_solver setting.
from discopt._relax.problem_classifier import classify_problem
# --- LP ---
m_lp = dm.Model("lp_demo")
x_lp = m_lp.continuous("x", shape=(3,), lb=0)
m_lp.minimize(2 * x_lp[0] + 3 * x_lp[1] + x_lp[2])
m_lp.subject_to(x_lp[0] + x_lp[1] + x_lp[2] >= 1)
print(f"LP : {classify_problem(m_lp)}")
# --- QP ---
m_qp = dm.Model("qp_demo")
x_qp = m_qp.continuous("x", shape=(3,), lb=0)
m_qp.minimize(x_qp[0] ** 2 + x_qp[1] ** 2 + x_qp[0] * x_qp[1])
m_qp.subject_to(x_qp[0] + x_qp[1] + x_qp[2] >= 1)
print(f"QP : {classify_problem(m_qp)}")
# --- NLP ---
m_nlp = dm.Model("nlp_demo")
x_nlp = m_nlp.continuous("x", shape=(3,), lb=-5, ub=5)
m_nlp.minimize(dm.sin(x_nlp[0]) + x_nlp[1] ** 2)
m_nlp.subject_to(x_nlp[0] + x_nlp[1] >= 1)
print(f"NLP : {classify_problem(m_nlp)}")
# --- MILP ---
m_milp = dm.Model("milp_demo")
y = m_milp.binary("y", shape=(2,))
x_milp = m_milp.continuous("x", shape=(2,), lb=0, ub=10)
m_milp.minimize(3 * y[0] + 2 * y[1] + x_milp[0])
m_milp.subject_to(x_milp[0] + x_milp[1] + y[0] >= 2)
print(f"MILP: {classify_problem(m_milp)}")
# --- MINLP ---
m_minlp = dm.Model("minlp_demo")
y2 = m_minlp.binary("y", shape=(2,))
x_minlp = m_minlp.continuous("x", shape=(2,), lb=0, ub=10)
m_minlp.minimize(dm.exp(x_minlp[0]) + y2[0] * x_minlp[1])
m_minlp.subject_to(x_minlp[0] + x_minlp[1] >= 1)
print(f"MINLP: {classify_problem(m_minlp)}")
LP : ProblemClass.LP
QP : ProblemClass.QP
NLP : ProblemClass.NLP
MILP: ProblemClass.MILP
MINLP: ProblemClass.MINLP
Auto-dispatch rules:
Classification |
Solver route |
|---|---|
LP |
HiGHS, with a discopt-verified certificate (#1229); opt out with |
QP |
POUNCE’s specialized QP path (augmented KKT) |
NLP |
The chosen |
MILP |
HiGHS (same route as LP) |
MIQP / MINLP |
Branch-and-Bound, or the MIP–NLP family when the model certifies convex at the root; the in-house Rust simplex solves each node’s LP relaxation |
Benchmark 1: Small NLP (n = 5)#
A small quadratic-plus-bilinear objective with a linear constraint. All three backends should return similar solutions quickly.
def build_small_nlp():
m = dm.Model("small_nlp")
x = m.continuous("x", shape=(5,), lb=-5.0, ub=5.0)
obj_expr = sum((x[i] - i) ** 2 for i in range(5)) + x[0] * x[1]
m.minimize(obj_expr)
m.subject_to(x[0] + x[1] >= 1.0)
return m
# Two engines, not three: "ipm" is an alias for "pounce". It is listed here so
# the alias is *demonstrated* -- the cell asserts below that the two spellings
# return the identical objective -- not so it can be read as a third backend.
backends = ["pounce", "ipm", "ipopt"]
# One warm-up solve per backend, discarded. Without it the first row absorbs
# every one-time import and setup cost and reads as a 30x difference that is not
# a property of the solver (the committed output of this cell showed exactly
# that: "ipm" 0.3374 s against "pounce" 0.0103 s, for the same engine).
for _b in backends:
build_small_nlp().solve(nlp_solver=_b)
objectives = {}
print(f"{'Backend':<12} {'Objective':>12} {'Wall time (s)':>14} {'Status'}")
print("-" * 52)
for backend in backends:
model = build_small_nlp()
t0 = time.perf_counter()
result = model.solve(nlp_solver=backend)
elapsed = time.perf_counter() - t0
print(f"{backend:<12} {result.objective:12.6f} {elapsed:14.4f} {result.status}")
objectives[backend] = result.objective
assert objectives["ipm"] == objectives["pounce"], (
'nlp_solver="ipm" is an alias for POUNCE, so the two must be bit-identical: '
f"{objectives['ipm']} vs {objectives['pounce']}"
)
Backend Objective Wall time (s) Status
----------------------------------------------------
pounce -0.250000 0.0041 optimal
ipm -0.250000 0.0042 optimal
ipopt -0.250000 0.0042 optimal
At this scale both engines converge in well under a second, on the identical
objective. "ipm" and "pounce" agree to the bit because they are the same
engine — the assertion in the cell is what keeps that from drifting back into a
three-way comparison. The timings are single un-replicated solves after a warm-up
and should not be read as a benchmark; see
Benchmarks by Problem Class for measured comparisons.
Benchmark 2: Medium NLP (n = 20)#
A quadratic objective with nonlinear (sine) coupling constraints.
def build_medium_nlp(n=20):
m = dm.Model("medium_nlp")
x = m.continuous("x", shape=(n,), lb=-5.0, ub=5.0)
obj_expr = sum((x[i] - 0.5 * i) ** 2 for i in range(n))
m.minimize(obj_expr)
# Coupling constraints
for i in range(n - 1):
m.subject_to(x[i] + x[i + 1] >= 0.5)
# One nonlinear constraint to keep it NLP
m.subject_to(dm.sin(x[0]) + x[1] >= -1.0)
return m
print(f"{'Backend':<12} {'Objective':>12} {'Wall time (s)':>14} {'Status'}")
print("-" * 52)
for backend in backends:
model = build_medium_nlp()
t0 = time.perf_counter()
result = model.solve(nlp_solver=backend)
elapsed = time.perf_counter() - t0
print(f"{backend:<12} {result.objective:12.6f} {elapsed:14.4f} {result.status}")
Backend Objective Wall time (s) Status
----------------------------------------------------
pounce 71.250000 0.2297 optimal
ipm 71.250000 0.2091 optimal
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
ipopt 71.250000 0.2213 optimal
At medium scale the backends remain competitive. Ipopt’s sparse linear algebra has not yet overtaken JAX’s dense JIT-compiled kernels.
Benchmark 3: Larger NLP (n = 50)#
A random quadratic objective with coupling constraints. As problem size grows, Ipopt’s sparse direct solvers (MA27/MA57) begin to show an advantage over dense JIT kernels.
def build_larger_nlp(n=50, seed=42):
rng = np.random.default_rng(seed)
m = dm.Model("larger_nlp")
x = m.continuous("x", shape=(n,), lb=-5.0, ub=5.0)
# Random positive-definite quadratic
targets = rng.standard_normal(n)
obj_expr = sum((x[i] - float(targets[i])) ** 2 for i in range(n))
m.minimize(obj_expr)
# Chain constraints
for i in range(n - 1):
m.subject_to(x[i] + x[i + 1] >= -1.0)
# Nonlinear constraint
m.subject_to(dm.sin(x[0]) + dm.cos(x[1]) >= -1.5)
return m
print(f"{'Backend':<12} {'Objective':>12} {'Wall time (s)':>14} {'Status'}")
print("-" * 52)
for backend in backends:
model = build_larger_nlp()
t0 = time.perf_counter()
result = model.solve(nlp_solver=backend)
elapsed = time.perf_counter() - t0
print(f"{backend:<12} {result.objective:12.6f} {elapsed:14.4f} {result.status}")
Backend Objective Wall time (s) Status
----------------------------------------------------
pounce 2.812503 0.6561 optimal
ipm 2.812503 0.6230 optimal
ipopt 2.812503 0.6220 optimal
For problems beyond roughly 200 variables, cyipopt with sparse linear algebra is
typically the fastest single-instance solver; the crossover depends on sparsity
structure. On this instance the three rows are within 2 % of each other, which is
what a single un-replicated solve can tell you and no more.
Benchmark 4: Batch of 16 Small NLPs#
Solving many instances of the same structure with different data is a serial loop
over solve() on every current backend. There is no vectorized batch entry point:
the jax.vmap-batched LP/QP IPM this section used to advertise
(lp_ipm_solve_batch / qp_ipm_solve_batch) went with the JAX interior-point
stack and does not exist — importing it raises ModuleNotFoundError. What
batching discopt does have is internal: the spatial B&B evaluates a batch of
nodes per iteration (batch_size, default 16).
batch_size = 16
rng = np.random.default_rng(123)
# Build 16 variants of the small NLP with different targets
models = []
for b in range(batch_size):
m = dm.Model(f"batch_{b}")
x = m.continuous("x", shape=(5,), lb=-5.0, ub=5.0)
targets = rng.standard_normal(5)
obj_expr = sum((x[i] - float(targets[i])) ** 2 for i in range(5))
m.minimize(obj_expr)
m.subject_to(x[0] + x[1] >= 0.5)
models.append(m)
# Serial solve with each backend
for backend in backends:
t0 = time.perf_counter()
results = [m.solve(nlp_solver=backend) for m in models]
elapsed = time.perf_counter() - t0
# `status` is the lowercase string "optimal". Comparing against "OPTIMAL"
# made this counter zero in every row of the committed output -- 0/16 three
# times, printed as though it were a result.
n_opt = sum(1 for r in results if r.status == "optimal")
assert n_opt == batch_size, f"{backend}: only {n_opt}/{batch_size} solved"
print(f"{backend:<12} serial {batch_size}x : {elapsed:.4f}s ({n_opt}/{batch_size} optimal)")
pounce serial 16x : 0.0535s (16/16 optimal)
ipm serial 16x : 0.0497s (16/16 optimal)
ipopt serial 16x : 0.0458s (16/16 optimal)
# The same batch again after a warm-up pass, so the numbers above are not
# dominated by one-time setup. This is a serial loop either way: there is no
# vectorized batch NLP entry point (see the note above).
_ = [m.solve(nlp_solver="pounce") for m in models]
t0 = time.perf_counter()
results_warm = [m.solve(nlp_solver="pounce") for m in models]
elapsed_warm = time.perf_counter() - t0
n_opt = sum(1 for r in results_warm if r.status == "optimal")
print(f"POUNCE (warm) serial {batch_size}x : {elapsed_warm:.4f}s ({n_opt}/{batch_size} optimal)")
assert n_opt == batch_size, f"only {n_opt}/{batch_size} solved"
POUNCE (warm) serial 16x : 0.0705s (16/16 optimal)
Warm and cold differ little here: POUNCE is compiled Rust with no per-shape compilation step to amortize, so repeating the batch buys nothing beyond process warm-up. If you need throughput over many independent instances, parallelize at the process level — there is no in-process vectorized batch solve.
Benchmark 5: Small MINLP#
For mixed-integer problems, the NLP backend is called at every node of the Branch-and-Bound tree [Belotti et al., 2013]. A faster NLP backend directly reduces total wall time.
def build_small_minlp():
m = dm.Model("small_minlp")
y = m.binary("y", shape=(2,))
x = m.continuous("x", shape=(3,), lb=0.0, ub=5.0)
obj = (x[0] - 1) ** 2 + (x[1] - 2) ** 2 + x[2] + 3 * y[0] + 2 * y[1]
m.minimize(obj)
m.subject_to(x[0] + x[1] + y[0] >= 2.0)
m.subject_to(x[1] + x[2] + y[1] >= 1.5)
m.subject_to(x[0] + y[0] + y[1] <= 4.0)
return m
print(f"{'Backend':<12} {'Objective':>12} {'Nodes':>8} {'Wall time (s)':>14} {'Status'}")
print("-" * 60)
for backend in backends:
model = build_small_minlp()
t0 = time.perf_counter()
result = model.solve(nlp_solver=backend)
elapsed = time.perf_counter() - t0
nodes = getattr(result, "node_count", "N/A")
print(f"{backend:<12} {result.objective:12.6f} {nodes:>8} {elapsed:14.4f} {result.status}")
OA: generating OA cuts for 0 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)
OA: generating OA cuts for 0 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)
OA: generating OA cuts for 0 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)
Backend Objective Nodes Wall time (s) Status
------------------------------------------------------------
pounce 0.000000 0 0.0415 optimal
ipm 0.000000 0 0.0436 optimal
ipopt 0.000000 0 0.0146 optimal
Node counts are identical across backends, as they must be — the tree is
deterministic and the backend only changes how each node’s NLP is solved. Here
every row reads 0 nodes: this MINLP certifies convex at the root, so it is routed
to the MIP–NLP family (outer approximation) and never opens a B&B node. The
comparison is then of the route’s NLP subsolves, not of per-node solves.
Decision Flowchart#
Is it LP or MILP?
|-- Yes --> Auto-dispatched to HiGHS, certificate verified by discopt (#1229)
|-- No --> Is it a QP?
|-- Yes --> Auto-dispatched to POUNCE's QP path
|-- No --> Large and very sparse (>200 vars)?
|-- Yes --> "ipopt" (sparse MA27 / MA57 / MUMPS)
|-- No --> "pounce" (the default; pure Rust, no
external toolchain)
There is no batch branch: nlp_solver chooses between two engines, and the
classification above is made for you regardless of which one you name.
Sparse IPM (retired alias)#
nlp_solver="sparse_ipm" was a JAX-IPM variant that paired JAX autodiff
with a SciPy sparse KKT factorization. The JAX interior-point stack has
since been removed; "sparse_ipm" is now a deprecated alias that resolves
to POUNCE. For medium-to-large sparse NLPs use nlp_solver="pounce" (the
default, pure-Rust) or nlp_solver="ipopt" (sparse MA27/MA57) instead.
Summary: Backend Recommendations#
Scenario |
Recommended backend |
Why |
|---|---|---|
LP / QP (any size) |
auto-dispatch |
Specialized IPM (POUNCE) / HiGHS oracle |
NLP (any size) |
|
Pure-Rust IPM, no external toolchain |
Large, very sparse NLP |
|
Sparse MA27/MA57 factorization |
MILP / MIQP |
auto-dispatch |
Warm-started Rust simplex / IPM at each B&B node |
MINLP |
|
Robust IPM node solves + spatial B&B |
Dependency-free deployment |
|
Pure Rust, no C/Fortran |
"ipm" / "sparse_ipm" remain accepted but alias to "pounce".
Exercise: Predict the Best Backend#
For each scenario below, predict which backend you would choose, then verify by timing.
(a) You need to solve 100 small LPs (5 variables each) as fast as possible.
(b) You have a single NLP with 100 continuous variables and sparse constraints.
(c) You are deploying a small MINLP solver to an embedded system with no C toolchain.
# Exercise skeleton -- fill in your predictions as comments, then run to verify.
# (a) 100 small LPs -- Prediction: _________
lp_models = []
rng_ex = np.random.default_rng(99)
for i in range(100):
m = dm.Model(f"lp_{i}")
x = m.continuous("x", shape=(5,), lb=0.0)
c = rng_ex.standard_normal(5)
m.minimize(sum(float(c[j]) * x[j] for j in range(5)))
m.subject_to(sum(x[j] for j in range(5)) >= 1.0)
lp_models.append(m)
t0 = time.perf_counter()
lp_results = [m.solve() for m in lp_models] # auto-dispatches to HiGHS (#1229)
print(f"(a) 100 LPs solved in {time.perf_counter() - t0:.4f}s (auto-dispatch to HiGHS)")
assert all(r.status == "optimal" for r in lp_results), "not every LP solved"
/Users/jkitchin/projects/discopt/python/discopt/debug/__init__.py:98: UserWarning: Variables with very large or infinite declared bounds: x[0] (lb=0, ub=9.999e+19 [default]), x[1] (lb=0, ub=9.999e+19 [default]), x[2] (lb=0, ub=9.999e+19 [default]), x[3] (lb=0, ub=9.999e+19 [default]), x[4] (lb=0, ub=9.999e+19 [default]). A bound marked [default] is the box applied to a column you declared with no bounds: 9.999e+19, which is FINITE (it sits just below the 1e+20 infinity sentinel), so the solve can return a certified 'optimal' sitting on that corner rather than 'unbounded'. NLP solvers may fail (NaN, iteration_limit) when bounds exceed ~1e15. Add tighter explicit bounds, e.g. m.continuous('x', lb=0, ub=1000).
return fn(*args, **kwargs)
(a) 100 LPs solved in 0.3693s (auto-dispatch to HiGHS)
# (b) Single NLP, n=100 -- Prediction: _________
def build_nlp_100():
m = dm.Model("nlp_100")
x = m.continuous("x", shape=(100,), lb=-5.0, ub=5.0)
rng_b = np.random.default_rng(7)
targets = rng_b.standard_normal(100)
obj_expr = sum((x[i] - float(targets[i])) ** 2 for i in range(100))
m.minimize(obj_expr)
for i in range(99):
m.subject_to(x[i] + x[i + 1] >= -2.0)
m.subject_to(dm.sin(x[0]) + x[1] >= -1.0)
return m
for backend in backends:
model = build_nlp_100()
t0 = time.perf_counter()
result = model.solve(nlp_solver=backend)
elapsed = time.perf_counter() - t0
print(f"(b) n=100 with {backend:<8}: {elapsed:.4f}s obj={result.objective:.4f}")
(b) n=100 with pounce : 1.8682s obj=3.0741
(b) n=100 with ipm : 1.8010s obj=3.0741
(b) n=100 with ipopt : 1.7825s obj=3.0741
# (c) Small MINLP for deployment -- Prediction: _________
# Answer: pounce -- pure Rust, no C dependencies, and the default. There is no
# fallback to name: it is one of only two engines, and the other one (ipopt) is
# the C/Fortran dependency this scenario rules out.
model = build_small_minlp()
t0 = time.perf_counter()
result = model.solve(nlp_solver="ipm") # change to "pounce" if available
elapsed = time.perf_counter() - t0
nodes = getattr(result, "node_count", "N/A")
print(f"(c) Small MINLP: {elapsed:.4f}s obj={result.objective:.4f} nodes={nodes}")
OA: generating OA cuts for 0 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)
(c) Small MINLP: 0.0149s obj=0.0000 nodes=0
Exercise answers:
(a) Auto-dispatch to the specialized LP IPM. No
nlp_solverargument needed – discopt detects the LP structure and routes accordingly.(b) For 100 variables,
"pounce"(the default) is the recommended choice;"ipopt"can be faster on very sparse problems via MA27.(c)
"pounce"– pure Rust with zero external C/Fortran dependencies. Ideal for containerized or embedded deployment.
Summary#
Choosing the right solver backend in discopt is straightforward:
Let auto-dispatch handle LP and QP problems – the specialized IPM solvers are fastest for these classes.
Use POUNCE (the default) for NLPs, MINLPs, and batch workloads – pure-Rust interior point with no external toolchain; batch node solves in Rust replace the old JAX
vmappath.Switch to cyipopt/Ipopt for large, very sparse NLPs – the MA27/MA57 sparse direct solvers scale well there.
"ipm"/"sparse_ipm"still work – they are back-compat aliases that route to POUNCE (the JAX IPM was retired).