Benchmarks by Problem Class

Benchmarks by Problem Class#

This notebook systematically benchmarks discopt against standard solvers for each problem class. For each class we:

Generate random problems of increasing size
Measure median-of-3 solve time
Verify solution accuracy
Visualize scaling behavior

LP comparisons use HiGHS [Huangfu and Hall, 2018] and scipy’s linprog. The LP IPM uses a Mehrotra predictor-corrector [Mehrotra, 1992]. MILP uses Branch & Bound [Land and Doig, 1960] with LP relaxations. NLP uses interior-point methods [Nocedal and Wright, 2006].

Class	discopt Solver	Comparison Solvers	Sizes
LP	LP IPM (auto)	HiGHS, scipy linprog	n=5,10,20,50,100
QP	QP IPM (auto)	scipy SLSQP	n=5,10,20,50,100
MILP	B&B + LP (auto)	HiGHS MIP	5-30 vars
MIQP	B&B + QP (auto)	Forced NLP path	5-12 vars
NLP	IPM / Ipopt	3 backends compared	2-20 vars
MINLP	B&B + NLP	batch_size comparison	5-12 vars

import time
import warnings

import discopt.modeling as dm
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import linprog, minimize

warnings.filterwarnings("ignore")

def median_time(fn, n_trials=3):
    """Run fn n_trials times, return (median_time, last_result)."""
    times = []
    result = None
    for _ in range(n_trials):
        t0 = time.perf_counter()
        result = fn()
        times.append(time.perf_counter() - t0)
    return np.median(times), result

print("Setup complete")

Setup complete

1. LP Benchmarks#

Random LP: min c’x s.t. Ax <= b, x >= 0. We compare discopt’s auto-detected LP IPM [Mehrotra, 1992] against HiGHS [Huangfu and Hall, 2018] (simplex) and scipy’s linprog.

import highspy

lp_sizes = [5, 10, 20, 50, 100]
lp_results = {"discopt": [], "highs": [], "scipy": []}

for n in lp_sizes:
    rng = np.random.RandomState(42 + n)
    m_cons = max(n // 2, 3)

    c_vec = rng.randn(n)
    A_mat = rng.randn(m_cons, n)
    b_vec = np.abs(A_mat @ np.ones(n)) + 1.0

    # --- discopt ---
    def run_discopt():
        m = dm.Model("lp_bench")
        x = m.continuous("x", shape=(n,), lb=0, ub=100)
        m.minimize(dm.sum(lambda i: c_vec[i] * x[i], over=range(n)))
        for k in range(m_cons):
            m.subject_to(dm.sum(lambda i: A_mat[k, i] * x[i], over=range(n)) <= b_vec[k])
        return m.solve()

    t_disc, r_disc = median_time(run_discopt)
    lp_results["discopt"].append(t_disc)

    # --- HiGHS ---
    # HiGHS is called directly through highspy. discopt used to ship a
    # `solvers.lp_highs` wrapper; today a pure LP passed to `m.solve()` is
    # itself routed to HiGHS with discopt-verified certificates (#1229), so this
    # baseline -- the bare library, called the same way as in the MILP cell
    # below -- measures discopt's modeling and certificate overhead.
    def run_highs_lp():
        h = highspy.Highs()
        h.setOptionValue("output_flag", False)
        for _ in range(n):
            h.addVar(0, 100)
        h.changeObjectiveSense(highspy.ObjSense.kMinimize)
        for i in range(n):
            h.changeColCost(i, c_vec[i])
        for k in range(m_cons):
            h.addRow(-highspy.kHighsInf, b_vec[k], n, list(range(n)), A_mat[k].tolist())
        h.run()
        return h.getInfoValue("objective_function_value")[1]

    t_highs, obj_highs = median_time(run_highs_lp)
    lp_results["highs"].append(t_highs)

    # --- scipy ---
    scipy_bounds = [(0, 100)] * n
    t_scipy, r_scipy = median_time(
        lambda: linprog(c_vec, A_ub=A_mat, b_ub=b_vec, bounds=scipy_bounds)
    )
    lp_results["scipy"].append(t_scipy)

    obj_disc = r_disc.objective if r_disc.objective is not None else float("nan")
    obj_scipy = r_scipy.fun if hasattr(r_scipy, "fun") else float("nan")
    print(
        f"n={n:3d} | discopt: {t_disc:.4f}s (obj={obj_disc:.4f}) | "
        f"HiGHS: {t_highs:.4f}s (obj={obj_highs:.4f}) | "
        f"scipy: {t_scipy:.4f}s (obj={obj_scipy:.4f})"
    )

n=  5 | discopt: 0.0186s (obj=-121.2780) | HiGHS: 0.0001s (obj=-121.2780) | scipy: 0.0006s (obj=-121.2780)

n= 10 | discopt: 0.0094s (obj=-20.5493) | HiGHS: 0.0001s (obj=-20.5493) | scipy: 0.0015s (obj=-20.5493)
n= 20 | discopt: 0.0142s (obj=-357.2774) | HiGHS: 0.0002s (obj=-357.2774) | scipy: 0.0005s (obj=-357.2774)

n= 50 | discopt: 0.1233s (obj=-1694.1677) | HiGHS: 0.0010s (obj=-1694.1677) | scipy: 0.0021s (obj=-1694.1677)

n=100 | discopt: 0.2644s (obj=-2262.3760) | HiGHS: 0.0016s (obj=-2262.3760) | scipy: 0.0021s (obj=-2262.3760)

2. QP Benchmarks#

Random QP: min 0.5 x’Qx + c’x s.t. sum(x) = 1, x >= 0 (portfolio-like). Compare discopt’s QP IPM against scipy SLSQP.

qp_sizes = [5, 10, 20, 50, 100]
qp_results = {"discopt": [], "scipy": []}

for n in qp_sizes:
    rng = np.random.RandomState(42 + n)
    # Random PSD matrix for covariance
    L = rng.randn(n, n) * 0.1
    Q = L.T @ L + 0.01 * np.eye(n)
    c_vec = rng.randn(n) * 0.01

    # --- discopt ---
    def run_discopt_qp():
        m = dm.Model("qp_bench")
        x = m.continuous("x", shape=(n,), lb=0, ub=1)
        obj = dm.sum(
            lambda i: dm.sum(lambda j: 0.5 * Q[i, j] * x[i] * x[j], over=range(n)), over=range(n)
        )
        obj = obj + dm.sum(lambda i: c_vec[i] * x[i], over=range(n))
        m.minimize(obj)
        m.subject_to(dm.sum(x) == 1)
        return m.solve()

    t_disc, r_disc = median_time(run_discopt_qp)
    qp_results["discopt"].append(t_disc)

    # --- scipy SLSQP ---
    def run_scipy_qp():
        def obj_fn(x):
            return 0.5 * x @ Q @ x + c_vec @ x

        x0 = np.ones(n) / n
        cons = {"type": "eq", "fun": lambda x: np.sum(x) - 1}
        bounds = [(0, 1)] * n
        return minimize(obj_fn, x0, method="SLSQP", bounds=bounds, constraints=cons)

    t_scipy, r_scipy = median_time(run_scipy_qp)
    qp_results["scipy"].append(t_scipy)

    obj_disc = r_disc.objective if r_disc.objective is not None else float("nan")
    obj_scipy = r_scipy.fun if hasattr(r_scipy, "fun") else float("nan")
    print(
        f"n={n:3d} | discopt: {t_disc:.4f}s (obj={obj_disc:.6f}) | "
        f"scipy: {t_scipy:.4f}s (obj={obj_scipy:.6f})"
    )

n=  5 | discopt: 0.0116s (obj=-0.009683) | scipy: 0.0014s (obj=-0.009683)
n= 10 | discopt: 0.0321s (obj=-0.001252) | scipy: 0.0021s (obj=-0.001251)

n= 20 | discopt: 0.1059s (obj=-0.001864) | scipy: 0.0046s (obj=-0.001862)

n= 50 | discopt: 0.6273s (obj=-0.005861) | scipy: 0.0112s (obj=-0.005860)

n=100 | discopt: 1.3726s (obj=-0.003202) | scipy: 0.0240s (obj=-0.003202)

3. MILP Benchmarks#

Random MILP: min c’x s.t. Ax <= b, x_int in {0,…,ub}, x_cont >= 0. Compare discopt’s B&B+LP [Land and Doig, 1960] against HiGHS MIP [Huangfu and Hall, 2018].

import highspy

milp_sizes = [5, 8, 12, 18, 25]
milp_results = {"discopt": [], "highs": []}

for n in milp_sizes:
    rng = np.random.RandomState(42 + n)
    n_int = max(n // 3, 2)
    n_cont = n - n_int
    m_cons = max(n // 2, 3)

    c_vec = rng.randn(n)
    A_mat = rng.randn(m_cons, n)
    b_vec = np.abs(A_mat @ np.ones(n)) + 2.0

    # --- discopt ---
    def run_discopt_milp():
        m = dm.Model("milp_bench")
        x_int = m.integer("x_int", shape=(n_int,), lb=0, ub=5)
        x_cont = m.continuous("x_cont", shape=(n_cont,), lb=0, ub=10)
        obj = dm.sum(lambda i: c_vec[i] * x_int[i], over=range(n_int)) + dm.sum(
            lambda i: c_vec[n_int + i] * x_cont[i], over=range(n_cont)
        )
        m.minimize(obj)
        for k in range(m_cons):
            lhs = dm.sum(lambda i: A_mat[k, i] * x_int[i], over=range(n_int)) + dm.sum(
                lambda i: A_mat[k, n_int + i] * x_cont[i], over=range(n_cont)
            )
            m.subject_to(lhs <= b_vec[k])
        return m.solve(time_limit=30, max_nodes=5000)

    t_disc, r_disc = median_time(run_discopt_milp)
    milp_results["discopt"].append(t_disc)

    # --- HiGHS MIP ---
    def run_highs_milp():
        h = highspy.Highs()
        h.setOptionValue("output_flag", False)
        h.setOptionValue("time_limit", 30.0)
        for i in range(n_int):
            h.addVar(0, 5)
            h.changeColIntegrality(i, highspy.HighsVarType.kInteger)
        for i in range(n_cont):
            h.addVar(0, 10)
        h.changeObjectiveSense(highspy.ObjSense.kMinimize)
        for i in range(n):
            h.changeColCost(i, c_vec[i])
        for k in range(m_cons):
            indices = list(range(n))
            values = A_mat[k].tolist()
            h.addRow(-highspy.kHighsInf, b_vec[k], len(indices), indices, values)
        h.run()
        info = h.getInfoValue("objective_function_value")
        return h.getModelStatus(), info[1]

    t_highs, r_highs = median_time(run_highs_milp)
    milp_results["highs"].append(t_highs)

    obj_disc = r_disc.objective if r_disc.objective is not None else float("nan")
    obj_highs = r_highs[1] if r_highs[0] == highspy.HighsModelStatus.kOptimal else float("nan")
    print(
        f"n={n:3d} (int={n_int}) | discopt: {t_disc:.4f}s (obj={obj_disc:.4f}, "
        f"status={r_disc.status}) | HiGHS: {t_highs:.4f}s (obj={obj_highs:.4f})"
    )

n=  5 (int=2) | discopt: 0.0044s (obj=-13.7409, status=optimal) | HiGHS: 0.0041s (obj=-13.7409)
n=  8 (int=2) | discopt: 0.0094s (obj=-16.3569, status=optimal) | HiGHS: 0.0015s (obj=-16.3569)
n= 12 (int=4) | discopt: 0.0088s (obj=-51.0009, status=optimal) | HiGHS: 0.0022s (obj=-51.0009)
n= 18 (int=6) | discopt: 0.0131s (obj=-57.0748, status=optimal) | HiGHS: 0.0018s (obj=-57.0748)

n= 25 (int=8) | discopt: 0.0181s (obj=-95.6228, status=optimal) | HiGHS: 0.0017s (obj=-95.6228)

4. MIQP Benchmarks#

Compare the specialized MIQP path (B&B + QP relaxation) against forcing the problem through the general NLP path.

Note

The cell below prints OA: nonlinear objective is not convex in the optimization sense; disabling master lower-bound updates and skipping objective OA cuts. The objective is convex; the outer-approximation route’s own convexity classifier false-negatives on a quadratic written as a nested dm.sum (#1352). It costs bound quality, not soundness — declaring the objective non-convex makes OA decline to update its master lower bound, which is the conservative direction — so the answer below is correct and certified. The line is left in rather than filtered, so the cost stays on the record. qp_solver works through it in full.

miqp_sizes = [4, 6, 8, 10]
miqp_results = {"miqp_path": [], "nlp_path": []}

for n in miqp_sizes:
    rng = np.random.RandomState(42 + n)
    n_int = max(n // 2, 2)
    n_cont = n - n_int

    L = rng.randn(n_cont, n_cont) * 0.1
    Q = L.T @ L + 0.01 * np.eye(n_cont)

    # --- MIQP path (auto-detected) ---
    def run_miqp():
        m = dm.Model("miqp_bench")
        x = m.continuous("x", shape=(n_cont,), lb=0, ub=1)
        z = m.binary("z", shape=(n_int,))
        obj = dm.sum(
            lambda i: dm.sum(lambda j: Q[i, j] * x[i] * x[j], over=range(n_cont)),
            over=range(n_cont),
        )
        obj = obj + dm.sum(lambda i: 0.5 * z[i], over=range(n_int))
        m.minimize(obj)
        m.subject_to(dm.sum(x) == 1)
        for i in range(min(n_cont, n_int)):
            m.subject_to(x[i] <= z[i])
        return m.solve(time_limit=30, max_nodes=1000)

    t_miqp, r_miqp = median_time(run_miqp)
    miqp_results["miqp_path"].append(t_miqp)

    # --- Forced NLP path (add dummy exp term) ---
    def run_nlp_forced():
        m = dm.Model("miqp_nlp")
        x = m.continuous("x", shape=(n_cont,), lb=0, ub=1)
        z = m.binary("z", shape=(n_int,))
        obj = dm.sum(
            lambda i: dm.sum(lambda j: Q[i, j] * x[i] * x[j], over=range(n_cont)),
            over=range(n_cont),
        )
        obj = obj + dm.sum(lambda i: 0.5 * z[i], over=range(n_int))
        obj = obj + 0 * dm.exp(x[0])  # force MINLP
        m.minimize(obj)
        m.subject_to(dm.sum(x) == 1)
        for i in range(min(n_cont, n_int)):
            m.subject_to(x[i] <= z[i])
        return m.solve(time_limit=30, max_nodes=1000)

    t_nlp, r_nlp = median_time(run_nlp_forced)
    miqp_results["nlp_path"].append(t_nlp)

    obj_miqp = r_miqp.objective if r_miqp.objective is not None else float("nan")
    obj_nlp = r_nlp.objective if r_nlp.objective is not None else float("nan")
    speedup = t_nlp / t_miqp if t_miqp > 0 else float("nan")
    print(
        f"n={n:2d} | MIQP path: {t_miqp:.3f}s (obj={obj_miqp:.4f}) | "
        f"NLP path: {t_nlp:.3f}s (obj={obj_nlp:.4f}) | speedup: {speedup:.1f}x"
    )

n= 4 | MIQP path: 0.012s (obj=0.5202) | NLP path: 0.044s (obj=0.5202) | speedup: 3.7x

n= 6 | MIQP path: 0.016s (obj=0.5265) | NLP path: 0.202s (obj=0.5265) | speedup: 12.5x

n= 8 | MIQP path: 0.021s (obj=0.5301) | NLP path: 0.748s (obj=0.5301) | speedup: 36.2x

n=10 | MIQP path: 0.030s (obj=0.5191) | NLP path: 1.097s (obj=0.5191) | speedup: 36.7x

5. NLP Benchmarks#

Compare discopt’s NLP backends on Rosenbrock-like problems of increasing dimension. The nlp_solver="ipm" selection below is a back-compat alias: the pure-JAX interior-point method has been retired, and on the NLP/MINLP path "ipm" now resolves to POUNCE, the pure-Rust Ipopt port. So this cell compares POUNCE against cyipopt [Wächter and Biegler, 2006] — two implementations of the same primal-dual interior-point algorithm.

nlp_sizes = [2, 5, 10, 15, 20]
nlp_results = {"pounce": [], "ipopt": []}

# Check which backends are available
has_ipopt = True
try:
    import cyipopt  # noqa: F401
except ImportError:
    has_ipopt = False
    print("cyipopt not available, skipping Ipopt comparison")

for n in nlp_sizes:
    # Rosenbrock: sum of 100*(x_{i+1} - x_i^2)^2 + (1 - x_i)^2
    def run_nlp_backend(backend):
        m = dm.Model("rosenbrock")
        x = m.continuous("x", shape=(n,), lb=-5, ub=5)
        obj = dm.sum(
            lambda i: 100 * (x[i + 1] - x[i] ** 2) ** 2 + (1 - x[i]) ** 2,
            over=range(n - 1),
        )
        m.minimize(obj)
        return m.solve(nlp_solver=backend)

    # IPM
    # "pounce" spelled out: `nlp_solver="ipm"` is an alias for it, so naming the
    # alias here made the column read as a third engine (see the note above).
    t_pounce, r_pounce = median_time(lambda: run_nlp_backend("pounce"))
    nlp_results["pounce"].append(t_pounce)

    # Ipopt
    if has_ipopt:
        t_ipopt, r_ipopt = median_time(lambda: run_nlp_backend("ipopt"))
        nlp_results["ipopt"].append(t_ipopt)
    else:
        nlp_results["ipopt"].append(float("nan"))

    obj_pounce = r_pounce.objective if r_pounce.objective is not None else float("nan")
    line = f"n={n:2d} | POUNCE: {t_pounce:.4f}s (obj={obj_pounce:.6f})"
    if has_ipopt:
        obj_ipopt = r_ipopt.objective if r_ipopt.objective is not None else float("nan")
        line += f" | Ipopt: {t_ipopt:.4f}s (obj={obj_ipopt:.6f})"
        # Two implementations of the same interior-point method on a convex-at-the
        # -solution problem: if they disagree, one of them is wrong and the timing
        # comparison is meaningless. Rosenbrock's optimum is the all-ones point.
        assert abs(obj_pounce - obj_ipopt) < 1e-4 * (1 + abs(obj_pounce)), (
            f"n={n}: POUNCE {obj_pounce} vs Ipopt {obj_ipopt}"
        )
    print(line)

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

n= 2 | POUNCE: 0.0300s (obj=0.000000) | Ipopt: 0.0309s (obj=0.000000)

n= 5 | POUNCE: 0.0496s (obj=0.000000) | Ipopt: 0.0495s (obj=0.000000)

n=10 | POUNCE: 0.0763s (obj=0.000000) | Ipopt: 0.0824s (obj=0.000000)

n=15 | POUNCE: 0.1198s (obj=0.000000) | Ipopt: 0.1225s (obj=0.000000)

n=20 | POUNCE: 0.1648s (obj=0.000000) | Ipopt: 0.1657s (obj=0.000000)

6. MINLP Benchmarks#

Self-comparison: how does batch_size affect MINLP solve time? We solve the same problem with batch_size=1, 8, and 16. batch_size is the number of B&B nodes exported from the tree per iteration [Belotti et al., 2013]; larger batches amortize the per-iteration marshaling cost across more nodes, so any effect shows up only when the tree has enough open nodes to fill a batch.

minlp_sizes = [4, 6, 8, 10]
batch_sizes = [1, 8, 16]
minlp_results = {f"batch={bs}": [] for bs in batch_sizes}

for n in minlp_sizes:
    rng = np.random.RandomState(42 + n)
    n_int = max(n // 3, 2)
    n_cont = n - n_int

    line_parts = [f"n={n:2d}"]
    for bs in batch_sizes:

        def run_minlp():
            m = dm.Model("minlp_bench")
            x = m.continuous("x", shape=(n_cont,), lb=0, ub=5)
            y = m.binary("y", shape=(n_int,))
            obj = dm.sum(lambda i: dm.exp(x[i] * 0.3) + x[i], over=range(n_cont)) + dm.sum(
                lambda i: 3.0 * y[i], over=range(n_int)
            )
            m.minimize(obj)
            m.subject_to(dm.sum(x) >= 2)
            for i in range(min(n_cont, n_int)):
                m.subject_to(x[i] <= 5 * y[i])
            return m.solve(batch_size=bs, time_limit=30, max_nodes=500)

        t, r = median_time(run_minlp)
        minlp_results[f"batch={bs}"].append(t)
        obj_val = r.objective if r.objective is not None else float("nan")
        line_parts.append(f"batch={bs}: {t:.3f}s (obj={obj_val:.3f}, nodes={r.node_count})")

    print(" | ".join(line_parts))

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

n= 4 | batch=1: 0.069s (obj=7.822, nodes=5) | batch=8: 0.064s (obj=7.822, nodes=5) | batch=16: 0.020s (obj=7.822, nodes=0)
n= 6 | batch=1: 0.025s (obj=6.700, nodes=1) | batch=8: 0.024s (obj=6.700, nodes=1) | batch=16: 0.011s (obj=6.700, nodes=0)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 3 rows (3 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

OA: generating OA cuts for 1 of 4 rows (4 convex; the rest are either nonconvex or already exact in the master)

n= 8 | batch=1: 0.025s (obj=8.647, nodes=1) | batch=8: 0.024s (obj=8.647, nodes=1) | batch=16: 0.012s (obj=8.647, nodes=0)

n=10 | batch=1: 0.029s (obj=9.647, nodes=1) | batch=8: 0.027s (obj=9.647, nodes=1) | batch=16: 0.012s (obj=9.647, nodes=0)

Scaling Visualizations#

A 2x3 grid showing solve time vs problem size for each problem class.

fig, axes = plt.subplots(2, 3, figsize=(15, 9))
fig.suptitle("discopt Solver Scaling by Problem Class", fontsize=14, fontweight="bold")

# LP
ax = axes[0, 0]
ax.semilogy(lp_sizes, lp_results["discopt"], "o-", label="discopt (HiGHS route)", color="#2196F3")
ax.semilogy(lp_sizes, lp_results["highs"], "s--", label="HiGHS", color="#FF9800")
ax.semilogy(lp_sizes, lp_results["scipy"], "^--", label="scipy linprog", color="#4CAF50")
ax.set_xlabel("Variables (n)")
ax.set_ylabel("Time (s)")
ax.set_title("LP")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# QP
ax = axes[0, 1]
ax.semilogy(qp_sizes, qp_results["discopt"], "o-", label="discopt (POUNCE QP path)", color="#2196F3")
ax.semilogy(qp_sizes, qp_results["scipy"], "s--", label="scipy SLSQP", color="#4CAF50")
ax.set_xlabel("Variables (n)")
ax.set_ylabel("Time (s)")
ax.set_title("QP")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# MILP
ax = axes[0, 2]
ax.semilogy(milp_sizes, milp_results["discopt"], "o-", label="discopt B&B+LP", color="#2196F3")
ax.semilogy(milp_sizes, milp_results["highs"], "s--", label="HiGHS MIP", color="#FF9800")
ax.set_xlabel("Variables (n)")
ax.set_ylabel("Time (s)")
ax.set_title("MILP")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# MIQP
ax = axes[1, 0]
ax.plot(miqp_sizes, miqp_results["miqp_path"], "o-", label="MIQP path (B&B+QP)", color="#2196F3")
ax.plot(miqp_sizes, miqp_results["nlp_path"], "s--", label="Forced NLP path", color="#F44336")
ax.set_xlabel("Variables (n)")
ax.set_ylabel("Time (s)")
ax.set_title("MIQP")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# NLP
ax = axes[1, 1]
ax.semilogy(nlp_sizes, nlp_results["pounce"], "o-", label="discopt POUNCE", color="#2196F3")
if has_ipopt:
    ax.semilogy(nlp_sizes, nlp_results["ipopt"], "s--", label="Ipopt", color="#FF9800")
ax.set_xlabel("Variables (n)")
ax.set_ylabel("Time (s)")
ax.set_title("NLP (Rosenbrock)")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# MINLP
ax = axes[1, 2]
colors_bs = ["#2196F3", "#FF9800", "#4CAF50"]
for idx, bs in enumerate(batch_sizes):
    ax.plot(
        minlp_sizes,
        minlp_results[f"batch={bs}"],
        "o-",
        label=f"batch_size={bs}",
        color=colors_bs[idx],
    )
ax.set_xlabel("Variables (n)")
ax.set_ylabel("Time (s)")
ax.set_title("MINLP (batch size comparison)")
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("benchmarks_by_class.png", dpi=150, bbox_inches="tight")
plt.show()
print("Plot saved to benchmarks_by_class.png")

../_images/9dc86c3bfe38bfd34200843313ac8a6fc8de4b2f77383e75f141546c7fac7a88.png

Plot saved to benchmarks_by_class.png

Summary#

Key observations:

LP/QP: discopt’s specialized IPM is competitive with HiGHS [Huangfu and Hall, 2018] and scipy for small-to-medium problems. HiGHS (C++ simplex) is faster for large LPs, as expected for a mature production solver.
MILP: HiGHS MIP is significantly faster due to presolve, cutting planes, and decades of engineering. discopt’s B&B+LP path [Land and Doig, 1960] is correct but not yet competitive on large instances.
MIQP: The specialized MIQP path (B&B + QP relaxation) offers significant speedup over forcing the problem through the general NLP path.
NLP: POUNCE (reached here through the "ipm" back-compat alias) and cyipopt [Wächter and Biegler, 2006] produce equivalent solutions, as they should — POUNCE is a Rust port of the same interior-point algorithm — and their timings are within noise of each other at these sizes.
MINLP: batch_size has little effect on these instances because the B&B tree closes in 1–5 nodes, so there are never enough open nodes to fill a batch [Belotti et al., 2013]. The knob matters only on trees deep enough to keep a batch populated.