POUNCE: Pure-Rust Ipopt Port

POUNCE: Pure-Rust Ipopt Port#

This tutorial introduces POUNCE, a pure-Rust port of the Ipopt interior-point NLP solver. POUNCE ships as its own compiled Python package, built with maturin and exposed to Python through PyO3 bindings; discopt drives it through a cyipopt-compatible interface.

By the end of this notebook you will be able to:

  1. Solve NLPs using the nlp_solver="pounce" backend.

  2. Understand the architecture: Rust solver with Python evaluation callbacks.

  3. Configure POUNCE solver options.

  4. Compare POUNCE performance against the Ipopt backend, and see why nlp_solver="ipm" is a compatibility alias rather than a third engine.

  5. Use POUNCE as the NLP subsolver inside Branch-and-Bound for MINLPs.

When to choose POUNCE:

  • Deployment without C dependencies. Unlike Ipopt (which requires BLAS, LAPACK, and optionally HSL), POUNCE is a single compiled extension with no external C/Fortran dependencies.

  • Single-threaded performance. Rust’s zero-cost abstractions and lack of garbage collection yield competitive single-threaded speed.

  • No Python GIL concerns. The Rust solver releases the GIL during computation, making it safe to use in multi-threaded Python applications.

References: The IPM algorithm is based on the primal-dual interior-point framework described in [Nocedal and Wright, 2006], Chapter 19.

import time

import discopt.modeling as dm
import numpy as np

Architecture#

POUNCE follows a callback-based design common to NLP solvers like Ipopt [Wächter and Biegler, 2006]. The architecture has four layers:

Python (discopt.modeling)
  |  formulates Model, builds NLPEvaluator
  v
Python (discopt.solvers.nlp_pounce)
  |  maps evaluator callbacks onto pounce.Problem
  |  (shape-compatible with cyipopt.Problem, so the cyipopt
  |   callback adapter is reused unchanged)
  v
PyO3 bindings (the `pounce` package's compiled extension)
  |  wraps the Python evaluator as a POUNCE NLP problem
  v
Rust solver (POUNCE)
  |  primal-dual IPM with line search
  |  sparse symmetric LDL^T factorization (the pure-Rust FERAL backend)
  |  calls back into Python for f(x), grad(x), Hessian(x)

Key points:

  • The Rust solver owns the IPM iteration loop, KKT system assembly, and line search logic.

  • Objective, gradient, constraint, Jacobian, and Hessian evaluations are performed by the Python NLPEvaluator via GIL-acquiring callbacks.

  • No external C libraries (BLAS, LAPACK, HSL) are required — all linear algebra is implemented in Rust.

  • Variable and constraint bounds are passed once at initialization; the Rust solver caches them as owned Vec<f64> data.

  • Nothing on this path touches JAX. The evaluator is numpy, and the factorization is Rust, so there is no tracing or compilation step between calling solve() and the first IPM iteration.

Simple NLP Examples#

Let us start with two small problems to verify that POUNCE produces correct results.

# Unconstrained quadratic: min (x - 3)^2 + (y + 1)^2
# Known optimal: x* = 3, y* = -1, f* = 0

m = dm.Model("unconstrained_quad")
x = m.continuous("x", lb=-10, ub=10)
y = m.continuous("y", lb=-10, ub=10)

m.minimize((x - 3) ** 2 + (y + 1) ** 2)

result = m.solve(nlp_solver="pounce")

print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}")
print(f"x = {result.x['x']:.6f}  (expected 3.0)")
print(f"y = {result.x['y']:.6f}  (expected -1.0)")
print(f"Time:      {result.wall_time:.4f}s")

Status:    optimal
Objective: 0.000000
x = 3.000000  (expected 3.0)
y = -1.000000  (expected -1.0)
Time:      0.0333s

# Constrained: min x^2 + y^2  s.t. x + y >= 1
# Known optimal: x* = y* = 0.5, f* = 0.5

m2 = dm.Model("constrained_simple")
x = m2.continuous("x", lb=-5, ub=5)
y = m2.continuous("y", lb=-5, ub=5)

m2.minimize(x**2 + y**2)
m2.subject_to(x + y >= 1.0, name="sum_lb")

result2 = m2.solve(nlp_solver="pounce")

print(f"Status:    {result2.status}")
print(f"Objective: {result2.objective:.6f}  (expected 0.5)")
print(f"x = {result2.x['x']:.6f}  (expected 0.5)")
print(f"y = {result2.x['y']:.6f}  (expected 0.5)")

Status:    optimal
Objective: 0.500000  (expected 0.5)
x = 0.500000  (expected 0.5)
y = 0.500000  (expected 0.5)

Constrained Rosenbrock#

The Rosenbrock function is a classic nonlinear test problem with a narrow curved valley. We add a disk constraint to make it more challenging:

\[ \min_{x,y} \quad (1 - x)^2 + 100\,(y - x^2)^2 \qquad \text{s.t.} \quad x^2 + y^2 \le 2, \quad x \in [-2, 2],\; y \in [-2, 2]. \]

The unconstrained optimum is \((1, 1)\) with \(f^* = 0\). The disk constraint \(x^2 + y^2 \le 2\) is inactive at that point, so the constrained optimum should also be \((1, 1)\).

mr = dm.Model("rosenbrock")
x = mr.continuous("x", lb=-2, ub=2)
y = mr.continuous("y", lb=-2, ub=2)

mr.minimize((1 - x) ** 2 + 100 * (y - x**2) ** 2)
mr.subject_to(x**2 + y**2 <= 2.0, name="disk")

rr = mr.solve(nlp_solver="pounce")

print(f"Status:    {rr.status}")
print(f"Objective: {rr.objective:.8f}  (expected ~0)")
print(f"x = {rr.x['x']:.6f}  (expected 1.0)")
print(f"y = {rr.x['y']:.6f}  (expected 1.0)")
print(f"Disk constraint: x^2 + y^2 = {rr.x['x'] ** 2 + rr.x['y'] ** 2:.4f} <= 2.0")

Status:    optimal
Objective: 0.00000002  (expected ~0)
x = 0.999862  (expected 1.0)
y = 0.999724  (expected 1.0)
Disk constraint: x^2 + y^2 = 1.9992 <= 2.0

Configuration#

POUNCE inherits solver options from the solve() call:

Option

Default

Description

time_limit

3600.0

Maximum wall-clock seconds

gap_tolerance

1e-4

Relative optimality gap for termination

abs_gap_tolerance

None

Absolute gap, when a relative one is not meaningful

Warning

max_iter, tol and print_level are not solve() keyword arguments, and this table used to list them as though they were (with time_limit given as 300, where the real default is 3600.0). solve() rejects an option it does not know rather than swallowing it, so following the old table did not silently do nothing — it raised:

TypeError: solve() got an unexpected keyword argument 'tol'. Did you mean 'abs_tol'? Unknown solver options are rejected rather than silently ignored

which is the right behaviour, and the reason this was worth correcting rather than leaving as harmless prose. The cell below asserts it.

Those three names are real POUNCE options — they just do not arrive through solve(). They go in an options dict on the entry points that take one, e.g. discopt.ml.trainable.train(model, options={"max_iter": 3000, "tol": 1e-8}) (see neural_dae).

# Solve with a time limit and iteration cap
mc = dm.Model("config_demo")
x = mc.continuous("x", lb=-10, ub=10)
y = mc.continuous("y", lb=-10, ub=10)
mc.minimize((x - 2) ** 2 + (y - 3) ** 2)

rc = mc.solve(nlp_solver="pounce", time_limit=10)
print(f"Status: {rc.status}, Objective: {rc.objective:.6f}")
print(f"x = {rc.x['x']:.4f}, y = {rc.x['y']:.4f}")
print(f"Wall time: {rc.wall_time:.4f}s (limit was 10s)")

# The table above, asserted: an option `solve()` does not know is rejected, not
# ignored. This is what makes a wrong options table a loud error and not a silent
# no-op -- worth demonstrating once.
try:
    mc.solve(nlp_solver="pounce", tol=1e-8)
except TypeError as exc:
    print(f"\nRejected as documented: {str(exc).splitlines()[0][:96]}")
else:
    raise AssertionError("an unknown solve() option was accepted")

Status: optimal, Objective: 0.000000
x = 2.0000, y = 3.0000
Wall time: 0.0025s (limit was 10s)

Rejected as documented: solve() got an unexpected keyword argument 'tol'. Did you mean 'abs_tol'? Unknown solver options

Performance Comparison#

discopt exposes two NLP engines for continuous solves:

  • "pounce" — the pure-Rust IPM described above; the default.

  • "ipopt" — Ipopt via cyipopt (requires BLAS/LAPACK, optionally HSL).

nlp_solver="ipm" and nlp_solver="sparse_ipm" are still accepted, but they are compatibility aliases, not engines. Each named the pure-JAX interior-point method that used to sit behind them; that solver has been retired, and both selectors now resolve to POUNCE. The cell after the comparison demonstrates this directly rather than asking you to take it on faith.

The sweep below therefore compares POUNCE against Ipopt, and only those two. Each problem is solved once per backend with no warm-up run, so read the absolute numbers as an order of magnitude rather than a benchmark — the first solve in the sweep absorbs one-off setup cost. For scaling measurements with repeated timing, see benchmarks_by_class.ipynb.

The interior-point framework is described in [Wright, 1997].

def make_unconstrained_quad(n=2):
    """min sum_i (x_i - i)^2"""
    m = dm.Model(f"quad_{n}")
    x = m.continuous("x", shape=(n,), lb=-10, ub=10)
    m.minimize(dm.sum(lambda i: (x[i] - i) ** 2, over=range(n)))
    return m


def make_constrained(n=5):
    """min sum x_i^2  s.t. sum x_i >= n/2"""
    m = dm.Model(f"constr_{n}")
    x = m.continuous("x", shape=(n,), lb=-5, ub=5)
    m.minimize(dm.sum(lambda i: x[i] ** 2, over=range(n)))
    m.subject_to(dm.sum(x) >= n / 2.0)
    return m


def make_rosenbrock():
    """min (1-x)^2 + 100*(y-x^2)^2"""
    m = dm.Model("rosenbrock")
    x = m.continuous("x", lb=-5, ub=5)
    y = m.continuous("y", lb=-5, ub=5)
    m.minimize((1 - x) ** 2 + 100 * (y - x**2) ** 2)
    return m


def make_sum_of_squares(n=10):
    """min sum_i (x_i - sin(i))^2  s.t. sum x_i <= n"""
    m = dm.Model(f"sos_{n}")
    x = m.continuous("x", shape=(n,), lb=-5, ub=5)
    targets = [float(np.sin(i)) for i in range(n)]
    m.minimize(dm.sum(lambda i: (x[i] - targets[i]) ** 2, over=range(n)))
    m.subject_to(dm.sum(x) <= float(n))
    return m


def make_portfolio(n=8):
    """min w' Sigma w  s.t. sum w = 1, mu'w >= 0.08"""
    np.random.seed(42)
    L = np.random.randn(n, n) * 0.05
    Sigma = L @ L.T + 0.01 * np.eye(n)
    mu = np.linspace(0.03, 0.15, n)
    m = dm.Model(f"portfolio_{n}")
    w = m.continuous("w", shape=(n,), lb=0, ub=1)
    m.minimize(
        dm.sum(
            lambda i: dm.sum(
                lambda j: Sigma[i, j] * w[i] * w[j],
                over=range(n),
            ),
            over=range(n),
        )
    )
    m.subject_to(dm.sum(w) == 1.0)
    m.subject_to(dm.sum(lambda i: mu[i] * w[i], over=range(n)) >= 0.08)
    return m


problems = [
    ("Quadratic (n=2)", make_unconstrained_quad(2)),
    ("Constrained (n=5)", make_constrained(5)),
    ("Rosenbrock (n=2)", make_rosenbrock()),
    ("Sum-of-squares (n=10)", make_sum_of_squares(10)),
    ("Portfolio (n=8)", make_portfolio(8)),
]

print(f"Prepared {len(problems)} test problems.")

Prepared 5 test problems.

backends = ["pounce", "ipopt"]
results_table = []

for name, model in problems:
    row = {"Problem": name}
    for backend in backends:
        try:
            t0 = time.perf_counter()
            r = model.solve(nlp_solver=backend)
            elapsed = time.perf_counter() - t0
            row[f"{backend}_time"] = f"{elapsed:.4f}s"
            row[f"{backend}_obj"] = f"{r.objective:.6f}"
            row[f"{backend}_status"] = r.status
        except Exception as e:
            # Ipopt is an optional dependency; if cyipopt is missing the table
            # says so instead of pretending the backend ran.
            row[f"{backend}_time"] = "unavailable"
            row[f"{backend}_obj"] = str(e)[:30]
            row[f"{backend}_status"] = "unavailable"
    results_table.append(row)

# Display results
header = f"{'Problem':<25s} | {'pounce':>12s} | {'ipopt':>12s}"
print(header)
print("-" * len(header))
for row in results_table:
    print(f"{row['Problem']:<25s} | {row['pounce_time']:>12s} | {row['ipopt_time']:>12s}")

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

Problem                   |       pounce |        ipopt
-------------------------------------------------------
Quadratic (n=2)           |      0.0028s |      0.0021s
Constrained (n=5)         |      0.0030s |      0.0026s
Rosenbrock (n=2)          |     14.1597s |     14.1705s
Sum-of-squares (n=10)     |      0.0039s |      0.0031s
Portfolio (n=8)           |      0.0266s |      0.0245s

# Verify both backends agree on objective values
print(f"{'Problem':<25s} | {'pounce obj':>14s} | {'ipopt obj':>14s}")
print("-" * 60)
for row in results_table:
    print(f"{row['Problem']:<25s} | {row['pounce_obj']:>14s} | {row['ipopt_obj']:>14s}")

Problem                   |     pounce obj |      ipopt obj
------------------------------------------------------------
Quadratic (n=2)           |       0.000000 |       0.000000
Constrained (n=5)         |       1.250000 |       1.250000
Rosenbrock (n=2)          |      -0.000000 |      -0.000000
Sum-of-squares (n=10)     |       0.000000 |       0.000000
Portfolio (n=8)           |       0.002064 |       0.002064

# `nlp_solver="ipm"` and `nlp_solver="sparse_ipm"` are not engines. Both named the
# retired pure-JAX IPM; both now resolve to POUNCE. Exact float equality below is
# what that looks like from the outside -- not "close", identical.
objectives = {}
print(f"{'selector':<12s} | {'status':<10s} | objective (repr, full precision)")
print("-" * 62)
for selector in ("pounce", "ipm", "sparse_ipm"):
    r = make_portfolio(8).solve(nlp_solver=selector)
    objectives[selector] = r.objective
    print(f"{selector:<12s} | {r.status:<10s} | {r.objective!r}")

print()
print("ipm == pounce         (exact):", objectives["ipm"] == objectives["pounce"])
print("sparse_ipm == pounce  (exact):", objectives["sparse_ipm"] == objectives["pounce"])

selector     | status     | objective (repr, full precision)
--------------------------------------------------------------
pounce       | optimal    | 0.0020639876561138674
ipm          | optimal    | 0.0020639876561138674
sparse_ipm   | optimal    | 0.0020639876561138674

ipm == pounce         (exact): True
sparse_ipm == pounce  (exact): True

Notes on timing:

  • There is no compilation step to amortize. The NLP path is numpy plus a compiled Rust solver, so nothing is traced or JIT-compiled between the solve() call and the first IPM iteration; a repeated solve of the same problem is not systematically cheaper than the first.

  • POUNCE times include the overhead of PyO3 callback crossings (Python to Rust and back for each evaluation). On problems this small that overhead, not the linear algebra, is most of the wall time.

  • Both solvers factor a sparse KKT system — Ipopt through MA57/MUMPS, POUNCE through its pure-Rust FERAL LDL^T backend. Ipopt’s linear algebra is the more mature of the two, and on large problems that maturity is where the difference shows up. Five problems of eight variables or fewer cannot tell you anything about that; see benchmarks_by_class.ipynb for problems at a size where the comparison means something.

Using POUNCE in Branch-and-Bound#

When solving a mixed-integer nonlinear program (MINLP), discopt uses Branch-and-Bound (B&B). At each tree node, a continuous NLP relaxation is solved. By setting nlp_solver="pounce", every node relaxation is solved with the Rust IPM (POUNCE), the default backend for single continuous solves.

# Small MINLP: min x^2 + y  s.t. x + y >= 2, y in {0, 1}
mm = dm.Model("minlp_pounce")
x = mm.continuous("x", lb=0, ub=5)
y = mm.binary("y")

mm.minimize(x**2 + y)
mm.subject_to(x + y >= 2.0, name="coupling")

rm = mm.solve(nlp_solver="pounce")

print(f"Status:     {rm.status}")
print(f"Objective:  {rm.objective:.6f}")
print(f"x = {rm.x['x']:.6f}")
print(f"y = {rm.x['y']:.0f}")
print(f"Node count: {rm.node_count}")
print(f"Gap:        {rm.gap:.2e}")
print(f"Wall time:  {rm.wall_time:.4f}s")

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

Status:     optimal
Objective:  2.000000
x = 1.000000
y = 1
Node count: 0
Gap:        0.00e+00
Wall time:  0.0048s

# Slightly larger MINLP with 3 binary variables
ml = dm.Model("minlp_larger")
x = ml.continuous("x", shape=(3,), lb=0, ub=5)
z = ml.binary("z", shape=(3,))

# Nonlinear objective
ml.minimize(dm.sum(lambda i: (x[i] - float(i + 1)) ** 2, over=range(3)) + dm.sum(z))

# Linking constraints: x[i] <= 5 * z[i]  (x active only if z selected)
for i in range(3):
    ml.subject_to(x[i] <= 5.0 * z[i], name=f"link_{i}")

# Budget: at most 2 active
ml.subject_to(dm.sum(z) <= 2, name="budget")

rl = ml.solve(nlp_solver="pounce")

print(f"Status:     {rl.status}")
print(f"Objective:  {rl.objective:.6f}")
print(f"x = {np.round(rl.x['x'], 4)}")
print(f"z = {np.round(rl.x['z'])}")
print(f"Node count: {rl.node_count}")
print(f"Wall time:  {rl.wall_time:.4f}s")

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

Status:     optimal
Objective:  3.000000
x = [0. 2. 3.]
z = [0. 1. 1.]
Node count: 0
Wall time:  0.0269s

When to Choose POUNCE#

Scenario

Recommended backend

Deployment in containers / minimal dependencies

POUNCE — no C toolchain needed

Batch solving many small NLPs

POUNCE — no compilation step to amortize

Large-scale NLPs (n > 1000)

ipopt — decades of tuning in HSL/MUMPS

Differentiable optimization / sensitivity

POUNCE — the default in discopt._relax.differentiable, which differentiates through the solve via the envelope theorem

Single-threaded, no-GIL performance

POUNCE — Rust execution

General-purpose robustness

ipopt — decades of development

Anything you cannot express as one of the two names above is not a third choice: "ipm" and "sparse_ipm" both land on POUNCE.

Drawbacks of POUNCE:

  • Evaluation overhead from PyO3 callbacks (each f(x), grad(x), etc. crosses the Python/Rust boundary).

  • Fewer algorithmic refinements than Ipopt (e.g., no adaptive barrier strategy, no inertia-free regularization).

  • Its sparse LDL^T backend is young next to HSL MA57, which is what the large-problem row of the table above is about.

Exercise#

Formulate a 10-variable constrained NLP and compare POUNCE to Ipopt. The problem:

\[ \min_{x} \quad \sum_{i=0}^{9} (x_i - t_i)^4 \qquad \text{s.t.} \quad \sum_{i=0}^{9} x_i = 5, \quad x_i \in [-3, 3]. \]

where \(t_i = \sin(i)\).

Both backends should report the same objective; the printed difference between them is the check.

n = 10
targets = [float(np.sin(i)) for i in range(n)]

# YOUR CODE HERE
# 1. Create a Model
# 2. Add continuous variables x with shape (n,), lb=-3, ub=3
# 3. Minimize sum_i (x[i] - targets[i])^4
# 4. Add equality constraint: sum(x) == 5
# 5. Solve with nlp_solver="pounce" and time it
# 6. Solve with nlp_solver="ipopt" and time it
# 7. Print both objectives and times, and check they agree

Solution#

Click to reveal
def build_quartic_model():
    m = dm.Model("quartic")
    x = m.continuous("x", shape=(n,), lb=-3, ub=3)
    m.minimize(dm.sum(lambda i: (x[i] - targets[i]) ** 4, over=range(n)))
    m.subject_to(dm.sum(x) == 5.0, name="sum_eq")
    return m


# Solve with POUNCE
m_pounce = build_quartic_model()
t0 = time.perf_counter()
r_pounce = m_pounce.solve(nlp_solver="pounce")
t_pounce = time.perf_counter() - t0

# Solve with Ipopt
m_ipopt = build_quartic_model()
t0 = time.perf_counter()
r_ipopt = m_ipopt.solve(nlp_solver="ipopt")
t_ipopt = time.perf_counter() - t0

print(f"{'Backend':<10s} | {'Status':<10s} | {'Objective':>12s} | {'Time':>8s}")
print("-" * 48)
print(f"{'pounce':<10s} | {r_pounce.status:<10s} | {r_pounce.objective:12.6f} | {t_pounce:8.4f}s")
print(f"{'ipopt':<10s} | {r_ipopt.status:<10s} | {r_ipopt.objective:12.6f} | {t_ipopt:8.4f}s")
print(f"\nObjective difference: {abs(r_pounce.objective - r_ipopt.objective):.2e}")

Backend    | Status     |    Objective |     Time
------------------------------------------------
pounce     | optimal    |     0.085947 |   0.0058s
ipopt      | optimal    |     0.085947 |   0.0047s

Objective difference: 4.16e-17

Summary#

In this tutorial you learned to:

  • Select the POUNCE backend with m.solve(nlp_solver="pounce").

  • Understand the Rust IPM architecture: a compiled pounce extension, PyO3 bindings, Python evaluation callbacks.

  • Configure solver options (time limit, tolerance, verbosity).

  • Compare POUNCE against Ipopt on small-to-medium NLPs.

  • Recognize "ipm" and "sparse_ipm" as compatibility aliases that resolve to POUNCE — shown above by bit-identical objectives, not asserted.

  • Use POUNCE as the NLP subsolver inside Branch-and-Bound for MINLPs.

POUNCE is the right choice when you need a dependency-free NLP solver for deployment, or when you want single-threaded Rust performance without the GIL. For large-scale problems, Ipopt’s linear algebra remains the more mature option [Wächter and Biegler, 2006]. Both implement the same primal-dual interior-point framework [Nocedal and Wright, 2006], which is why they agree so closely on the problems above.

See also: