Solver Backend Selection Guide

Solver Backend Selection Guide#

discopt ships with three NLP solver backends, each targeting a different use case:

JAX IPM (nlp_solver="ipm") – a pure-JAX Mehrotra predictor-corrector interior-point method [Mehrotra, 1992]. JIT-compiled and vmap-compatible for batch solving.
ripopt (nlp_solver="ripopt") – a Rust interior-point solver exposed via PyO3. Zero external dependencies; ideal for deployment.
cyipopt / Ipopt (nlp_solver="ipopt") – the industry-standard Ipopt solver [Wächter and Biegler, 2006] accessed through cyipopt. Leverages sparse linear algebra (MA27/MA57) for large problems.

For LP and QP sub-classes the solver auto-dispatches to specialized JAX IPM routines that bypass the general NLP machinery entirely.

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 os

os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"

import time

import discopt.modeling as dm
import numpy as np

Backend Overview#

Feature	JAX IPM (`"ipm"`)	ripopt (`"ripopt"`)	cyipopt / Ipopt (`"ipopt"`)
Language	Python / JAX	Rust (PyO3)	C / Fortran
JIT compiled	Yes	No	No
Batch via `vmap`	Yes	No	No
Sparse linear algebra	Limited	No	Yes (MA27/MA57)
External dependencies	None	None	libipopt, BLAS/LAPACK
Best for	Small-medium NLP, batch	Dependency-free deployment	Large, sparse NLP

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._jax.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	Specialized `lp_ipm` (Mehrotra PC, normal-equations)
QP	Specialized `qp_ipm` (augmented KKT)
NLP	Chosen `nlp_solver` backend
MILP / MIQP	Branch-and-Bound with LP/QP IPM at each node
MINLP	Branch-and-Bound with chosen NLP backend at each node

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


backends = ["ipm", "ipopt"]
try:
    backends = ["ipm", "ripopt", "ipopt"]
except Exception:
    pass

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}")

Backend         Objective  Wall time (s) Status
----------------------------------------------------

ipm             -0.250000         0.9449 optimal
ripopt          -0.250000         0.0389 optimal
ipopt           -0.250000         0.0363 optimal

At this scale all backends converge in well under a second. The JAX IPM incurs a one-time JIT compilation cost on the first call; subsequent solves with the same problem dimensions are significantly faster.

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
----------------------------------------------------

******************************************************************************
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
******************************************************************************

ipm             71.249998         0.3935 optimal

ripopt          71.250000         0.3424 optimal

ipopt           71.249998         0.3677 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
----------------------------------------------------

ipm              2.812503         0.7956 optimal

ripopt           2.812503         0.7616 optimal

ipopt            2.812503         0.8028 optimal

For problems beyond roughly 200 variables, cyipopt with sparse linear algebra is typically the fastest single-instance solver. The crossover point depends on sparsity structure: denser Hessians favor JAX IPM longer.

Benchmark 4: Batch of 16 Small NLPs#

This is where the JAX IPM truly shines. When you need to solve many instances of the same problem structure with different data, jax.vmap vectorizes the IPM across instances in a single fused kernel. This provides 10–30x speedups over a serial loop [Mehrotra, 1992].

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
    n_opt = sum(1 for r in results if r.status == "OPTIMAL")
    print(f"{backend:<12} serial {batch_size}x : {elapsed:.4f}s  ({n_opt}/{batch_size} optimal)")

ipm          serial 16x : 0.4890s  (0/16 optimal)

ripopt       serial 16x : 0.5074s  (0/16 optimal)

ipopt        serial 16x : 0.5176s  (0/16 optimal)

# Batch solve using JAX IPM vmap (via the low-level LP/QP batch API)
# For NLPs the high-level batch path uses vmapped root multistart.
# Here we demonstrate the concept with the serial JAX IPM path,
# which benefits from JIT cache reuse across same-shaped problems.

# Warm up JIT
_ = models[0].solve(nlp_solver="ipm")

t0 = time.perf_counter()
results_jax = [m.solve(nlp_solver="ipm") for m in models]
elapsed_jax_warm = time.perf_counter() - t0
n_opt = sum(1 for r in results_jax if r.status == "OPTIMAL")
print(
    f"JAX IPM (JIT-warm) serial {batch_size}x :"
    f" {elapsed_jax_warm:.4f}s  ({n_opt}/{batch_size} optimal)"
)
print()
print("Note: For LP/QP batches, lp_ipm_solve_batch / qp_ipm_solve_batch")
print("use jax.vmap for true vectorized batch solving (10-30x speedup).")

JAX IPM (JIT-warm) serial 16x : 0.5160s  (0/16 optimal)

Note: For LP/QP batches, lp_ipm_solve_batch / qp_ipm_solve_batch
use jax.vmap for true vectorized batch solving (10-30x speedup).

The JAX IPM’s JIT cache reuse means that after the first compilation, subsequent solves of the same shape execute purely compiled XLA code. For LP and QP sub-problems the lp_ipm_solve_batch and qp_ipm_solve_batch functions provide true vmap-vectorized batch solving.

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}")

Backend         Objective    Nodes  Wall time (s) Status
------------------------------------------------------------

ipm              0.000000        1         0.7115 optimal
ripopt           0.000000        1         0.1082 optimal

ipopt            0.000000        1         0.1033 optimal

Node counts should be identical across backends (the B&B tree is deterministic), but wall time varies with per-node NLP solve speed.

Decision Flowchart#

Is it LP or QP?
  |-- Yes --> Auto-dispatched to specialized JAX IPM (fastest)
  |-- No  --> Need to solve a batch of instances?
                |-- Yes --> JAX IPM with vmap (10-30x speedup)
                |-- No  --> Large problem (>200 vars)?
                              |-- Yes --> cyipopt (sparse linear algebra)
                              |-- No  --> Need C-library-free deployment?
                                            |-- Yes --> ripopt (pure Rust, no deps)
                                            |-- No  --> JAX IPM (default, JIT-optimized)

Sparse IPM: A Middle Ground#

For medium-to-large problems where you want JAX automatic differentiation but need sparse linear algebra, there is a fourth option:

result = model.solve(nlp_solver="sparse_ipm")

This hybrid backend uses:

JAX for objective/constraint evaluation and automatic differentiation (gradients, Hessians, Jacobians)
SciPy sparse direct solvers for the KKT system factorization

It is a good middle ground for problems in the 50–500 variable range where the dense JAX IPM starts to slow down but you do not want the external dependency burden of Ipopt.

Summary: Backend Recommendations#

Scenario	Recommended backend	Why
LP / QP (any size)	auto-dispatch	Specialized Mehrotra PC solver
Small NLP (< 50 vars)	`"ipm"` (default)	JIT-compiled, no overhead
Medium NLP (50–200 vars)	`"ipm"` or `"sparse_ipm"`	Dense JIT or sparse hybrid
Large NLP (> 200 vars)	`"ipopt"`	Sparse MA27/MA57 factorization
Batch of same-shape NLPs	`"ipm"`	JIT cache reuse; `vmap` for LP/QP
MILP / MIQP	auto-dispatch	LP/QP IPM at each B&B node
MINLP (small)	`"ipm"`	Fast per-node NLP solves
MINLP (large subproblems)	`"ipopt"`	Sparse subproblem solves
Dependency-free deployment	`"ripopt"`	Pure Rust, no C libraries

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 lp_ipm
print(f"(a) 100 LPs solved in {time.perf_counter() - t0:.4f}s (auto-dispatch to LP IPM)")

(a) 100 LPs solved in 0.0380s (auto-dispatch to LP IPM)

/Users/jkitchin/Dropbox/projects/discopt/python/discopt/modeling/core.py:1903: UserWarning: Variables with very large or infinite bounds: x[0] (lb=0, ub=1e+20), x[1] (lb=0, ub=1e+20), x[2] (lb=0, ub=1e+20), x[3] (lb=0, ub=1e+20), x[4] (lb=0, ub=1e+20). NLP solvers may fail (NaN, iteration_limit) when bounds exceed ~1e15. Add tighter explicit bounds, e.g. m.continuous('x', lb=0, ub=1000).
  result = solve_model(

# (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 ipm     : 3.8593s  obj=3.0741

(b) n=100 with ripopt  : 2.8141s  obj=3.0741

(b) n=100 with ipopt   : 3.0426s  obj=3.0741

# (c) Small MINLP for deployment -- Prediction: _________
# Answer: ripopt -- pure Rust, no C dependencies.
# If ripopt is not available in this environment, JAX IPM is the fallback.
model = build_small_minlp()
t0 = time.perf_counter()
result = model.solve(nlp_solver="ipm")  # change to "ripopt" 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}")

(c) Small MINLP: 0.1239s  obj=0.0000  nodes=1

Exercise answers:

(a) Auto-dispatch to the specialized LP IPM. No nlp_solver argument needed – discopt detects the LP structure and routes to the Mehrotra predictor-corrector LP solver with JIT cache reuse.
(b) For 100 variables, "ipopt" is likely fastest due to sparse MA27 factorization. "sparse_ipm" is a good second choice if Ipopt is unavailable.
(c) "ripopt" – 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 always fastest for these classes.
Use the JAX IPM (default) for small-to-medium NLPs and batch workloads – JIT compilation and vmap vectorization provide significant performance advantages.
Switch to cyipopt/Ipopt for large, sparse NLPs – the MA27/MA57 sparse direct solvers scale far better than dense linear algebra.
Choose ripopt for dependency-free deployment – pure Rust with no C/Fortran toolchain required.
Consider sparse_ipm as a middle ground – JAX autodiff combined with SciPy sparse factorization.

For further reading on interior-point methods see [Nocedal and Wright, 2006], [Wright, 1997], and [Mehrotra, 1992]. The Ipopt implementation is detailed in [Wächter and Biegler, 2006], and the broader MINLP landscape is surveyed in [Belotti et al., 2013].