Introduction

POUNCE is a general interior-point method, implemented in pure Rust — one numerical backbone that now spans nonlinear, conic/quadratic, and polynomial global optimization rather than a single problem class. Its nonlinear-programming core began as a faithful port of the Ipopt filter line-search method — the algorithm, console output, and option semantics follow upstream Ipopt closely enough that anyone used to reading ipopt logs can drop in pounce without relearning where the numbers live — and it has since grown into a family of solvers sharing that backbone:

• Nonlinear programming — the filter line-search interior-point method (the Ipopt port) plus an active-set SQP path, for general smooth problems

  min  f(x)
  s.t. g_L <= g(x) <= g_U
       x_L <=   x  <= x_U
  

  where f and g are twice-continuously-differentiable.

• Conic & quadratic — LP, convex QP, second-order (SOCP), positive-semidefinite (SDP), and the non-symmetric exponential and power cones, each solved to the global optimum.

• Global optimization — certified global optima for nonconvex polynomial problems via SOS / Lasserre relaxations. Nonconvex problems that are not polynomial get a local answer; POUNCE has no spatial branch-and-bound solver.

See Choosing a Solver for which solver fits which problem.

Pure Rust by default

The default build is pure Rust — no Fortran, no commercial solver, no system BLAS required. The bundled FERAL backend provides a sparse symmetric LDLᵀ factorization. The HSL MA57 backend is available behind the optional ma57 feature for users who have a license for libcoinhsl and have it installed (see Installation).

Status

Production-ready for the core IPM workflow. The algorithm-side core, NLP interface, line search, filter, barrier update (monotone + Mehrotra adaptive), KKT solve, restoration phase, AMPL .nl reader, the C ABI (pounce-cinterface), the Python wrapper (pounce-solver), and the CLI all solve a wide range of NLPs from the standard test suites (Hock-Schittkowski, CUTEst, Mittelmann ampl-nlp, CHO parameter estimation, gas/water network design). Sensitivity analysis (sIPOPT port), reduced-Hessian computation, the auxiliary-equality + FBBT presolve, and the active-set SQP path are all wired in and available behind option keys. Existing PyIpopt / cyipopt / JuMP / AMPL clients link against libpounce_cinterface in place of libipopt unchanged.

The conic and global solvers are wired end-to-end alongside the NLP core: the convex interior-point solver (pounce-convex) handles LP / QP, SOCP, exponential / power cones, and small SDPs — with a Conic Benchmark Format (.cbf) reader cross-checked against the CBLIB tier — and adds SOS / Lasserre polynomial global optimization (sos_minimize). These are reachable from the CLI, the Python package, and the JSON solve report. There is no spatial branch-and-bound solver for general factorable nonconvex problems — outside the polynomial case, nonconvex models are solved locally.

License

EPL-2.0, the same license as upstream Ipopt.

Where to go next

• Installation — build and install POUNCE.
• Quick Start — solve your first problem.
• Running Solves — the command-line driver in depth.
• Acknowledgments — the papers behind the algorithm.

Installation

Three routes, in the order most people want them:

With pip

pip install pounce-solver

Prebuilt wheels for Linux, macOS, and Windows (CPython 3.9+). No Rust toolchain is involved. This installs both interfaces at once:

pounce problem.nl                       # the CLI
python -c "import pounce; print(pounce.__version__)"

For Pyomo models:

pip install pyomo-pounce

import pyomo.environ as pyo
results = pyo.SolverFactory("pounce").solve(model)

Optional extras — none needed for a normal solve:

pip install "pounce-solver[jax]"     # pounce.jax autodiff frontend
pip install "pounce-solver[torch]"   # pounce.torch autodiff frontend
pip install "pounce-solver[viz]"     # debugger plots (pounce-dbg-viz)
pip install "pounce-solver[gams]"    # GAMS solver link — see gams.md

If the CLI will not start (GLIBC_2.39 not found)

Releases up to and including 0.9.0 bundled a Linux CLI built against a newer glibc than the wheel advertised, so pounce fails to exec on older distributions (Debian 12, Ubuntu 22.04, RHEL/Rocky/Alma 8 and 9, and most HPC images) with:

pounce/bin/pounce: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.39' not found

import pounce still works; only the CLI (and therefore Pyomo, which shells out to it) is affected. This is fixed for releases after 0.9.0. If you hit it, use the container or build from source below.

With a container

No toolchain and nothing installed on the host:

docker run --rm -v "$PWD:/work" ghcr.io/jkitchin/pounce:latest problem.nl
apptainer pull pounce.sif docker://ghcr.io/jkitchin/pounce:0.9.0

Both images carry the CLI, the Python API, and the Pyomo plugin. See Docker & Containers for tags, bind mounts, and a Slurm example.

From source

Prerequisites

A stable Rust toolchain. Nothing else is needed for the default pure-Rust build. Install Rust via rustup:

curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh
source "$HOME/.cargo/env"

Verify the install:

rustc --version && cargo --version

Build

From the repository root:

make            # release build of the workspace
make test       # run all tests
make clippy     # lint
make doc        # rustdoc for the Rust API

Install

make install                          # installs to $HOME/.local
sudo make install PREFIX=/usr/local   # or system-wide

This drops the pounce binary into $PREFIX/bin and the libpounce_cinterface shared library into $PREFIX/lib. Make sure $HOME/.local/bin is on your PATH, then verify:

pounce --version

HSL MA57 backend (optional)

The default FERAL backend needs no external libraries. To build with the HSL MA57 linear solver instead, you need a CoinHSL install whose lib/ directory holds libcoinhsl. Point the COINHSL_DIR environment variable at it and build with the ma57 feature:

export COINHSL_DIR=/path/to/CoinHSL
cargo build -p pounce-cli --release --features ma57

The feature makes MA57 available; selecting it is a separate step, because linear_solver defaults to feral in every build:

pounce problem.nl linear_solver=ma57

Build CoinHSL from https://www.hsl.rl.ac.uk/ipopt/. MA57 is primarily useful for benchmarking against upstream Ipopt; the FERAL backend is the supported default for everyday use, and a build without --features ma57 never touches COINHSL_DIR.

Using POUNCE as a Rust library

The workspace is a set of library crates (see Algorithm & Workspace for the layout). To browse the Rust API, build and open the rustdoc:

make doc        # generates target/doc

Docker & Containers

Prebuilt images carry the whole POUNCE surface — the pounce CLI, the Python API, and the Pyomo plugin — with no Rust toolchain and no pip install on your side. This is the path of least resistance on a cluster, where you often cannot install a toolchain anyway, and the quickest way to try POUNCE without touching your environment.

docker pull ghcr.io/jkitchin/pounce:latest
docker run --rm ghcr.io/jkitchin/pounce:latest --version

Which tag

Pin X.Y.Z for anything you intend to re-run months later. latest and edge both move under you.

Running solves

The entrypoint is the pounce CLI, so arguments after the image name go straight to it:

docker run --rm ghcr.io/jkitchin/pounce:latest --list-problems
docker run --rm ghcr.io/jkitchin/pounce:latest --problem rosenbrock

Your own problems live on the host, so mount a directory. The working directory inside the image is /work:

docker run --rm -v "$PWD:/work" ghcr.io/jkitchin/pounce:latest \
  problem.nl print_level=5 tol=1e-10

That writes problem.sol next to problem.nl in the mounted directory — see Running Solves for the full option surface. The image runs as UID 1000 rather than root, so files it creates in a bind mount are not root-owned on the host. If your host UID differs, add --user "$(id -u):$(id -g)".

Python and Pyomo

Override the entrypoint to get a shell or an interpreter:

docker run --rm -it --entrypoint python ghcr.io/jkitchin/pounce:latest
docker run --rm -it --entrypoint bash -v "$PWD:/work" ghcr.io/jkitchin/pounce:latest

Both import pounce (with numpy and scipy) and Pyomo’s SolverFactory('pounce') work out of the box:

docker run --rm -v "$PWD:/work" --entrypoint python \
  ghcr.io/jkitchin/pounce:latest my_model.py

The optional extras are not installed — no JAX, no PyTorch, no plotly, no GAMS bindings. Add what you need in a derived image:

FROM ghcr.io/jkitchin/pounce:0.10.0
USER root
RUN pip install --no-cache-dir "pounce-solver[jax]==0.9.0"
USER pounce

On a cluster (Apptainer / Singularity)

Most HPC sites run Apptainer (formerly Singularity) rather than Docker, because it needs no daemon and no root. It pulls Docker images directly:

apptainer pull pounce.sif docker://ghcr.io/jkitchin/pounce:0.10.0
apptainer run pounce.sif problem.nl

Two differences from docker run are worth knowing before you write a job script:

• You are yourself. Apptainer runs the container as your own user, not the image’s, so output files land with your ownership and no --user flag is needed.
• $HOME and $PWD are already there. Apptainer bind-mounts them by default, so a .nl file in your submit directory is usually visible with no -B flag at all. Add -B /scratch:/scratch (or your site’s equivalent) for anything outside them.

Build the .sif once on a login node and reuse it — pulling on every array task hammers the registry and will get rate-limited. In a Slurm script:

#!/bin/bash
#SBATCH --job-name=pounce
#SBATCH --cpus-per-task=4

apptainer exec $HOME/images/pounce.sif \
  pounce $SLURM_SUBMIT_DIR/problem.nl --json-output result.json

Pin the digest rather than the tag if the run has to be reproducible years later — a tag can be re-pushed, a digest cannot:

apptainer pull pounce.sif docker://ghcr.io/jkitchin/pounce@sha256:<digest>

Building your own

You do not need the images to be published — both Dockerfiles are in the repository and build from a clone. From the repository root:

make docker          # compiles the current working tree -> pounce:dev
make docker-release  # installs the released wheels -> pounce:<version>

make docker is the one to reach for when testing a branch: it compiles whatever is checked out, and stamps the commit into pounce --about so the image can say what it contains. make docker-release needs no Rust toolchain and takes seconds. See docker/README.md for build arguments and the .dockerignore caveat.

What is not in the image

The HSL MA57 backend is absent, and cannot be added by us: CoinHSL is license-restricted and not redistributable. Passing linear_solver=ma57 will not work in a container. The pure-Rust FERAL backend is the default everywhere and needs no external libraries — see Installation if you have a CoinHSL license and want a local build with it.

The GAMS link is likewise absent, since it needs your own GAMS install and license on the host. See GAMS.

Quick Start

This page assumes POUNCE is built and on your PATH (see Installation).

Solve an AMPL .nl file

pounce problem.nl

This solves the problem and writes a sibling problem.sol next to the input, following the AMPL solver convention. The console output mirrors upstream ipopt’s banner, per-iteration table, and final summary.

Append KEY=VALUE pairs to override options — the syntax and semantics match the upstream Ipopt CLI:

pounce problem.nl print_level=8 max_iter=500 tol=1e-10

See Solver Options for details.

Try a built-in problem

POUNCE ships several self-contained test problems that exercise the full pipeline without parsing a .nl file (run pounce --list-problems for the full set):

pounce --list-problems
pounce --problem rosenbrock
pounce --problem quadratic

From Python

import numpy as np
from pounce import minimize

res = minimize(lambda x: ((x - 1) ** 2).sum(), x0=np.zeros(3))
print(res.fun, res.x)

See the Python API chapter for the full cyipopt-compatible interface.

From Pyomo

import pyomo_pounce  # registers 'pounce'
from pyomo.environ import SolverFactory

SolverFactory('pounce').solve(model)

See the Pyomo chapter for details.

Full help

pounce --help

Choosing a Solver

POUNCE is not a single solver but a small family of them sharing one numerical backbone. This page is the map: what each solver is, when to reach for it, and how they fit together.

The one-sentence version: convex and conic problems are solved to the global optimum; nonconvex problems are solved locally by default, or to a certified global optimum via the SOS (polynomial) and spatial branch-and-bound (general) paths. Every solver, whatever its flavor, ultimately factorizes a symmetric KKT system through the shared pounce-linsol layer, which in turn drives a pluggable backend (FERAL by default, HSL MA57 optionally).

The solvers at a glance

When to reach for the active-set QP. auto never selects it: a cold, one-shot convex QP goes to the interior-point path, which is materially more robust on that workload (137 of the 138 Maros-Mészáros problems, against substantially fewer for a cold active-set solve). That is the character of the method rather than a defect — an active-set iteration count is combinatorial in the size of the active set, while an interior-point count is nearly independent of problem size. Choose solver_selection=qp-active-set when you want an exact vertex solution, or when you are solving a sequence of similar QPs — MPC steps, branch-and-bound nodes, continuation — where the working set carries across solves and solve_parametric can trace the homotopy from the previous solution instead of starting over.

POUNCE has no spatial branch-and-bound solver for general factorable nonconvex NLPs — no solver_selection=global CLI route, no minimize_global Python entry point. The only certified-global path for a nonconvex problem is SOS / Lasserre, and it covers polynomials only.

When to choose each

General nonlinear program (the common case) → NLP filter-IPM

If your model has nonlinear objective or constraints and you don’t know (or can’t assume) convexity, this is the default and the most mature path. It is POUNCE’s port of Ipopt’s filter line-search interior-point method: robust on nonconvex problems, with a feasibility restoration phase for hard starts and exact or limited-memory Hessians. It returns a local KKT point — for a nonconvex problem there is no global guarantee.

• CLI: pounce model.nl (or a built-in problem).
• Python: the cyipopt-style Problem class, or the scipy-style minimize facade.
• Reach for limited-memory Hessians (hessian_approximation=limited-memory) when second derivatives are unavailable or expensive.

A sequence of related NLPs, or a stable active set → NLP active-set SQP

Selected with algorithm=active-set-sqp. It solves the NLP as a sequence of quadratic subproblems (handed to pounce-qp), which warm-starts extremely well when the active set is stable across solves — e.g. a parametric sweep or a control loop. For a single cold solve of a general NLP, prefer the filter-IPM.

Linear or convex quadratic program → Convex IPM (LP/QP)

If P ⪰ 0 (or P = 0 for an LP), use the convex interior-point solver: it returns the global optimum, detects primal/dual infeasibility, and offers warm-starting, batched and multiple-RHS solving, a build-once / solve-many QpFactorization handle, and post-optimal sensitivity (QpSensitivity — the sIPOPT analog). The CLI’s auto routing classifies an .nl and sends LP/convex-QP problems here automatically.

• Python: pounce.qp.solve_qp (and solve_qp_batch, solve_qp_multi_rhs).

Second-order, exponential, or power cones → Convex IPM (conic)

The same convex solver handles conic programs: second-order cones, the exponential and power cones that express geometric programming, entropy / log-sum-exp, logistic models, and p-norm constraints, and the positive-semidefinite cone for small dense SDPs. Also global. This is the path to use when you can cast a nominally-nonconvex problem into a convex cone — you trade modeling effort for a global guarantee. (The PSD cone is self-scaled and runs on the symmetric driver; the exp/power cones run on the non-symmetric HSDE driver, so the two families can’t yet be mixed in one problem.)

A common special case routes here automatically: a convex quadratically-constrained QP (QCQP). When auto routing finds a convex-quadratic inequality ½xᵀHx + aᵀx + b ≤ 0 (H ⪰ 0), it reformulates each such constraint to one second-order cone (H = FᵀF) and sends the whole problem to the conic solver — no .cbf and no manual cone bookkeeping needed. This works from a .nl/Pyomo model on the CLI and from minimize() in Python (which probes each constraint’s Hessian and only routes when it can prove the feasible set is convex). See LP / QP Solver Routing.

• Python: pounce.qp.solve_socp(..., cones=[("exp", 3), ("pow", 0.5), ...]) for an explicit cone program, or just minimize(...) for a convex QCQP.
• CLI: a Conic Benchmark Format file, pounce model.cbf (see the CBLIB benchmark tier), or any convex-QCQP .nl under auto routing.

Nonconvex problem, global optimum required → SOS (polynomials only)

When the problem is genuinely nonconvex and a local optimum is not good enough, the one path to a certified global optimum is for polynomials:

• Polynomial objective/constraints → SOS / Lasserre (sos_minimize, or pounce.sos_minimize). A single semidefinite program certifies the global minimum (the largest γ with p − γ in the Putinar cone), and the global minimizers are recovered from the moment matrix — even multiple ones, via a facial-reduction step. Best for modest degree and dimension; the SDP grows with the relaxation order.

If the problem is nonconvex and not polynomial (exp/ln/trig), POUNCE cannot certify a global optimum. Reformulate into the convex cone library if you can; otherwise multistart the local NLP solver and accept that the result is uncertified.

See Global Optimization for the SOS path in depth, and for the multistart fallback.

Indefinite QP, or a QP inner-solver → Active-set QP

pounce-qp is a sparse parametric active-set solver that accepts an indefinite Hessian (via inertia control), with two-sided bounds and factorization-reuse across a homotopy. It is the engine behind the active-set SQP path, and is the right choice for MPC-style problems or any setting where you re-solve a slowly-changing QP many times. Use the convex IPM instead when P ⪰ 0 and you want a single robust solve with infeasibility certificates.

How to override the automatic routing

The CLI classifies each .nl problem and picks a solver, but you can force the choice:

pounce model.nl solver_selection=auto          # default: classify, then route
pounce model.nl solver_selection=nlp           # filter-IPM (or active-set-sqp via algorithm=)
pounce model.nl solver_selection=lp-ipm        # convex LP interior-point
pounce model.nl solver_selection=qp-ipm        # convex QP interior-point
pounce model.nl solver_selection=socp          # conic interior-point (convex QCQP)
pounce model.nl solver_selection=qp-active-set # active-set QP

solver_selection is an ordinary POUNCE option, not a command-line flag: it is passed as a trailing KEY=VALUE pair (the ipopt CLI convention), and so also works from an options file, the pounce_options environment variable, or Pyomo’s solver.options. Forcing a value the problem class does not support is rejected with a message rather than silently ignored.

See LP / QP Solver Routing for how classification works and when it falls back to the more general solver.

The shared backbone

Every interior-point and active-set solver above assembles a symmetric KKT system and factorizes it through pounce-linsol. That trait layer is backend-agnostic:

• FERAL (pounce-feral) — a pure-Rust sparse symmetric LDLᵀ factorization. The default; no external dependencies.
• HSL MA57 (pounce-hsl) — the well-known Harwell solver via libcoinhsl, enabled with the ma57 build feature for large or ill-conditioned systems.

Because the backend is pluggable, the same solver code runs on either without change.

Cross-cutting layers

These are not solvers you select, but stages and tools the solvers share:

• Presolve (pounce-presolve) — an optional front-end that tightens bounds (feasibility-based bound tightening), removes redundant rows, and repairs LICQ degeneracies before the solve.
• Restoration (pounce-restoration) — the feasibility-recovery phase the filter-IPM enters when a step cannot reduce both infeasibility and the objective; pounce-l1penalty offers an ℓ₁-exact penalty reformulation for degenerate / LICQ-violating problems.
• Sensitivity — pounce-sensitivity gives sIPOPT-style parametric steps and reduced Hessians for the NLP; QpSensitivity does the same for the convex QP. See Sensitivity Analysis.
• Cone library (pounce-convex) — nonnegative, second-order, exponential, power, and (for small dense problems) positive-semidefinite cones, so small SDPs solve as a convex class. The PSD cone cannot yet be mixed with the exponential/power cones in one problem (they use different drivers).
• Solve report — every path can emit the machine-readable pounce.solve-report/v1 JSON (status, iterations, residuals, timing). See JSON Solve Report.

Global vs. local — the honest summary

POUNCE settles a problem globally along two routes, and locally along one:

• Global by convexity — LP, convex QP, SOCP, and the exponential / power / PSD cone classes. Local is global, so a convex or conic reformulation buys the guarantee outright.
• Global by certificate (polynomials) — the SOS / Lasserre optimizer certifies the global minimum of a nonconvex polynomial from a single SDP; see Global Optimization.
• Local for general NLP — the filter-IPM and SQP paths converge to a KKT point, which for a nonconvex problem carries no global guarantee.

There is no third route: a nonconvex, non-polynomial problem gets a local answer, and no certificate.

Two practical levers for a “global” answer: modeling (cast as much as you can into the convex cone library) and, when that is not possible, the SOS / Lasserre optimizer for polynomials.

Running Solves

The pounce command-line driver solves built-in TNLPs and AMPL .nl files. Its console output mirrors upstream ipopt’s banner, per-iteration table, and final summary, so anyone used to reading ipopt logs can read pounce logs unchanged.

Basic usage

pounce problem.nl
pounce problem.nl print_level=8 max_iter=500 tol=1e-10
pounce problem.nl linear_solver=ma57            # with --features ma57
pounce problem.nl --options-file tuned.opt      # upstream-format options file
pounce problem.nl option_file_name=tuned.opt   # the Ipopt spelling; same thing
pounce problem.nl --no-options-file            # ignore ./pounce.opt, ./ipopt.opt

Trailing KEY=VALUE pairs follow the same syntax and semantics as the upstream Ipopt CLI; they override values loaded from the options file. With no options file named, ./pounce.opt or ./ipopt.opt is read if present, as ipopt reads ./ipopt.opt. See Solver Options.

Built-in problems

pounce --list-problems
pounce --problem quadratic
pounce --problem rosenbrock

• quadratic — min (x[0]-3)² + (x[1]-4)² (unconstrained, optimum (3, 4)).
• rosenbrock — min 100·(x[1]-x[0]²)² + (1-x[0])² (unconstrained, optimum (1, 1)).
• bounded-quadratic — quadratic with box bounds 0 ≤ x ≤ 2 (optimum at the upper corner (2, 2)).
• eq-quadratic — min x[0]² + x[1]² s.t. x[0] + x[1] = 1 (a single equality).
• circle — min x[0] s.t. x[0]² + x[1]² = 1 (a nonlinear equality).
• infeasible-eq — two contradictory equalities (x[0]+x[1]=1 and =2); exercises the infeasibility-detection path.

Run pounce --list-problems for the authoritative list.

Built-in problems have no .nl stub, so they only write a .sol file when --sol-output is given explicitly.

Degenerate / MPCC NLPs — the ℓ₁-exact penalty-barrier wrapper

For problems where the standard IPM thrashes in restoration because LICQ fails at the iterate (degenerate equalities, MPCC-like complementarity), enable the Thierry–Biegler ℓ₁-exact penalty-barrier wrapper:

pounce problem.nl l1_exact_penalty_barrier=yes

The wrapper turns every equality row c_i(x) = g_i into a slack-relaxed c_i(x) − p_i + n_i = g_i with (p_i, n_i) ≥ 0, augments the objective by ρ · Σ(p + n), and runs a Byrd–Nocedal–Waltz outer loop that escalates ρ until the slacks collapse (constraints satisfied) or saturate (locally infeasible problem detected). The user-visible (x*, λ*) are reported in the original variable space.

For everyday use, the simpler form is an auto-fallback:

pounce problem.nl l1_fallback_on_restoration_failure=yes

POUNCE first runs the standard solve. If it terminates in Restoration_Failed, Infeasible_Problem_Detected, Solved_To_Acceptable_Level, Maximum_Iterations_Exceeded, or Not_Enough_Degrees_Of_Freedom, the wrapper is invoked transparently and the result is promoted to Solve_Succeeded only if the retry succeeds. Otherwise the original status is preserved.

The tuning knobs are listed under Solver Options.

AMPL / Pyomo solver mode

AMPL drivers — and Pyomo’s ASL interface — invoke a solver as solver problem.nl -AMPL. Pass -AMPL to run pounce that way:

pounce problem.nl -AMPL

It changes nothing about the solve itself; it switches the process to the AMPL exit-code contract (see below), so the driver reads the termination from the .sol file rather than the exit status. The pyomo-pounce package builds on top of this mode.

Dual sign conventions

Two different quantities are in play, and they differ by a sign. Getting them confused is easy, so pounce names them differently:

They are related by marginal = −λ. The .sol writer performs this negation, so a .sol (and therefore Pyomo’s model.dual) carries shadow prices with the same sign Ipopt, glpk, CBC and CONOPT report. For min x₀²+x₁² s.t. x₀+x₁ = b the optimum is b²/2, so a .sol written at b = 2 reports +2.

The Python API’s mult_g keeps the Lagrange-multiplier convention, which matches cyipopt and satisfies the stationarity identity above directly. If you are comparing a mult_g against a model.dual for the same solve, expect them to differ by a sign — that is by design, not a bug.

Before v0.9.1 the .sol writer emitted the Lagrange multiplier without converting it, so every AMPL/Pyomo dual came back negated relative to every other solver (#271). Objectives and primal solutions were never affected.

Exit codes

• 0 — Solve_Succeeded (or Solved_To_Acceptable_Level).
• non-zero — any other ApplicationReturnStatus.

In AMPL solver mode (-AMPL) the exit code instead follows the AMPL contract: 0 for any solve that ran and produced a .sol file — limit-reached, infeasible, even a failed solve — since the termination is carried by the file’s solve_result_num. Genuine startup failures (unreadable .nl, bad option) still exit non-zero.

Diagnostics & introspection

pounce --about                                   # version, build info, features, backends
pounce problem.nl --dump kkt:5-10 --dump iterate # dump per-iteration diagnostics
pounce problem.nl --dump kkt --dump-dir /tmp/d   # override the dump root

• --about — print version, build info, enabled features, and linear-solver backends, then exit.
• --dump <cat>[:<spec>] — write the diagnostic category to per-iteration files (JSONL). Wired categories are kkt and iterate; an optional :<spec> selects iterations (e.g. kkt:5, kkt:2-10, iterate:all).
• --dump-dir <path> — override the dump root (default ./pounce-dump-<timestamp>).
• --dump-format <fmt> — dump format (default jsonl).

Help

pounce --help
pounce --version          # also -v, -V

Solver Options

POUNCE accepts options the same way upstream Ipopt does. Option names and semantics follow Ipopt’s, so an existing Ipopt options file or KEY=VALUE invocation works unchanged.

Setting options

On the command line — append KEY=VALUE pairs after the input:

pounce problem.nl tol=1e-10 max_iter=500 print_level=8

From an options file — upstream ipopt.opt format (one name value pair per line, # comments):

pounce problem.nl --options-file tuned.opt
pounce problem.nl option_file_name=tuned.opt   # the Ipopt spelling; same thing

From an options file nobody named. With neither of the above, POUNCE looks in the working directory for pounce.opt, then ipopt.opt, and reads the first one it finds — the way ipopt picks up ./ipopt.opt:

printf 'max_iter 5\n' > ipopt.opt
pounce problem.nl        # → Using option file "ipopt.opt".

The run says which file configured it (on the same sb gate as the banner). If both default names are present, pounce.opt wins and the other is named in a warning rather than left to look applied. --no-options-file skips the lookup entirely — the escape hatch for a directory holding an options file written for some other run.

Command-line KEY=VALUE pairs — and $pounce_options, AMPL’s <solver>_options channel — override values loaded from the options file, never the reverse.

An options file that was named but cannot be read is an error, not a shrug:

$ pounce problem.nl option_file_name=typo.opt
pounce: failed to load options file: options file "typo.opt" does not exist. …

Upstream opens a named file with a bare ifstream and reads nothing if that fails, so a typo there runs at stock defaults without a word. That silence is the one thing this path exists to remove (#518): a run configured entirely through an options file that quietly ran at defaults still reported success, which invalidates any benchmark set up that way.

One upstream quirk is inherited: option_file_name set inside an options file chains nowhere, because the file has already been chosen by the time it is read. POUNCE warns about that rather than ignoring it.

Commonly used options

For the full upstream option catalogue, see the Ipopt options reference; POUNCE reuses those names.

For scaling-specific options (nlp_scaling_method, target-gradient overrides, linear_system_scaling), see the Scaling reference page. For nonlinear bound tightening (presolve_fbbt, fbbt_tol, fbbt_max_iter, fbbt_max_constraints), see the FBBT reference page.

Options POUNCE does not implement

POUNCE’s option registry is a faithful port of Ipopt’s: every name Ipopt registers is registered here, so an ipopt.opt written for Ipopt parses unchanged. Registering an option is not the same as implementing it, though — and for a long time setting an unimplemented one did nothing at all, silently.

Options naming a feature POUNCE does not have now fail the solve, naming the option, the feature, and what to use instead:

$ pounce model.nl dependency_detector=mumps
pounce: `dependency_detector` configures linear-dependency detection on the
equality constraints, which pounce does not implement. It is registered so an
ipopt.opt written for Ipopt still parses, but setting it used to do nothing at
all — silently — so it is refused instead. Instead: pounce's presolve removes
structurally redundant rows; see `presolve`. Remove it to run.

The features in question: the Chen-Goldfarb (CG-penalty) / inexact-Newton line search, derivative approximation by finite differences, linear-dependency detection, the per-iteration NaN/Inf derivative check, multiplier recalculation by least squares, a selectable constraint-violation norm, magic steps, bound replacement, the L-BFGS augmented-system variants, skipping the finalize callback, the dynamic HSL loader, and suppress_all_output / debug_print_level.

option_file_name was on that list until #518 implemented it; refusing an option is the cheap half of “implement it or fail loudly”, and an entry leaves this table by getting the other half.

Two deliberate exceptions:

• Setting an option to its registered default is allowed. A generated ipopt.opt spells out defaults, and dependency_detector=none asks for nothing. Only a value that differs from the default is a request POUNCE cannot honour.
• Caching hints warn instead of failing. grad_f_constant, hessian_constant, jac_c_constant and jac_d_constant tell the solver a quantity does not change between iterations. POUNCE re-evaluates regardless, so ignoring them costs evaluations and never correctness — failing the solve would be a worse trade.

Options whose feature runs and whose value simply is not read yet are not in this category; they still solve, with the default in effect. Wiring those is tracked on #191 and #483.

Derivative checker

Wrong analytic derivatives are the most common cause of an NLP that stalls, cycles, or converges to something that is not a solution — and they are invisible from the iteration log. derivative_test compares what your TNLP returns against finite differences at the (bound-projected) starting point, before the solve:

pounce problem.nl derivative_test=first-order

first-order checks eval_grad_f and eval_jac_g; second-order adds eval_h; only-second-order checks the Hessian alone. The Hessian is checked one multiplier block at a time — obj_factor = 1, λ = 0 against differences of eval_grad_f, then obj_factor = 0, λ = eⱼ against differences of row j of eval_jac_g.

Entries that look wrong are marked *:

Derivative checker: first derivatives at the starting point (perturbation 1.0e-8, tolerance 1.0e-4).
* grad_f[    1]       =    3.5000000000000000e0    ~    3.0000000119209290e0  [  1.667e-1]
1 suspicious derivative(s) and 0 missing sparsity entrie(s) out of 6 checked (8 evaluations).

Two checks beyond upstream Ipopt’s, because both catch a class of bug no value-by-value comparison can:

• A Jacobian or Hessian entry whose finite difference is nonzero but which the sparsity structure omits (! in the report). A missing structural entry is not a wrong number — it is a derivative the solver can never see.
• The perturbation is taken downward when stepping up would leave a variable’s box, so a model using sqrt, log, or 1/x is not evaluated outside its own domain by the checker.

The test is advisory: it reports and the solve continues. It is written to stderr, so it survives print_level=0 and never mixes into --json-output’s stdout. It is slow — the second-order test costs roughly (m+1)·n evaluations — so leave it off for production runs.

check_derivatives_for_naninf is a separate upstream option, for a per-iteration NaN/Inf guard, and is not implemented.

Choosing a linear solver

POUNCE implements two KKT backends:

• feral — pure-Rust sparse symmetric indefinite solver. The effective default; no Fortran toolchain, no HSL licence.
• ma57 — HSL MA57, available only in a cargo build --features ma57 build.

The option’s registered value list is a faithful port of upstream Ipopt’s (ma27, ma77, ma86, ma97, mumps, pardiso, pardisomkl, spral, wsmp, custom), so an ipopt.opt written for Ipopt parses here unchanged. Selecting one of those fails the solve with a message naming it. They used to fall through to FERAL silently, which meant linear_solver=ma97 “worked” and a benchmark comparing backends compared FERAL with itself.

The registered default is feral, which diverges from upstream’s ma57 on purpose: a default has to name a solver the binary actually contains. Under the upstream default a pure-Rust build advertised MA57 to every print_user_options dump while running FERAL, and an HSL-enabled build used MA57 without being asked. If you build --features ma57 and want it, select it explicitly — that is the one behavioural change here.

Not a failure: explicit ma57 on a build without the feature falls back to FERAL and says so in the banner (FERAL (ma57 requested but not compiled)). That substitution is reported rather than hidden, and failing a portable ipopt.opt over a build flag would cost more than it buys.

The per-backend tuning options (ma97_scaling, mumps_pivtolmax, pardiso_*, wsmp_*, spral_*, …) remain registered for the same ipopt.opt-compatibility reason. They are unreachable now that their backend cannot be selected.

Bound relaxation and honor_original_bounds

Before the solve, POUNCE widens every variable and constraint bound by bound_relax_factor (default 1e-8, capped by constr_viol_tol), exactly as upstream Ipopt does — it keeps the interior-point iterates strictly feasible without the user’s bounds becoming numerically degenerate. The consequence is that a solution pinned to a bound is reported just past it:

min (x − 3)²  s.t.  0 ≤ x ≤ 1     →     x = 1.00000000937

honor_original_bounds=yes projects the reported point back into the bounds you declared, so that solve returns exactly x = 1. Reach for it whenever the value flows somewhere that cares about the domain — a sqrt(1 − x), a domain assertion, or a Pyomo Var the value is loaded back into.

The default is no, matching upstream. As upstream also documents, the constraint-violation and complementarity figures in the end-of-run summary are for the non-projected point; only the reported x (and the objective and constraint values evaluated at it) move.

Large constraint values and primal_noise_floor_kappa

On a model whose constraint values run to ~1e7 and beyond, a converged solve could exit Search_Direction_Becomes_Too_Small while holding the correct optimum. The cause is arithmetic, not the model.

The KKT error the convergence test compares against tol is

max( ‖∇L‖∞ / s_d ,  max(‖c‖∞, ‖d − s‖∞) ,  ‖compl‖∞ / s_c )

The dual and complementarity terms are normalised; the primal one — like upstream Ipopt’s — is a bare absolute residual. But c_i = g_i(x) − b_i and d_i − s_i are each a difference of quantities the row’s own size, so they are quantised in units of eps · |b_i|. At |b| ~ 1e8 the smallest nonzero value the primal term can take is one ulp, 1.5e-8 — already larger than the default tol = 1e-8. Asking for nlp_err <= tol there is asking the residual to land on a bitwise-exact 0 rather than on one ulp, which is arithmetic luck rather than a property of the iterate.

POUNCE therefore judges the primal term in the strict test against each row’s own floating-point resolution: a row’s residual counts only where it exceeds max(placement floor, kappa · eps · |row magnitude|), with kappa = primal_noise_floor_kappa (default 64). Verdicts are flat across kappa from 8 to 1024 on the measured set.

Three things bound what this can do:

• Only the strict test reads it. constr_viol is still checked against constr_viol_tol (default 1e-4) on the full, unfloored residual, so nothing the floor forgives can exceed the feasibility tolerance you set — however large your data grows.
• The acceptable-level band keeps the raw error. It sits two decades above tol, clear of any realistic quantum.
• It cannot rescue an infeasible model. On a model with no feasible point the filter and restoration phase reach a verdict on their own criteria; the floor only ever participates at a point the rest of the algorithm already believes is converged.

Set primal_noise_floor_kappa = 0 to switch the floor off and restore upstream Ipopt’s bare-absolute primal term exactly.

When the floor changes the reported picture, the end-of-run summary says so — a large-|b| solve prints the tested value under the raw one:

Overall NLP error.......:   2.3841857910156250e-07    2.3841857910156250e-07
  ...above the per-row floating-point noise floor:   0.0000000000000000e+00

Solves where the two agree — every model whose data is O(1) — print the usual block unchanged.

One case this does not paper over: tightening constr_viol_tol below a row’s own ulp (say 1e-8 on data at 1e8) still will not certify. That is the tolerance gate doing what you asked — the residual you requested is not representable at that scale.

Solved_To_Acceptable_Level and acceptable_progress_kappa

Solved_To_Acceptable_Level is the fallback verdict for a solve that cannot reach tol: after acceptable_iter (default 15) consecutive iterates with an NLP error under acceptable_tol (default 1e-6), the solver stops and hands back the point it has. That criterion is a count of iterates inside a band, and on its own it asks only is the error small — never has anything stopped moving.

Those come apart. An interior-point iterate can be near-stationary for the current barrier subproblem — a much weaker statement than near-KKT for the NLP — for fifteen iterations running while the solve is still descending. Two measured cases: the kissing model stopped with objective 1.00000108 where continuing reaches 0.84544259 and a strict certificate, 18% lower; NARX_CFy stopped with both residuals near 1e-7 where sixty more iterations collapse them by five orders.

POUNCE therefore also requires the streak to have flattened. Across the acceptable_iter iterates that made it up:

• the spread (max − min) of the NLP error must be within acceptable_progress_kappa · acceptable_tol; and
• the spread of the objective within the same fraction of acceptable_tol · max(1, |f|).

acceptable_progress_kappa defaults to 0.1, so at default tolerances both quantities must have stayed inside a tenth of the acceptable band over the whole streak.

It is a spread, not a trend, and either signal alone is enough to keep solving. Both choices are deliberate: kissing’s error was an order of magnitude worse at the iterate it stopped on than at one it had already reached inside the same streak — it was wandering across the band, not converging inside it — while its objective was flat to all eight printed figures over the same iterates.

Three things bound what this can do:

• It cannot lose a verdict. The refused termination is recorded, and a run that fails to do better ends at exactly that iterate under exactly that status. A misfire costs iterations, never the answer — you will not see Maximum_Iterations_Exceeded where the count alone would have said Solved_To_Acceptable_Level.
• It never looks at a solve that converges. A solve that reaches tol never completes an acceptable-level streak, so nothing here runs.
• A genuine stall flattens. When the iterate, the objective and the error are all pinned — the case the acceptable-level exit exists for — the window is flat and termination happens as before.

Set acceptable_progress_kappa = 0 to switch the progress test off and restore upstream Ipopt’s bare consecutive-count criterion. Widening acceptable_tol widens the flat bar with it, so asking for a looser band still gets you the early exit.

Big models that start feasible and the theta_max ceiling

The filter has a hard ceiling. Any trial iterate whose constraint violation θ exceeds

theta_max = theta_max_fact · max(1, θ₀)

is rejected outright, before any of the filter’s usual tests run. It is a global-convergence safeguard: it keeps the line search from wandering arbitrarily far from feasibility.

The trouble is the 1. POUNCE’s θ is a 1-norm over constraint rows — ‖c‖₁ + ‖d − s‖₁, a sum of m residuals — so a ceiling of T really says “a mean per-row violation of T/m”, and that allowance shrinks as the model grows. And on a problem started at a feasible point, θ₀ = 0, the max collapses and the ceiling is the bare constant theta_max_fact however large the model is.

robot_a is the measured case: 52 013 constraint rows, a feasible start, so theta_max locked at 1e4 — a mean per-row allowance of 0.19 — while the route to the optimum passes through θ ≈ 9.4e7. Every step toward the solution was refused at the gate, and the solve ground to its iteration limit at objective 8.173304 instead of the true 1.0431952.

POUNCE’s answer is theta_max_adaptive_trigger (default 3), described below. theta_max_row_scale_kappa (default 0, off) is an earlier, static attempt at the same problem, kept because it is occasionally the more direct lever; it floors the reference at the row count instead:

theta_max = theta_max_fact · max(θ₀, theta_max_row_scale_kappa · rows, 1)

so the ceiling means a mean per-row violation of theta_max_fact regardless of m. Measured under defaults otherwise, against Ipopt 3.14 on the same machine:

Ipopt has the same defect and no correction for it; all three solve under theta_max_fact = 1e8 set by hand, which is the blunt version of the same move.

The adaptive rule: theta_max_adaptive_trigger

On by default. Rather than guessing from problem size whether a model needs headroom, POUNCE measures whether the ceiling is what is refusing the line search, and raises it only then.

A trial refused because θ_trial > theta_max takes a distinct early exit, before the filter and Armijo tests run at all. So the acceptor can count those refusals and compare them against the number of trials attempted. When every trial of a line search was refused at the gate, for theta_max_adaptive_trigger consecutive line searches, the ceiling is demonstrably the binding constraint — not the filter — and it is multiplied by theta_max_adaptive_factor (default 100), at most theta_max_adaptive_max_raises times per solve (default 4).

Three properties follow, and they are what the static floor could not offer:

• A converging model cannot trip it. Converging means trials are getting past the gate; a model accepting steps never accumulates the streak. brainpc1/3/5/7 — the family the static floor damaged at every kappa — are untouched by construction, not by a lucky constant.
• A blocked model trips it immediately. robot_a is refused at the gate from its first line search onward.
• The ceiling stays finite. Wächter–Biegler’s global-convergence argument (Thm. 2) needs theta_max finite, not fixed. A bounded number of bounded raises keeps it finite, so a solve cannot ratchet the safeguard away one line search at a time.

Requiring a streak rather than a single line search is deliberate: one Newton direction that overshoots into a huge θ can legitimately have all its trials refused, and backtracking is the right response to that. Only a model that cannot get past the gate repeatedly is one whose route needs the headroom.

Measured, defaults otherwise:

Across the whole Vanderbei corpus (733 problems) the rule changes four outcomes: britgas goes from its iteration limit to Optimal in 16 iterations, catenary solves to the same objective in 50 iterations instead of 56, and coshfun and brainpc0 fail either way — coshfun now reporting diverging iterates, which is what Ipopt 3.14 also does on it. Net Optimal count 702 → 703.

Note brainpc0 does trip the rule while brainpc1/3/5/7 do not, despite identical row counts. That is precisely the distinction a size-based floor cannot draw.

The restoration sub-IPM always runs with the rule disabled. Upstream already corrects the resto phase’s instance of this degeneracy by hard-coding resto.theta_max_fact = 1e8 (IpRestoMinC_1Nrm.cpp:91), so a rule that ratchets further would be compounding a correction already made.

Set theta_max_adaptive_trigger = 0 to restore upstream Ipopt’s fixed ceiling exactly.

When to reach for the static floor instead

Symptoms, all three together:

• a large number of constraint rows (thousands upward);
• a feasible or near-feasible starting point — the iteration log’s first inf_pr is 0 or very small;
• the solve stalls with inf_pr flat and alpha_pr tiny, and raising max_iter does not help.

The quick confirmation is to set theta_max_fact = 1e8 by hand. If that unsticks the model, theta_max_row_scale_kappa = 1 is the principled version of it — it scales the ceiling to the model rather than to a constant you picked.

Why the static floor is off by default

Because raising the ceiling unconditionally is not free. It relaxes a global-convergence safeguard, and a model that was not being blocked by it can wander instead. On the Vanderbei corpus, brainpc1/3/5/7 (m = 6900, θ₀ = 1e-2) all regress — brainpc1 from Optimal in 64 iterations to divergent, objective 3.7e3 against the correct 4.4e-04.

A scan over kappa showed why this cannot be tuned away. The damage is a step function, not a gradient:

brainpc3 and brainpc7 land on the identical worse answer at every nonzero kappa, even 0.01 — where the ceiling moves only from 1e4 to 6.9e5. The instant it rises at all, they break. robot_a meanwhile improves monotonically all the way to kappa = 1. There is no separating value.

That is a verdict on the design, not on the tuning: the real question is whether a model’s route to the optimum needs the extra headroom, and the row count does not answer it. A static floor cannot know — which is why the adaptive rule above, which asks the question directly, is the default and this one is not.

What still bounds the option when you do turn it on:

• It only ever raises the reference. A model whose own θ₀ already exceeds kappa · rows gets exactly upstream’s ceiling.
• It reduces to upstream on a single-row problem, where the floor is max(kappa · 1, 1) = 1 — upstream’s constant.
• theta_max is still finite, and still fixed for the whole solve after its first line search. This rescales the safeguard; it does not remove it.
• The restoration sub-IPM is untouched, at any kappa. Upstream already fixes its own instance of this by hard-coding resto.theta_max_fact = 1e8 (IpRestoMinC_1Nrm.cpp:91) — the resto NLP is also initialised feasible, so it hit the same degeneracy. Stacking the row floor on top would push that inner ceiling to 1e8 · m, i.e. remove it, so the sub-IPM always runs with kappa = 0.

Large gradients and dual_inf_scale_kappa

The dual side of the same story. dual_inf_tol (default 1.0) is a bare absolute bound on ‖∇L‖∞ — but the aggregate above normalises that quantity, dividing it by s_d, which grows with the mean magnitude of the multipliers. On a model whose gradients live at 1e10 the two gates are judging one number by standards ten orders apart.

Vanderbei’s orthrds2 is the reported case: s_d ≈ 1.6e10 with ‖∇L‖∞ = 89.7, so the aggregate’s dual term is 5.6e-09 — comfortably inside the default tol = 1e-8, i.e. stationary to nine digits relative to the size of the gradients involved — while the component gate refused it against 1.0. The solve exited Solved_To_Acceptable_Level holding the answer, and dual_inf_tol=1e3 alone turned it into Optimal Solution Found at the same objective.

The simplest statement of the defect: multiply an objective by a positive constant. Same feasible set, same solution, same active set, same Newton step — and every multiplier, s_d and ‖∇L‖∞ scale with it, so a large enough constant costs the certificate.

The strict test therefore judges the unscaled dual infeasibility against

max( dual_inf_tol ,  kappa · tol · dual_scale )

with kappa = dual_inf_scale_kappa (default 1) and dual_scale the magnitude of the largest single term ∇L is assembled from (∇f, Jᵀy, the bound multipliers). Since ∇L is the sum of those terms, ‖∇L‖∞ / dual_scale is the fraction of them that failed to cancel — a scale-invariant statement of stationarity, and the thing the absolute bound was standing in for.

What bounds it:

• It cannot forgive non-stationarity. A point where nothing cancelled has ‖∇L‖∞ ≈ dual_scale, a ratio of 1 against a bar of 1e-8. min −exp(x) s.t. x >= 0 running away to inf_du = 8.8e+47 is refused by eight orders, because its ∇f runs away by exactly the same factor.
• The aggregate still has to pass. nlp_err <= tol is tested on the same iterate; this only removes the second, inconsistent standard.
• It is inert on ordinary models. At the defaults the floor does not rise above dual_inf_tol until dual_scale exceeds dual_inf_tol / tol = 1e8, so every model with O(1) gradients keeps upstream’s comparison bit for bit.
• Only the strict gate reads it. acceptable_dual_inf_tol (1e10) is untouched.

Set dual_inf_scale_kappa = 0 to switch the floor off and restore upstream Ipopt’s bare-absolute bound. That is also the setting to reach for if you tighten dual_inf_tol and want that absolute standard honoured unconditionally — the floor is a floor, so it can override a tightened dual_inf_tol on a large-gradient model.

Objective sense and obj_scaling_factor

obj_scaling_factor multiplies the objective the IPM minimizes, so a negative value maximizes — upstream’s documented spelling for a maximization problem stated as a minimization. Because it changes what is being optimized rather than just its conditioning, it is honored only by the general NLP interior-point path: a model that would otherwise route to the specialized convex solvers (LP / convex QP / SOCP, see LP/QP Routing) is re-routed under solver_selection=auto, and an explicit convex solver_selection is refused rather than silently answering with the minimizer.

A positive factor is a pure conditioning knob; the convex path reports natural units either way, so it keeps the fast path.

Barrier-parameter (μ) strategy

The barrier parameter μ controls the inner subproblem’s relaxation of complementarity. The two strategies are monotone (default — geometric schedule) and adaptive (quality-function oracle picks each μ from the current iterate’s complementarity). See μ-strategy for when to switch.

Quality-function oracle (adaptive-μ details)

These are only consumed when mu_strategy=adaptive and mu_oracle=quality-function. Defaults mirror upstream IpQualityFunctionMuOracle::RegisterOptions.

Adaptive-μ globalization

Tuning the safeguards that fall back to monotone-μ mode when the adaptive oracle stops making progress. Defaults mirror upstream IpAdaptiveMuUpdate::RegisterOptions.

ℓ₁ penalty-barrier wrapper options

These tune the degenerate-NLP wrapper described in Running Solves. All are default-tuned and rarely need overriding:

NLP Presolve

POUNCE’s TNLP-wrapper presolve pipeline runs before the IPM starts. It tightens variable bounds, drops redundant rows, and (optionally) eliminates square auxiliary-equality sub-systems structurally. All are off by default — set the master switch first:

presolve=yes applies equally to CLI solves and to every IpoptApplication::optimize_tnlp library solve; callers no longer need to wrap a callback TNLP manually. The wrapper postsolves before finalize_solution, so callback payloads remain in the submitted TNLP’s original variable and constraint space. Bare callback TNLPs do not expose an expression provider, so presolve_fbbt=yes remains a no-op for this library entry point.

Linear-equality variable elimination (Phase 6)

presolve_linear_eq_reduction=yes is the only pass that removes columns. It reads the model’s linear equality rows and eliminates the variables they determine, iterating to a fixed point so chains propagate:

• a variable whose declared bounds are equal becomes a constant;
• a singleton row a·x = b pins its variable at b/a;
• a two-variable row a₁·x + a₂·y = b substitutes one variable for the other, x := α·y + β. There is no anchoring requirement: a row linking two otherwise-free interior variables — an arc equality, a Reference alias, a unit-conversion link — aggregates away, which is the case the auxiliary-equality pass cannot reach because it only solves determined square blocks.

Rows that collapse to 0 = 0 under the accumulated substitutions are dropped as structurally redundant.

A row written with a constant on the left — x0 − 2·x1 + 3 = 3 — is eligible on the same terms as x0 − 2·x1 = 0. The .nl reader folds a constant row body into the row’s bounds when the file is read, so the pass sees an ordinary linear equality.

Every eliminated variable’s bounds are transferred onto its survivor, so the reduced box is never looser than the original. finalize_solution lifts the primal back to the original variable order and recovers a multiplier for each consumed row, so .sol / JSON solution blocks keep the original model’s shape and can still be read positionally by AMPL or Pyomo.

Three things to know before turning it on:

• Dual attribution. A transferred bound’s multiplier comes back on the variable that declared the bound, not on the survivor that inherited it. The plan records where each reduced bound came from, and postsolve rescales the multiplier by the substitution’s coefficient α — and moves it to the other side of the box when α < 0, since a negative coefficient turns a lower bound into an upper one. On a model with a single active transferred bound the reported duals match a no-presolve solve exactly.

  Where the survivor’s own bound and a transferred bound are active at the same point, the split between the two multipliers is genuinely non-unique — the reduced problem has one where the full problem has two — and the pass leaves the whole multiplier on the survivor. That is a valid KKT point, but it is not the split a no-presolve solve happens to report. The same holds for a variable pinned by a singleton row a·x = b whose value lands on one of its own bounds: the row multiplier absorbs it.

  A practical consequence for .sol readers, unchanged by any of the above: the writer omits exact zeros from suffix blocks (it always has), so a variable whose bound multiplier is zero gets no ipopt_zL_out / ipopt_zU_out entry at all rather than an entry of zero. Code that indexes those suffixes must treat a missing index as zero — as it already must for any variable whose bound multiplier lands exactly on zero. Row multipliers are unaffected: the dual block is dense and comes back at the original row count.

• Bounds the reduced problem never saw. Re-attribution can only move a multiplier the solver reported, and sometimes there is none. The transfers can leave a survivor’s reduced box as a single point, and a variable with equal bounds is a fixed variable, which the solver drops — so it comes back with no bound multiplier at all even though the cluster it stands for is sitting on a bound that needs one. Postsolve fills that in: whatever stationarity residual the recovered row multipliers cannot close is a bound multiplier that was never reported, and it goes on the declared bound the point is actually resting on — the survivor’s own where that is the active one, otherwise the column the bound was borrowed from, through the same α rescale as above. A residual with no active declared bound to carry it is left alone rather than parked somewhere that would break complementarity.

  One column is deliberately outside this: one the model declares fixed (x_l == x_u). The solver drops those as parameters whether or not the reduction runs, and reports no multiplier for them either way, so nothing here changes what they report.

• Failing closed. If the equality system is contradictory, the pass stands down entirely and hands the model to the solver untouched, rather than being the first and only voice to call a model infeasible. The same goes for a model whose every column is determined: a zero-variable problem is not a shape worth handing the IPM.

It is off by default because it changes the variable count, which the sensitivity and reduced-Hessian paths index against the original .nl. (The CLI already disables presolve entirely when those are requested.)

LP and convex QP take a different route to the same reduction. Those models never reach Phase 6 — the CLI dispatches them to pounce-convex before any presolve wrapper is built — but they are not left unreduced. pounce-convex has its own presolve, on by default, and it now performs the two-variable aggregation as part of that catalog, sharing this planner rather than restating it. So the reduction is the same; only the switch differs (qp_presolve=no / presolve=no turns it off there, and presolve_linear_eq_reduction does not apply). See LP/QP routing.

The two agree on dual attribution as well: a transferred bound’s multiplier is reported on the column that declared the bound, not on the survivor that inherited it. They get there differently — Phase 6 records during planning where each reduced bound came from, while the convex path reads the leftover reduced cost at postsolve and hands it to whichever column is sitting on its own bound, because it also has inequality rows and its own bound-tightening layer to account for. Where the survivor’s own bound is active as well, both leave the multiplier on the survivor; the split is genuinely non-unique there and either answer is a valid KKT point.

Feasibility-based bound tightening (Phase 1b)

Interval-arithmetic propagation through nonlinear constraint expression DAGs (see FBBT). Available today for .nl-loaded problems via NlTnlp; other TNLP sources opt out silently.

Auxiliary-equality preprocessing (Phase 0)

A separate set of options controls the structural elimination pass documented in Auxiliary-Equality Preprocessing:

FERAL backend tuning

linear_solver=feral (the default — see Commonly used options) is configurable through seven feral_* options. Defaults are tuned for the IPM workload and rarely need changing; reach for these when profiling a specific problem. Each also falls back to a matching POUNCE_FERAL_* environment variable when left unset on the OptionsList (see Environment overrides).

feral_ordering variants

All six concrete and adaptive options live under the same string option. feral_ordering also falls back to the POUNCE_FERAL_ORDERING environment variable when not set on the OptionsList.

When in doubt: leave feral_ordering at the default. When a hard problem looks linear-solver-bound, try feral_ordering auto_race before per-variant manual sweeping — it’s the safe choice when the per-problem winner is uncertain.

Caller-supplied ordering (External)

Beyond the string variants above, a structure-aware caller can inject a precomputed permutation the generic AMD/METIS pass cannot see — a block-triangular / Schur ordering (Parker, Garcia & Bent, arXiv:2602.17968) or a tearing ordering from equation-oriented decomposition. Because a permutation is a vector it cannot travel through the string feral_ordering option; supply it programmatically instead:

• Python: Problem.set_ordering(perm) (and get_ordering() / clear_ordering()) — see the Python guide.
• Rust: IpoptApplication::set_external_ordering(perm).

perm is a 0-based, new-to-old permutation (perm[k] is the original index that becomes index k) whose length must equal the augmented KKT system dimension (variables + slacks + constraint duals), not the problem’s n. FERAL validates it as a bijection and fails the factorization with an error on a wrong length or duplicate — a valid but poor ordering only costs fill/time, never correctness. This maps to FERAL’s OrderingMethod::External (feral#107) and honors only the default FERAL backend.

Environment overrides (FERAL and debug gates)

A handful of knobs are reachable through environment variables. The feral_* numerics knobs read their POUNCE_FERAL_* variable only as a fallback when the matching option is left unset on the OptionsList — set the option (per solve, recordable, discoverable via the debugger’s opt command) in preference to the env var (process-wide, invisible to the solve report). The debug gates below have no option equivalent; they exist purely to switch on extra diagnostic output.

FERAL numerics fallbacks

Each maps one-to-one to a registered option in FERAL backend tuning. Prefer the option; the env var is the fallback for callers with no OptionsList (some tests, legacy embeddings).

FERAL_PARALLEL (legacy, no POUNCE_ prefix) forces feral’s internal factor serial or parallel process-wide — 0/off/false/no to force serial, 1/on/true/yes to force parallel, and unset to leave feral’s own platform-derived default alone. The force-on direction is the only override available to CLI, Python and NL callers on a host where that autodetection is wrong (feral falls back to sequential when the rayon pool fails to build); the first-class per-backend lever, FeralConfig.parallel, is the Rust solver API, not an option.

Debug and diagnostic gates

These switch on extra diagnostic emission for a specific subsystem. Most emit at debug level under a pounce::* tracing target, so setting the gate alone is not enough — pair it with a matching RUST_LOG (e.g. RUST_LOG=pounce::mu=debug) or the output stays filtered. Presence-only unless a value is noted; they are diagnostic aids, not part of the stable interface, and may change between releases.

Two already-documented gates round out the set: POUNCE_DBG_LLM and POUNCE_DBG_VIEWER (see the debugger guide).

Logging and colored output

POUNCE emits structured logs and a colored iteration table through the tracing ecosystem. Behavior is governed by environment variables (not solver options), so they apply to the pounce CLI, the C/Python frontends, and anything embedding the library.

Filtering by subsystem. Solver internals log under namespaced targets — pounce::algorithm, pounce::linsol, pounce::mu, pounce::sqp, pounce::linesearch, pounce::restoration, pounce::presolve, pounce::py. For example, to trace only the restoration phase:

RUST_LOG=pounce::restoration=debug pounce problem.nl

Program output vs. logs. The iteration table, the final summary, and --dump diagnostics are program output on stdout; diagnostic and progress messages are logs on stderr. Redirecting one does not affect the other:

pounce problem.nl > result.txt 2> solve.log

Color. The iteration table is colored with a tiger/rust theme: restoration lines take a background that varies by restoration kind (soft-stay → tan, soft-exit → amber, hard → deep rust), and the row text shades from black toward red as the primal step length alpha shrinks (stalling). Color is emitted only when stdout is a terminal; redirected output and NO_COLOR get plain text with identical column alignment.

Machine-readable iterations. POUNCE_LOG_FORMAT=json turns the per-iteration records into JSON on stderr:

POUNCE_LOG_FORMAT=json pounce problem.nl 2> iters.jsonl

LP / QP Solver Routing

POUNCE can route linear programs (LP), convex quadratic programs (QP), and convex quadratically-constrained QPs (QCQP) to a specialized interior-point solver (pounce-convex) instead of the general nonlinear (NLP) filter-IPM. The specialized path uses Mehrotra predictor-corrector and reaches the solution in materially fewer iterations on these problem classes — typically 30–50% fewer than the general NLP path on bound- or inequality-constrained convex QPs.

Routing is automatic and transparent: you do not change how you call POUNCE. The same pounce problem.nl, the same SolverFactory('pounce') in Pyomo, and the same AMPL solve all work unchanged — POUNCE inspects the problem and picks the solver.

How routing works

When POUNCE loads a problem it classifies it into one of:

The classifier is conservative: a problem is sent to the convex solver only when POUNCE can prove it is convex — an LP or convex QP (degree-≤2 objective with a positive-semidefinite Hessian, linear constraints), or a convex QCQP (additionally allowing convex-quadratic inequality constraints, each with a positive-semidefinite Hessian and a one-sided ≤ bound, which are reformulated to second-order cones). Anything it cannot prove convex — transcendental terms, an indefinite objective Hessian, a quadratic equality, or a quadratic inequality whose feasible set is nonconvex — falls back to the general NLP solver, which always produces a correct (locally optimal) answer. You never get a wrong “optimum” from a misclassification.

Note on QP detection. The AMPL .nl format has no dedicated quadratic section: a QP’s quadratic terms are written into the nonlinear expression tree. POUNCE walks that tree to recover the Hessian and test convexity, the same way QP-capable AMPL solvers do.

Note on row constants. A .nl writer may leave a constant on the left of a constraint — x0 + x1 + 3 <= 6 rather than x0 + x1 <= 3 — and it too lands in the nonlinear expression tree. The reader folds such a constant into the row’s bounds when the file is read, so a model that is otherwise an LP still classifies as one. The shift is exact: body and bound move together, so the solution and every multiplier are the same as for the hand-folded model.

Choosing the solver explicitly

The solver_selection option overrides the automatic choice. It is a normal POUNCE option, so it works on the command line, in an options file, or through Pyomo’s solver.options.

# Let POUNCE decide (default):
pounce model.nl

# Force the NLP path even on a convex QP (e.g. to compare):
pounce model.nl solver_selection=nlp

# Insist the problem is a convex QP — fail loudly if it is not:
pounce model.nl solver_selection=qp-ipm

# Solve that same QP with the active-set engine instead of the IPM:
pounce model.nl solver_selection=qp-active-set

A forced value that does not match the detected class is rejected with a clear message rather than silently ignored:

pounce: problem class NLP does not match forced solver qp-ipm
        (expected an LP or convex QP)

qp-active-set hands the QP directly to pounce-qp’s ParametricActiveSetSolver, through the same convex driver the IPM uses — so it inherits presolve, postsolve, dual recovery, .sol writing, timing and the convex status vocabulary. It is not the same route as algorithm=active-set-sqp, which wraps the QP in the full SQP outer loop; that option still exists and is the right one for a genuine NLP.

Choose it deliberately. For a cold, one-shot convex QP the interior-point path (qp-ipm, and what auto selects) is materially more robust: on the 138-problem Maros-Mészáros set the IPM solves 137 while the active-set engine solves substantially fewer, mostly by exhausting its iteration budget on large degenerate instances. That is the expected character of a cold active-set method rather than a defect — its iteration count is combinatorial in the size of the active set, where an interior-point count is nearly independent of problem size. The active-set engine earns its keep on warm-started sequences — MPC steps, branch-and-bound nodes, continuation — where consecutive QPs differ little and the working set carries over; see solve_parametric.

What it will not do is lie: it reports Maximum_Iterations_Exceeded rather than a wrong answer, and every claimed optimum is re-verified against the original problem’s KKT conditions before being reported.

From Pyomo

solver = SolverFactory('pounce')
solver.options['solver_selection'] = 'qp-ipm'   # or 'auto', 'nlp', ...
solver.solve(model)

What you get back

Before solving, POUNCE prints a one-line routing banner naming the detected class, the solver it selected, and the effective solver_selection — so it is always clear which of POUNCE’s solvers ran and why:

Problem class: LP. Selected solver: convex QP interior-point (pounce-convex) [solver_selection=auto].

(The banner is suppressed alongside the startup banner — sb yes or JSON-debug protocol mode — to keep stdout clean for machine consumers.)

The convex IPM then reports the same way as the NLP path: an optimal-status line, the objective value (in your original sense — a maximize objective and any constant term are reported correctly), and a .sol file with the primal solution when one is requested.

POUNCE (LP IPM, pounce-convex): Optimal Solution Found.
        obj=2.00000000  iters=2

Driver. The convex path uses the homogeneous self-dual embedding (HSDE) interior-point driver — the same self-dual formulation Clarabel/ECOS use. It is self-starting, returns verified infeasibility/unboundedness certificates, and conditions the KKT system internally through its per-cone scaling, so it solves even badly-scaled LPs (e.g. NETLIB nl, ‖c‖ ~ 1e6) without external pre-scaling.

Presolve

Before the convex interior-point solve, POUNCE runs a presolve pass that shrinks the problem and can detect trivial infeasibility or unboundedness without solving. It removes empty, duplicate, and activity-redundant rows; fixes and substitutes structural columns (singleton-row fixings, free columns, free column singletons); folds away two-variable equality rows (below); and recovers both the primal and dual of the eliminated pieces so the reported solution is for your original problem. When it reduces the model, it logs a one-line summary:

Presolve: 40 → 24 vars, 12 → 4 rows (fixed 3, free-fixed 2, substituted 3, aggregated 8, ...)

Two-variable equality rows (aggregation)

A row a₁·x + a₂·y = b linking two variables says one of them is the other, up to a scale and a shift — an arc equality between two units, a Reference alias, a unit conversion. Neither variable is determined by it, so nothing in the older catalog could act on it, and on a flowsheet these rows are most of the model. POUNCE now substitutes one variable for the other and drops the row, iterating to a fixed point so chains of aliases collapse to a single column. Any bound on the eliminated variable is carried across onto the one that survives, so the reduced problem describes exactly the same feasible set.

Two things this deliberately does not do:

• It never calls your model infeasible. A contradictory alias system — x = y and x = y + 1 — makes the pass stand down and hand the model over untouched, for the rest of presolve or the solver itself to judge.
• It does not run on the conic path (SOCP, exponential/power cones, SDP, SOS). Those rows are structurally coupled in fixed-size blocks that a substitution would rewrite.

The aggregation shares its planner with the NLP path’s Phase 6, so the two agree on what can be eliminated (see NLP Presolve).

Infeasibility verdicts are re-derived before they are reported

A presolve infeasibility comes back in milliseconds with no iteration behind it, so when it is wrong it is the most expensive answer the solver can give. Two reductions — forcing constraints and dominated columns — fix a variable at a value they choose from a tolerance judgment, and a fixing that is wrong is substituted into every row that variable appears in until some row reads as contradictory: a false infeasibility, reported against a row nowhere near the reduction that caused it.

So presolve does not report an infeasibility on the strength of the pass that found it. It re-derives the verdict from your original model with those two reductions switched off, and reports Infeasible_Problem_Detected only if that pass reaches the same conclusion on its own. If it does not, the model is solved normally and presolve says so:

Presolve: discarded an unconfirmed infeasibility claim — <screen> (<detail>); solving normally

A confirmed verdict now names the screen that proved it and the row, column, or bound it tripped on, rather than exiting silently:

Presolve: proved primal infeasible — empty equality row (equality row 7 is `0 = 3e0`)

Nothing that only reports is withheld from the re-derivation — empty rows, activity ranges, parallel rows, and emptied-row residuals all still apply — so no infeasibility presolve could detect before goes undetected now. What the guard costs, in the rare case it fires, is a handful of eliminations.

When the reduction is truncated

The reductions are iterated to a fixpoint — each one can expose work for the next, so presolve keeps going until nothing fires. It also carries a cap on how many layers that may take, and on a model with a long bound-propagation chain the cap is what stops it. When that happens the summary line says so:

Presolve: 315 → 128 vars, 233 → 77 rows (fixed 61, ..., tightened 158, cap-truncated after 32 layers)

This is common and it is not a problem. Measured across the LP and QP suites, the cap binds on 46% of LP models and 25% of QP models — and on every one of the 394 models that presolve at all, it changed only how tightly variable boxes were narrowed, never the structural reduction: same variables, same rows, same fixings, aggregations, forcing rows and dominated columns as running the iteration to convergence. Bound propagation is the one reduction that can keep going indefinitely, so it is what the cap ends up trimming.

What you get is still a correct problem — every reduction applied is a sound transform with its own dual recovery, and your solution is postsolved back to the original either way. The suffix is there so a reduction that came out of a truncated loop is distinguishable from one that converged, which matters when you are comparing two runs or reporting a bug against presolve. There is no option to turn it up.

Presolve is on by default. Turn it off with qp_presolve=no (e.g. to compare timings or isolate a solver issue):

pounce model.nl qp_presolve=no

Scope and limitations

• Convex problems only. Nonconvex (indefinite-Hessian) QPs, quadratic equalities, and quadratic inequalities whose feasible set is nonconvex are solved by the NLP path to a local minimum; POUNCE does not do global optimization.
• Convex QCQP (convex-quadratic constraints) routes to the conic IPM: each convex-quadratic inequality ½xᵀQx + aᵀx + b ≤ 0 (with Q ⪰ 0) is reformulated to one second-order cone (Q = FᵀF, so ‖Fx‖² = xᵀQx) and solved alongside the QP objective and linear constraints.

Both the primal solution and the constraint duals are written to the .sol file, in the same sign convention as POUNCE’s NLP path (so Pyomo and AMPL read them identically regardless of which solver ran).

Requests the convex path does not implement

The convex solvers are a specialized fast path, not a drop-in for every option the NLP path honors. Where a request would be dropped rather than merely unused, routing gives way rather than answering a different question:

A positive obj_scaling_factor is not in this table: it only rescales conditioning, and the convex path reports natural units either way, so both paths give the same answer.

When the convex path cannot certify an LP

Routing gives way one more time, and this one is decided after the solve rather than before it. Under auto, an LP whose convex solve finishes without a KKT certificate — Solved to acceptable level (reduced accuracy) or Maximum iterations exceeded — is re-solved on the general NLP interior-point path, which owns the whole verdict. Nothing from the declined convex solve is printed or written, so a rerouted run still reports exactly one status.

The case this exists for is the NETLIB gen / gen1 family. They are highly degenerate and rank-deficient, strict complementarity fails, and a pure interior-point method cannot certify the optimal vertex: the convex IPM spends its whole 200-iteration budget (190.8 s) and stops at a primal residual of 1.4e-7 against tol = 1e-8. The NLP filter-IPM — the same binary, the default for every other class — solves the same model in 19 iterations and 0.98 s to a strict certificate, matching Ipopt-3.14.20/MA57 to four figures. Rerouting is also the faster answer here: a second solve of one second is nothing against the three minutes the first one costs.

The fallback is narrow by construction, and does not fire when:

A tightened tol is deliberately not in that list: that is an accuracy request, so trying the engine that can meet it is the right response.

Infeasible and unbounded problems

The convex solver detects infeasibility and unboundedness directly, reporting a clean status instead of exhausting the iteration budget:

• Primal infeasible — no point satisfies the constraints. Reported with AMPL solve_result_num 200.
• Unbounded (dual infeasible) — the objective decreases without bound along a feasible direction. Reported with solve_result_num 300.

Each verdict is backed by a verified certificate (a Farkas infeasibility proof or an unbounded recession direction that is checked, not merely inferred), so these statuses are never reported in error; a problem the solver cannot certify simply runs to the iteration limit — and, if it is an LP under auto, is then handed to the NLP path (above).

solver_selection=qp-active-set follows the same contract. Its inner QP certifies the recession ray of the linearization, which on a nonlinear model is not yet a statement about the problem, so the ray is re-tested against the true objective and constraints before the 300 is reported; a ray that does not survive yields Search_Direction_Becomes_Too_Small, never an unboundedness claim.

The design and roadmap live in dev-notes/lp-qp-routing.md.

Convex Solver: LP, QP, and SOCP

POUNCE ships a specialized convex conic interior-point solver (pounce-convex) alongside the general NLP filter-IPM. It solves the standard-form convex program

minimize    ½ xᵀP x + cᵀx
subject to  A x = b
            G x ⪯_K h
            lb ≤ x ≤ ub

where P ⪰ 0 and the inequality block lies in a product cone K of nonnegative orthants and second-order cones. P = 0 is an LP; an all-orthant K is an LP/QP; second-order blocks make it an SOCP.

The method is a Mehrotra predictor–corrector primal–dual interior-point algorithm with Nesterov–Todd scaling for the cones, sharing the pure-Rust feral sparse LDLᵀ backend with the NLP path. It reaches optimality in materially fewer iterations than routing the same problem through the general NLP solver (≈30–50% fewer on bound/inequality QPs).

Inspiration. The conic interior-point design follows Clarabel (Goulart & Chen) — handling a quadratic objective directly and a product of symmetric cones — and the presolve follows PaPILO (the presolving library of SCIP). POUNCE does not wrap either (the pure-Rust guarantee) but ports their ideas; see Acknowledgments.

This chapter covers the Python API (pounce.qp and the differentiable pounce.jax layers). For automatic CLI/Pyomo routing of .nl LPs/QPs, see LP / QP Solver Routing. Runnable, progressive notebooks live in python/notebooks/: 15_convex_qp.ipynb, 16_socp.ipynb, 17_differentiable_convex.ipynb.

Quadratic programs

import numpy as np
from pounce.qp import solve_qp

# min ½·2‖x‖² − 3x₀ − 4x₁  s.t.  x₀ + x₁ ≤ 1,  0 ≤ x ≤ 1
r = solve_qp(
    P=np.diag([2.0, 2.0]),
    c=[-3.0, -4.0],
    G=[[1.0, 1.0]], h=[1.0],
    lb=[0, 0], ub=[1, 1],
)
r.status   # 'optimal'
r.x        # primal solution
r.y, r.z   # equality / inequality multipliers
r.z_lb, r.z_ub  # bound multipliers (≥ 0)
r.obj, r.iters

P (lower triangle used, assumed symmetric), A, and G accept dense arrays or scipy-sparse matrices; any of them may be omitted. The result is a QpResult dataclass with a .success property. The solver reports verified infeasibility / unboundedness ('primal_infeasible' / 'dual_infeasible') backed by a Farkas / recession certificate rather than an iteration-limit guess.

Second-order cone programs

A second-order (Lorentz) cone is { (t, x) : t ≥ ‖x‖₂ }. Partition the inequality rows of Gx ⪯_K h with cones — a list of (kind, dim) specs ("nonneg" or "soc"; a bare int means a second-order cone). Each slack block s = h − Gx must lie in its cone.

from pounce.qp import solve_socp

# minimize ‖x − x*‖  ⇔  min t s.t. (t, x − x*) ∈ SOC
r = solve_socp(
    c=[1.0, 0.0, 0.0],                 # minimize t
    G=-np.eye(3), h=[0.0, -2.0, 1.0],  # s = (t, x₀−2, x₁+1) ∈ SOC(3)
    cones=[("soc", 3)],
)
r.x   # ≈ [0, 2, -1]:  t* = 0, x = x*

Mixed cones compose — e.g. cones=[("nonneg", 1), ("soc", 2)] puts the first slack in ℝ₊ and the next two in a 2-D second-order cone. Large cones use a sparse diagonal-plus-rank-1 KKT representation (one auxiliary variable per cone, the ECOS/Clarabel “sparse SOC” trick) so the factorization stays sparse.

Warm starting

Feed a previous (or nearby) solution back to seed the interior-point iteration — useful for parametric sweeps, receding-horizon MPC, and branch-and-bound subproblems:

base = solve_qp(P=P, c=c, G=G, h=h, lb=lb, ub=ub)
nxt  = solve_qp(P=P, c=c2, G=G, h=h, lb=lb, ub=ub, warm_start=base)

The warm start only affects the iteration count, never the solution (a mismatch is ignored). The recentering is adaptive for the orthant (sized to the warm point’s KKT residual, so it exploits a nearby problem’s duals yet self-corrects when the active set moves) and re-centers the cone duals for second-order blocks (a converged conic point sits on the cone boundary, where the scaling is singular).

The step length is what makes it pay off

A warm start lowers the starting duality measure μ₀; whether that turns into fewer iterations depends on how much of each Newton step the solver is allowed to take. With a static fraction-to-boundary parameter τ, every step covers at most a τ fraction of the distance to the cone boundary, so μ falls by a fixed factor per iteration and the count is log₁/₍₁₋τ₎(μ₀/tol) however good the start was — a logarithm of the perturbation, not the one or two Newton steps a nearby problem deserves.

So on orthant blocks the step follows the Mehrotra tail τ = clamp(1 − μ, tau, tau_max): as the solve converges τ approaches 1 and a near-optimal iterate takes a near-full Newton step. On the QP families in the warm-start benchmark this is worth 35–60% of the warm iterations. Both ends are tunable, and both are method="ipm" only:

r = solve_qp(P=P, c=c2, G=G, h=h, warm_start=base,
             tau=0.95,        # floor: the flat τ far from the solution
             tau_max=0.999)   # ceiling on the tail (default: just under 1)

Passing tau_max=tau pins τ flat — the most conservative setting, and the one to reach for if a badly-conditioned sequence starts producing numerical_failure. Two scopes are deliberate and not tunable: second-order and PSD blocks always keep the static tau (their boundary is curved, and an iterate that close to it breaks the Nesterov–Todd scaling), and cold solves are unaffected because they run the homogeneous self-dual embedding, a different loop.

Batching and factorization reuse

from pounce.qp import solve_qp_batch, QpFactorization

# Solve many independent QPs in parallel (rayon, across instances).
results = solve_qp_batch([dict(P=P, c=c_k, G=G, h=h) for c_k in cs])

# Build the KKT symbolic factor once, solve many same-structure problems.
fac = QpFactorization(P=P, c=c0, G=G, h=h, lb=lb, ub=ub)
for c_k in cs:
    rk = fac.solve(P=P, c=c_k, G=G, h=h, lb=lb, ub=ub)  # reuses the factor

solve_qp_batch parallelizes across instances (outer-parallel / inner-serial) and QpFactorization reuses the AMD ordering and symbolic factorization across solves that share a structure — the two compose with warm starting.

Presolve (PaPILO-inspired)

Before the interior-point solve, POUNCE can apply a transaction-stack presolve with full primal and dual postsolve, modeled on PaPILO. The catalog:

• empty / duplicate / parallel (scalar-multiple) rows,
• fixed-variable elimination (singleton equalities),
• free columns and free-column singletons,
• activity-based redundancy and infeasibility detection,
• forcing constraints (a row at its activity extreme pins its variables),
• dominated columns (sign-definite columns optimal at a bound),
• bound tightening (domain propagation), with the active-bound multiplier re-attributed to its source row in postsolve,

iterated to a fixpoint so reductions cascade. Each reduction carries the data to reverse itself, and the postsolve reconstructs a valid KKT point of the original problem — the dual recovery is the contract, and is verified by KKT-residual tests. A cone-aware variant (presolve_conic) gates the ≤-row reductions off second-order-cone blocks (which are coupled) and recovers the reduced cone partition.

The iteration also carries a layer cap, and on a model with a long bound-propagation chain — commonly, on roughly half the LP corpus — the cap is what stops it rather than the fixpoint. That distinction is visible: presolve reports which of the two happened and the CLI says so on its summary line (see LP / QP Solver Routing). A truncated reduction is still correct — every reduction it did apply is a sound transform with its own dual recovery — and measured across the LP and QP suites the truncation costs only box tightness, never a structural reduction.

Presolve is applied automatically on the CLI LP/QP route; it lives in pounce-convex::presolve for Rust callers. See LP / QP Solver Routing.

Differentiable convex layers (JAX)

pounce.jax exposes the solve as a differentiable JAX op via the implicit-function theorem on the KKT system at the optimum (Amos & Kolter, OptNet, 2017). The forward calls the solver; the backward is a single linear solve through the same KKT matrix.

import jax, jax.numpy as jnp
from pounce.jax import solve_qp, solve_socp, QpLayer

# x*(c) for a parametric QP, differentiable w.r.t. all of P, c, G, h, A, b.
def loss(c):
    x = solve_qp(P=P, c=c, G=G, h=h)
    return jnp.sum((x - target) ** 2)

grad_c = jax.grad(loss)(c0)        # exact gradient via implicit diff
J = jax.jacrev(lambda c: solve_qp(P=P, c=c, G=G, h=h))(c0)

• Gradients are provided w.r.t. every parameter that enters through the optimum: c, b, h, and the matrices P, G, A (the full OptNet matrix derivatives; ∇P is the symmetric gradient).
• solve_socp differentiates SOCPs too — the complementarity row uses the cones’ arrow operators in place of the orthant’s diagonal.
• QpLayer captures a fixed P/G/A structure for use inside a larger JAX model, with jax.grad / jacrev / vmap and a parallel .batch.
• A warm start may be passed through (non-differentiated — it cannot change the solution or its gradients, only the iteration count).

All gradients are validated against finite differences in the test suite.

Global Optimization

Most of POUNCE settles a problem at a local optimum (the NLP filter-IPM and SQP) or exploits convexity so that local is global (the convex/conic IPM). For a genuinely nonconvex problem, POUNCE offers one certified-global route, and it is for polynomials:

• The SOS / Lasserre hierarchy (pounce-convex) — for polynomial problems, via a single semidefinite program. Callable from Rust (sos_minimize) and Python (pounce.sos_minimize).

It returns a result that is certified: a lower bound together with a moment certificate that, when exact, pins the global minimum and recovers its minimizer(s).

There is no general-purpose spatial branch-and-bound solver in POUNCE. For a nonconvex problem that is not polynomial — anything with exp/ln/trig — POUNCE has no certified-global path. Use the local NLP solver from several starting points (see the multistart notebooks below), or reformulate into the convex cone library. A pounce-global crate was prototyped and removed from main before release; its design is recorded in dev-notes/spatial-bnb-design.md.

The SOS / Lasserre path (polynomials)

When the objective and constraints are polynomials, the sum-of-squares / moment approach in pounce-convex certifies the global minimum from a single semidefinite program — no branching — by searching for the largest γ such that p(x) − γ lies in the Putinar cone (a sum of squares plus constraint multipliers). The SDP is solved by POUNCE’s own convex conic interior-point method; flat truncation of the resulting moment matrix certifies when the bound is exact, and a facial-reduction step recovers every global minimizer — even when the optimum is attained at several points.

From Python, a polynomial is a dict mapping an exponent tuple to its coefficient (the all-zeros key is the constant term):

from pounce.sos import sos_minimize

# x**4 - 2 x**2 + 3  ->  global minimum 2, attained at BOTH x = +1 and x = -1
r = sos_minimize({(4,): 1.0, (2,): -2.0, (0,): 3.0})
r.lower_bound       # ≈ 2.0
r.is_exact          # True — flat-truncation certificate: the bound is the minimum
r.minimizers        # both x = +1 and x = -1

Constraints are polynomials too, passed as inequalities (g_i(x) ≥ 0) and equalities (h_j(x) = 0); raise the relaxation order to tighten the bound (the Lasserre hierarchy) at the cost of a larger SDP. A runnable walkthrough — double well, a constrained problem, and a 2-D example — is in 18_sos_global_optimization.ipynb.

The same solver from Rust, via the pounce-rs facade with the convex feature on (pounce-rs = { version = "0.9", features = ["convex"] }):

#![allow(unused)]
fn main() {
use pounce_rs::convex::{sos_minimize, PolyProblem, Polynomial};
use pounce_rs::linsol::backend;      // the sparse LDLᵀ factory the solver takes

// x⁴ − 2x² + 3 → global minimum 2 at x = ±1.
let p = Polynomial::new(1, vec![(vec![4], 1.0), (vec![2], -2.0), (vec![0], 3.0)]);
let sol = sos_minimize(&PolyProblem::new(p), None, backend);
// sol.lower_bound ≈ 2; when the moment matrix is flat, sol.minimizers holds
// the global minimizer(s) — here both x = +1 and x = −1.
}

The full treatment lives in the pounce_convex::sos module documentation — reachable without a second dependency, since pounce_rs::convex re-exports the pounce_convex crate itself for anything outside its curated surface.

When SOS fits: polynomials of modest degree and dimension — one SDP, recovers all global minimizers, but the SDP grows with the relaxation order.

When SOS does not fit

For a general factorable problem (exp/ln/trig), or a polynomial whose SDP would be too large, the textbook tool is spatial branch-and-bound — and POUNCE does not have one. Two things you can do:

• Reformulate into the cone library. If the model can be cast as an LP, convex QP, SOCP, or an exponential / power / PSD cone program, local is global and the guarantee comes for free. See Choosing a Solver.
• Multistart the local solver. Running the NLP filter-IPM from many starting points finds the low minima in practice, but certifies nothing — there is no bound proving you have the global one. Three notebooks work through the tactics: repulsion-based sampling (19_find_minima_repulsion.ipynb), random restarts (20_find_minima_restart.ipynb), and basin hopping (21_find_minima_hopping.ipynb).

Solution Output

The .sol file

Following the AMPL solver convention, solving a positional .nl file writes a sibling <stub>.sol next to it — pounce problem.nl produces problem.sol. The file carries the primal x and dual lambda blocks plus an objno line with the AMPL solve_result_num, so AMPL (or any .sol reader) can pull the solution back:

pounce problem.nl                       # writes problem.sol
pounce problem.nl --sol-output out.sol  # write to an explicit path
pounce problem.nl --no-sol              # skip the .sol write

A .sol is written even when the solve fails, so the solve_result_num is always recoverable. Built-in problems (--problem …) have no .nl stub, so they only produce a .sol when --sol-output is given explicitly.

Reading solve_result_num

The objno line carries an AMPL solve_result_num (Gay 2005, Hooking Your Solver to AMPL §5). Consumers key on the band, not the exact number:

Pyomo maps each band to a TerminationCondition, so anything in 200–299 arrives as TerminationCondition.infeasible.

Infeasible: proved vs. local

Within the infeasible band POUNCE distinguishes how it knows:

The difference is real, not cosmetic. 201 is a structural detection made on the model’s bounds before iterating, not a certified proof — it is subject to the same floating-point limits as any interval computation, and is withheld whenever the violation is smaller than the feasibility tolerance. 200 is different in kind — on a nonconvex problem a positive local minimum of the violation does not rule out a feasible point elsewhere, which is why the console message says “Problem may be infeasible.”

Because 200 is an inference rather than a proof, it is withdrawn when POUNCE holds a point that contradicts it. Before any numerical path reports 200, the model’s own starting point is evaluated against every constraint; if it satisfies them all, the feasible set is demonstrably non-empty and the verdict becomes Error_In_Step_Computation (500) — an honest “the solve broke down” rather than a wrong answer. This can only ever withdraw a verdict: a model with no feasible point cannot produce such a point, so a correct 200 is unaffected. Supplying a feasible starting point is therefore worth doing on a model you believe is feasible but POUNCE reports otherwise.

When the region is found empty the solve is skipped entirely and the message names how it was found, so the claim is checkable:

POUNCE 0.9.0: InfeasibleProblemDetected (detected by presolve: bound propagation)
objno 0 201

201 requires presolve to be enabled (presolve=yes); it is off by default. A presolve-derived infeasibility is only reported when the contradiction holds on the original box — one produced by presolve’s own auxiliary elimination is re-checked after rollback and never certified.

One more route to 200: over-determined systems

An over-determined model — more equality rows than free variables, such as x == 0.2 with x == 0.8 — cannot be solved at all: it fails a structural gate before the first iteration. That used to be reported as Not_Enough_Degrees_Of_Freedom (504, the failure band), which says “cannot attempt this” for a model whose answer is already decided.

POUNCE now checks such a model for a bound-propagation contradiction on that failure path and reports 200 when it finds one. This does not need presolve=yes — nothing is transformed and no solve runs through the check — so it is the one way to reach the infeasible band with the default options and no iterations. A consistent over-determined system is unaffected and still reports 504.

Because the solve provably cannot run here, this route measures constraint residuals against each row’s declared magnitude rather than an absolute tolerance, so the verdict does not change when every row is multiplied by a constant. Elsewhere — wherever a solve can run — an infeasibility smaller than the feasibility tolerance is still withheld, as described above.

Choosing an output format

The .sol and JSON outputs are not exclusive — you can request both in the same run.

JSON Solve Report

Pass --json-output PATH to write a structured solve report alongside the regular console output:

pounce problem.nl --json-output result.json
pounce problem.nl --json-output result.json --json-detail full

The report carries everything an AMPL .sol file holds — status, primal x, dual lambda, suffix blocks — plus FAIR-aligned provenance metadata (Wilkinson et al. 2016, DOI 10.1038/sdata.2016.18) and, optionally, the per-iteration trajectory.

Detail levels

Choose summary for production logs and batch runs; full for debugging — it is the JSON equivalent of upstream’s print_level=8.

Reproducibility: recorded environment overrides

Solve-affecting environment variables — the POUNCE_FERAL_* linear-solver knobs and the legacy FERAL_PIVTOL / FERAL_PARALLEL — are captured into fair_metadata.environment when set, so a run that differs because one was exported in a shell profile says so instead of differing silently:

"environment": [
  { "name": "POUNCE_FERAL_PIVTOL", "value": "1e-6" }
]

The block is omitted entirely when no such variable is set (the common case). Debug-only gates (POUNCE_DBG_*) are not captured. See the schema reference for the full field contract.

Schema stability

The schema is versioned (pounce.solve-report/v1) so downstream tooling can pin against a major version:

• Adding fields is non-breaking — consumers must tolerate unknown fields.
• Removing or renaming a field bumps the major version (v1 → v2).

The Schema v1 Reference documents every field, the FAIR mapping, and the versioning policy in full.

POUNCE solve-report schema, v1

Schema tag: pounce.solve-report/v1

This document is the canonical reference for the JSON solve report emitted by pounce --json-output PATH and pounce_sens --json-output PATH. The report carries everything an AMPL .sol file holds — status, primal x, dual lambda, suffix blocks — plus FAIR-aligned provenance metadata and (optionally) the per-iteration trajectory.

Implementation: the serde structs live in crates/pounce-solve-report/src/lib.rs (per-iteration IterRecord in crates/pounce-nlp/src/solve_statistics.rs); crates/pounce-cli/src/solve_report.rs wires them to the CLI.

Why a structured solve report?

Production NLP workflows often need to (a) capture which solve produced which numbers for audit / reproducibility, (b) feed solver output into downstream tooling (notebooks, dashboards, ML pipelines) that don’t want to parse a free-form .sol file, and (c) compare runs across versions of pounce. Both upstream Ipopt’s stdout summary and AMPL’s .sol were designed for human consumption and AMPL’s reader respectively — neither carries provenance metadata, neither is schema-versioned, and neither is trivially machine-parseable across ecosystems.

A versioned JSON schema with FAIR-aligned provenance solves all three.

FAIR alignment

The fair_metadata block maps onto the four FAIR principles (Wilkinson et al. 2016, “The FAIR Guiding Principles for scientific data management and stewardship”, Scientific Data 3, 160018, DOI 10.1038/sdata.2016.18; citation verified via Crossref on 2026-05-14):

Versioning policy

schema is the version tag. Compatibility rules:

• Adding fields is non-breaking. Consumers MUST tolerate unknown fields. New optional fields land between versions; the major version doesn’t bump.
• Removing or renaming fields bumps the major version (v1 → v2). Consumers should pin against a major version (schema starts_with "pounce.solve-report/v1").
• Changing field semantics without a rename is forbidden. If semantics need to change, add a new field and deprecate the old.

The pre-1.0 phase of POUNCE itself does NOT relax this rule for the schema. Once a solve-report version ships, its field set is frozen even while the rest of the solver is under churn.

Top-level shape

{
  "schema": "pounce.solve-report/v1",
  "fair_metadata": { ... },
  "problem":       { ... },
  "solution":      { ... },
  "statistics":    { ... },
  "iterations":    [ ... ],  // optional, omitted when empty
  "linear_solver": { ... }   // optional, omitted when backend did not report
}

Fields

schema (string, required)

Identifier for this schema version. Always "pounce.solve-report/v1" for v1. Major-version bumps change the prefix; minor / patch (additive) changes do not.

fair_metadata (object, required)

solver sub-object

Input descriptor (input)

Tagged enum keyed on kind. Possible shapes:

{ "kind": "nl-file", "path": "/path/to/foo.nl", "size_bytes": 366 }
{ "kind": "builtin", "name": "rosenbrock" }
{ "kind": "tnlp-direct" }

• nl-file — the input came from .nl file at path. size_bytes is present when the file’s metadata is readable; consumers that want bit-exact provenance can hash the file themselves.
• builtin — the input was a built-in problem named by name (e.g. pounce --problem rosenbrock).
• tnlp-direct — used by library callers building a TNLP in-process without a .nl round-trip.

Environment overrides (environment)

An array of { "name", "value" } objects, one per solve-affecting environment variable set in the process at report time:

"environment": [
  { "name": "POUNCE_FERAL_PIVTOL", "value": "1e-6" }
]

The whole array is omitted when none are set (the common case). Only the variables that change pounce’s numerics or parallelism are captured — the POUNCE_FERAL_* linear-solver knobs and the legacy FERAL_PIVTOL / FERAL_PARALLEL. These alter the factorization and can otherwise silently differ a run between two machines (e.g. one with POUNCE_FERAL_PIVTOL exported in a shell profile) with nothing in the report saying so. The POUNCE_DBG_* debug gates are deliberately not captured — they only add diagnostic output and never change the result.

Presence records that the variable was set, not that it took effect: an explicit OptionsList setting (e.g. feral_pivtol in an options file) takes precedence over the env fallback. See Options › Environment overrides for the option each variable maps to.

problem (object, required)

Problem dimensions reported by the TNLP at get_nlp_info().

solution (object, required)

Suffix entries

{
  "name": "sens_sol_state_1",
  "target": "var",
  "kind": "real",
  "values": [0.576..., 0.378..., -0.046..., 4.5, 1.0]
}

statistics (object, required)

Projection of pounce_nlp::solve_statistics::SolveStatistics minus the per-iteration history (which lives at the top level when present).

null values. The four residuals are produced by the convergence check at the end of a solve. A solve the solver refused — rejected during setup (NotEnoughDegreesOfFreedom, InvalidProblemDefinition), aborted, or caught by the batch panic handler — never reaches it, and these slots are emitted as null rather than 0.0. A zero there is indistinguishable from a perfect solve, and consumers acted on it: it was enough to make pounce.minimize report success=True for a problem the solver had declined to attempt.

The two objective fields follow the same rule for the same reason: 0.0 is an ordinary objective value, so it cannot signal “never evaluated”. They are seeded from the current iterate whenever one exists, so they are null only when the solve produced no point at all.

Consumers should treat null as “not computed”, not as zero. pounce’s own readers map it to NaN, which fails closed against any value <= tol test. | num_obj_evals | integer | eval_f call count. | | num_constr_evals | integer | eval_g call count. | | num_obj_grad_evals | integer | eval_grad_f count. | | num_constr_jac_evals | integer | eval_jac_g count. | | num_hess_evals | integer | eval_h count. | | total_wallclock_time_secs | float | Wall time spent inside optimize_*. | | restoration_calls | integer | Number of restoration-phase entries (pounce#12). | | restoration_inner_iters | integer | Cumulative inner-IPM iterations across all restoration calls. | | restoration_outer_iters | integer | Outer iterations that ran in restoration mode (R-line equivalents). | | restoration_wall_secs | float | Wall time spent inside perform_restoration. |

Eval counters (num_*_evals) populate only on the .nl-file path because the pounce binary’s CountingTnlp wrapper tracks them. Library callers using IpoptApplication::optimize_tnlp directly see zeros there; the underlying counts are still available through upstream’s IpoptCalculatedQuantities if needed.

iterations (array of object, optional)

Per-iteration trajectory. Emitted only at --json-detail full (when IpoptApplication::enable_iter_history() was called). Omitted from JSON entirely when empty.

Each row maps to one line of the upstream-formatted console iter table. Fields:

linear_solver (object, optional)

Aggregate post-mortem from the symmetric-indefinite linear backend that solved the KKT systems. Populated only when the backend self-instruments (the default FERAL backend does; HSL MA57 and custom backends plugged through set_linear_backend_factory do not). Omitted from JSON when no backend reported.

Detail levels

The --json-detail LEVEL flag selects how much detail is emitted. Levels map to verbosity in the same spirit as upstream’s print_level (0 silent → 12 maximum debug):

summary is the right choice for production logs and batch runs. full is the debugging equivalent of upstream’s print_level=8.

Worked example

pounce_sens crates/pounce-cli/tests/fixtures/parametric.nl out.sol --json-output result.json --json-detail full produces (truncated for brevity):

{
  "schema": "pounce.solve-report/v1",
  "fair_metadata": {
    "result_id": "1778777029606881000-76543",
    "created_at_iso": "2026-05-14T16:43:49.606Z",
    "created_at_unix_nanos": 1778777029606881000,
    "elapsed_seconds": 0.011,
    "solver": {
      "name": "pounce",
      "version": "0.1.0",
      "target_triple": "x86_64-apple-darwin"
    },
    "license": "EPL-2.0",
    "input": {
      "kind": "nl-file",
      "path": "crates/pounce-cli/tests/fixtures/parametric.nl",
      "size_bytes": 366
    }
  },
  "problem": { "n_variables": 5, "n_constraints": 4, "n_objectives": 1, "minimize": true },
  "solution": {
    "status": "SolveSucceeded",
    "solve_result_num": 0,
    "objective": 0.5510204081632656,
    "x":      [0.6326530575201161, 0.3877551079678144, 0.020408165487930466, 5.0, 1.0],
    "lambda": [-0.16326530000405073, -0.28571431357898697, -0.16326530000405073, 0.18075803406303625],
    "suffixes": [{
      "name": "sens_sol_state_1",
      "target": "var",
      "kind": "real",
      "values": [0.5765305974643309, 0.3775510440570709, -0.04591835847859835, 4.5, 1.0]
    }]
  },
  "statistics": { "iteration_count": 9, "final_dual_inf": 2.89e-14, "...": "..." },
  "iterations": [
    { "iter": 0, "objective": 0.0451, "inf_pr": 5.0, "inf_du": 0.407, "mu": 0.1,
      "d_norm": 0.0, "regularization": 0.0, "alpha_dual": 0.0, "alpha_primal": 0.0,
      "alpha_primal_char": " ", "ls_trials": 0 },
    { "iter": 1, "objective": 0.957, "inf_pr": 0.212, "...": "..." }
  ]
}

Consumer guidance

• Pin the major version. Check schema.startswith("pounce.solve-report/v1") before consuming.
• Tolerate unknown fields. New optional fields will land between minor versions of pounce. Use serde(default) / equivalent.
• Distinguish “no solve” from “solve produced zero”. Pre-solve, scalar fields are 0.0 (not NaN, because JSON has no NaN literal). statistics.iteration_count == 0 is the signal that no solve occurred.
• solution.x / solution.lambda may be empty. When the binary couldn’t capture the final iterate (currently: the pounce binary on its newton_driver fast-path for m=0, n≤1000 problems), the arrays are empty and the keys are omitted from JSON entirely. pounce_sens always populates them.

References

• Wilkinson et al. (2016). “The FAIR Guiding Principles for scientific data management and stewardship.” Scientific Data 3, 160018. DOI 10.1038/sdata.2016.18. (Verified via Crossref 2026-05-14.)
• Gay (2005). “Hooking Your Solver to AMPL.” https://ampl.com/REFS/hooking2.pdf. §5 (Returning Results to AMPL) for the .sol baseline this schema is structured around.
• SPDX license identifiers: https://spdx.org/licenses/.

Verifying Solutions

pounce verify <problem.nl> <claim.sol> [OPTIONS]

pounce verify independently checks that the solution in a .sol file actually satisfies the constraints and bounds of a .nl problem. It re-derives feasibility from the model itself — it does not trust the .sol’s status line, and it does not rerun the solver. This makes it the trust anchor when pounce is a tool an agent calls: the agent proposes a solution, and a small, deterministic checker disposes.

Optimization is unusually well-suited to this because a solution is far cheaper to verify than to produce: a claimed x* is just numbers, and feasibility is a single constraint evaluation — g_l ≤ g(x*) ≤ g_u, x_l ≤ x* ≤ x_u — O(nnz) work with no resolve and no dense linear algebra.

Status. The verify check itself — recompute feasibility against the canonical model, with a content-addressed receipt — is solid and ready to use; it needs no secrets and is the recommended default. The signing and remote-service trust layer layered on top of it (HMAC receipts, the signer_service.py reference, running the MCP server as a remote authority) is a proof of concept: it demonstrates the architecture but is not hardened for production. If you want to rely on the signed/remote path for real, see Status and hardening at the end for the checklist of what that would take.

What it defends against

In an agent workflow, three things can go wrong with “here is a solution”:

The key design rule: always verify against the canonical problem, never against whatever the agent claims it solved. If the agent loosened a bound to manufacture feasibility, the returned x* still violates the canonical bound, and verify reports it.

Output and exit codes

$ pounce verify gaslib40_steady.nl good.sol
pounce verify — independent solution check
  problem : gaslib40_steady.nl  (1694 vars, 1682 cons)
            sha256:4bb435a3…
  solution: good.sol
            sha256:b77d9e7b…
  claimed solve_result_num: 0

  feasibility (tol 1.0e-6):
    max constraint violation: 1.407e-12  at c[114] (value 1.4e-12, bounds [0, 0])
    max bound violation     : 9.775e-9   at x[24]  (value 1.05, bounds [1.05, 2.0])
  objective at x*: 1.2899875310e0

  optimality (tol 1.0e-6, duals + bound multipliers supplied):
    KKT stationarity residual (bound-projected)  : 2.675e-3  (dual sign +1)
    dual infeasibility (with z_L/z_U suffixes)   : 6.248e-14
    constraint complementarity (rows, |λ|·slack) : 0.000e0
    bound complementarity (vars, |z|·slack)      : 9.091e-10

  VERDICT: VERIFIED — solution is feasible for the canonical problem

A consumer (CI step, agent harness, Makefile) gates on the exit code.

Options

Feasibility is the gate; optimality is reported

By default only feasibility gates the exit code. Feasibility is rigorous and sign-convention-independent — it is the guarantee that matters when the claim is “this solution meets the constraints.”

When the .sol carries constraint duals, verify also reports a KKT stationarity residual (the bound-projected “dual infeasibility”: the part of ∇f + Jᵀλ that a valid sign-constrained bound multiplier cannot absorb) and the complementarity residuals below. These are informational unless you pass --require-optimal. The AMPL dual-sign convention can differ from pounce’s, so verify computes the residual for both signs and reports the better one plus the sign it used.

The two complementarity residuals

Two distinct quantities answer to “complementarity”, and at the same point they can disagree by many orders of magnitude. verify names both by what they range over, and never prints an unqualified complementarity residual:

The bound one is what a solver prints as Complementarity — Ipopt’s and pounce’s own end-of-solve report alike. Do not compare a solver’s Complementarity against the row line; they measure different things and neither is wrong for what it measures.

When the .sol carries no ipopt_zL_out / ipopt_zU_out suffixes, bound complementarity is reported as not checked, never as 0.0.

Bound multipliers sharpen the stationarity check

Without the suffixes, z_L and z_U are inferred from which bounds are active, and the reported stationarity residual is bound-projected: it projects out exactly the component a bound multiplier would carry, so it reads 0.0 on a point whose bound multiplier is missing or wrong.

When the .sol does carry them (pounce always writes them; so does Ipopt’s AMPL interface), verify additionally reports the exact dual infeasibility ‖∇f + Jᵀλ − (zL_out + zU_out)‖∞, using the multipliers the file actually claims. That number is directly comparable to the solver’s Dual infeasibility, and it is what --require-optimal gates on whenever it is available — the projected residual can only understate it.

The JSON receipt

--json-output writes a machine-readable receipt that content-addresses both inputs by SHA-256 — so a downstream consumer can confirm exactly which problem and which solution were checked:

{
  "pounce_verify_version": 1,
  "solver": "pounce 0.4.0",
  "problem":  { "path": "…", "sha256": "4bb435a3…", "n_vars": 1694, "n_cons": 1682 },
  "solution": { "path": "…", "sha256": "b77d9e7b…", "duals_present": true },
  "tolerances": { "feasibility": 1e-6, "optimality": 1e-6 },
  "feasibility": {
    "max_constraint_violation": 1.4e-12,
    "worst_constraint": { "index": 114, "name": "c[114]", "value": 1.4e-12,
                          "lower": 0.0, "upper": 0.0, "violation": 1.4e-12 },
    "max_bound_violation": 9.77e-9,
    "worst_bound": { "index": 24, "name": "x[24]", … },
    "feasible": true
  },
  "optimality": {
    "available": true,
    "stationarity_residual": 2.6e-3,
    "stationarity_residual_with_bound_multipliers": 6.2e-14,
    "constraint_complementarity_residual": 0.0,
    "bound_complementarity_residual": 9.1e-10,
    "bound_multipliers_present": true,
    "complementarity_residual": 0.0,
    "optimal": true, "note": "…"
  },
  "verdict": "VERIFIED",
  "verified": true
}

bound_complementarity_residual and stationarity_residual_with_bound_multipliers are null when the .sol carries no ipopt_zL_out / ipopt_zU_out suffixes — null means not checked, not zero. complementarity_residual is a deprecated alias of constraint_complementarity_residual, kept so v1 consumers keep parsing; its bare name is the one that invited the wrong comparison, so read the qualified field instead.

A consumer should accept a solution iff:

1. verified == true, and
2. problem.sha256 equals the SHA-256 of its own canonical .nl, and
3. (when signing is used) the signature validates — see below.

Checking the hash in step 2 is what closes the “solved the wrong problem” gap at the receipt layer: the receipt is only meaningful for the exact problem bytes it names.

The default: recompute, don’t trust a receipt

The strongest and simplest design uses no key and no signature at all: the consumer runs pounce verify itself, against its own copy of the canonical .nl.

# the consumer does this — not the agent
pounce verify ./canonical/problem.nl ./from-agent/claim.sol || reject

Because verification is keyless, deterministic, and cheap (O(nnz), no resolve), the consumer can afford to just do it rather than trust someone else’s word. In this design the agent is never in the trust path: it hands over x*, and the consumer believes its own arithmetic. There is no key to steal, so the question “what if the agent gets the key?” does not arise. Forgery is impossible because nothing is being trusted on faith — feasibility is decided by evaluating g(x*), not by matching fields in a document.

This is the recommended default. Prefer it whenever the consumer can run pounce verify (or call a verifier it controls). Reach for signatures only when it genuinely cannot — see below.

Signed receipts — trust transport, conditional on key isolation

Signing addresses a narrower situation: the consumer won’t or can’t recompute — a remote or expensive verifier, or an audit log you want to trust later without re-solving — and instead wants to trust a receipt produced elsewhere. A signature lets that receipt be checked without redoing the work.

When the POUNCE_VERIFY_KEY environment variable is set (non-empty), the receipt gains:

"signature_alg": "HMAC-SHA256",
"signed_fields": ["verify_version","nl_sha256","sol_sha256",
                  "n_vars","n_cons","feasible","verified","verdict"],
"signature": "5bdcc146bf60754e…"

The signature is HMAC-SHA256(key, preimage), where preimage is a deliberately float-free byte string — only hex hashes, integer counts, and the verdict — so any language reproduces it byte-for-byte without float-formatting parity problems. The exact preimage is:

pounce-verify-receipt/v1
verify_version=1
nl_sha256=<hex>
sol_sha256=<hex>
n_vars=<int>
n_cons=<int>
feasible=<true|false>
verified=<true|false>
verdict=<VERIFIED|REJECTED>

(eight lines, \n-joined, with a trailing newline; booleans lowercase.) A holder of the key recomputes the HMAC over this preimage and compares it to signature.

What the signature does and does not guarantee

HMAC gives existential unforgeability under chosen-message attack — but only while the key stays secret. That single condition carries the entire weight of the scheme:

• A party without the key cannot mint a receipt that validates. It cannot produce a correct MAC for any message it has not already seen one for. This is the real guarantee.
• A party with the key can sign anything — including a hand-written receipt that falsely says verified: true. If the agent obtains the key, signing buys you nothing, and you must fall back to recomputation.

So the security reduces to a question that is not cryptographic: can the agent reach the key? If the keyholder and the agent share a process, a user account, an environment, or a filesystem, the agent simply reads it — echo "$POUNCE_VERIFY_KEY", cat keyfile, or introspecting shared memory — and the signature becomes theater. An agent with a shell on the same user/host as the signer is not kept out by HMAC.

Setting POUNCE_VERIFY_KEY for a pounce verify that runs in the agent’s own environment is the weakest form and provides no real protection against that agent. It is only meaningful when the key lives behind a privilege boundary the agent cannot cross. See Out-of-process signing below.

Two further limits HMAC does not address, even with a perfectly isolated key:

• Replay. A previously, legitimately-signed receipt can be re-presented. Here that is benign: a receipt attests “this specific x* (by sol_sha256) is feasible for this specific problem (by nl_sha256),” and that statement stays true. Replay cannot manufacture a valid-looking attestation for a new (problem, solution) pair the signer never blessed. If you ever need freshness (“the work was done this time”), add a nonce/timestamp to the signed preimage — but for “does this solution meet the constraints,” freshness is irrelevant.
• Wrong model. The signature says nothing about whether the .nl encodes the problem you meant. That spec-ownership question must live outside the agent regardless of any crypto. verify certifies feasibility against a given model; it does not certify that the model is the right one. Likewise, feasibility is fully checkable but global optimality is not — the stationarity residual certifies a first-order/KKT point, not a global minimum.

Out-of-process signing

For a signature to mean anything against a capable agent, the key must live in a separate trust boundary — a different process and user/container/ host that the agent cannot exec into or read. The signer exposes “verify-and-sign” as an RPC; the key never enters the agent’s address space.

 ┌─ agent boundary (no key) ────────────┐      ┌─ signer boundary (holds key) ─────┐
 │ agent proposes x*                     │      │ POUNCE_VERIFY_KEY in THIS env only │
 │   │                                   │ RPC  │   1. run `pounce verify` on the    │
 │   └── POST /verify {nl, sol} ─────────┼─────▶│      CANONICAL .nl (+ the key)     │
 │                                       │      │   2. binary signs the receipt      │
 │   signed receipt ◀────────────────────┼──────┤   3. return receipt JSON           │
 └───────────────────────────────────────┘      └────────────────────────────────────┘
        │
        └── relays receipt to the consumer
consumer: accept iff  verified==true  ∧  problem.sha256==canonical  ∧  signature valid

What each party can do under this split:

The boundary is only real if the agent cannot run code as the signer’s user or on its host. Running the signer as a separate user, container, or host (or behind a KMS/HSM that signs without exposing the key) is what turns “signed” from theater into a guarantee. An MCP server is already a separate process from the model, which helps — but only achieves isolation if the agent also lacks a shell on the same user/host.

A minimal reference signer is in studio/mcp/examples/signer_service.py: a stdlib HTTP service that holds the key in its own environment, shells out to pounce verify, and returns the signed receipt. The agent calls it; the agent’s environment never contains the key.

Use in an agent workflow

Putting it together — recompute by default, sign only to transport trust:

agent ── proposes x* ──▶  consumer / verifier-it-controls
                            1. pin + hash the canonical .nl
                            2. pounce verify .nl .sol   (against the CANONICAL model)
                          ◀─ accept iff verified==true ∧ problem.sha256==canonical

When the verifier must be remote and the consumer won’t recompute, insert an out-of-process signer (above) and add ∧ signature valid to the consumer’s acceptance test — remembering that the last clause is only as strong as the signer’s key isolation.

The pounce-studio MCP server exposes verify_solution so an agent can request a check but cannot fake its result. Deploy that server as a distinct boundary from the agent (separate user/container) for the signature to carry weight; otherwise rely on the consumer recomputing.

Status and hardening

What is ready to use as-is:

• The feasibility check (pounce verify, and the consumer-recomputes pattern). It is deterministic, keyless, content-addressed, and rigorous — this is the part to build on.

What is a proof of concept — demonstrates the shape, not hardened:

• HMAC signing via POUNCE_VERIFY_KEY, the signer_service.py reference, and treating a remotely-deployed MCP server as a signing authority.

If you ever want to depend on the signed/remote path in production, these are the gaps to close. None are implemented here.

Key management

• Don’t keep the key in a plain environment variable or file. Use a KMS/HSM (or sealed secret) that signs without exposing the key to the process — then even a compromised signer host can’t exfiltrate it.
• Add key rotation and a key id in the receipt (kid) so a consumer knows which key to check against and old receipts stay verifiable across rotations.
• Consider an asymmetric scheme (e.g. Ed25519) instead of HMAC when more than one party must verify without also being able to sign — HMAC’s symmetric key means every verifier is also a forger. Public-key signatures give public verifiability with a single private signer.

Transport / service (the moment it leaves stdio)

• TLS on the endpoint; never plaintext for a service that holds a key.
• Authn/authz — bearer token or OAuth on every request (MCP’s HTTP transport supports this). An unauthenticated endpoint that runs solves and shells out is effectively remote code execution.
• Resource limits — request-size caps, solve timeouts (there is a timeout_seconds, but also wall/CPU/memory limits at the OS level), concurrency caps, and rate limiting.
• Sandbox the solve — treat every .nl as untrusted input. Parsing and evaluating an arbitrary model is attacker-controlled computation; run it in a locked-down container/user with no network and a constrained filesystem.

Input handling

• Over a network the path-based tools (nl_file/sol_file) assume a shared filesystem. Prefer content upload (the server receives and hashes the exact .nl/.sol bytes) so there’s no path-traversal surface and the receipt binds what was actually sent. The reference signer’s POUNCE_SIGNER_ROOT allowlist is a stopgap, not a substitute.

Freshness / replay

• The current preimage has no nonce or timestamp, so a signed receipt is replayable. That is benign for “is this x* feasible” (a timeless fact), but if a consumer needs “this was checked recently” or “in response to my request,” add a nonce/timestamp (and a receipt expiry) to the signed preimage and bump the pounce-verify-receipt version.

Auditability

• Log every verification (problem hash, solution hash, verdict, key id, caller identity) to an append-only store, so a disputed result can be reconstructed. Keep the key out of the logs.

Standing non-goals (true regardless of hardening)

• verify certifies feasibility against a given model — it does not certify the model is the right one. Model/spec correctness must be owned outside the agent.
• Feasibility is fully checkable; global optimality is not. The stationarity residual certifies a first-order/KKT point, not a global minimum.

Sensitivity Analysis

POUNCE includes a parametric sensitivity capability compatible with upstream Ipopt’s contrib/sIPOPT/ (Pirnay, López-Negrete & Biegler 2012, DOI 10.1007/s12532-012-0043-2). It computes the first-order change in the optimal primal solution with respect to a problem parameter, reusing the KKT factorization from the converged solve. Four entry points cover the common workflows.

AMPL CLI

The main pounce driver auto-detects the sIPOPT suffixes (sens_state_1, sens_state_value_1, sens_init_constr) in an input .nl, runs a post-optimal sensitivity step after the solve, and writes the perturbed primal back as a sens_sol_state_1 suffix — no separate binary or flag needed:

pounce problem.nl                   # writes problem.sol
pounce problem.nl out.sol --json-output result.json --json-detail full

pounce_sens is retained as a thin backward-compatibility alias: pounce_sens in.nl out.sol is identical to pounce in.nl out.sol, so existing AMPL / solver scripts keep working unchanged.

Related flags:

• --sens-boundcheck / --sens-bound-eps EPS — clamp the perturbed primal x* + Δx onto the declared [x_l, x_u] box.
• --compute-red-hessian / --rh-eigendecomp — compute the reduced Hessian (and its eigendecomposition) over the variables tagged by the red_hessian integer var-suffix.

Rust library

Reach the sensitivity path through the pounce-rs facade, with the sensitivity feature on:

[dependencies]
pounce-rs = { version = "0.9", features = ["sensitivity"] }

SensSolve is a builder that wraps the on_converged callback plumbing into a single call:

#![allow(unused)]
fn main() {
use pounce_rs::sensitivity::SensSolve;

let result = SensSolve::new(vec![2, 3])
    .with_deltas(vec![0.05, 0.0])
    .with_reduced_hessian()
    .run(&mut app, tnlp);
// result.dx, result.reduced_hessian, result.status
}

with_reduced_hessian_eigen() adds the eigendecomposition; with_boundcheck(eps) enables the bound projection.

Eigenvector sign convention

Every eigendecomposition POUNCE hands back — the reduced Hessian’s here and through the CLI and Python wrappers, the QP one from QpSensitivity.reduced_hessian, and covariance().eigen() / information().eigen() in pyomo-pounce — returns sign-pinned eigenvectors: the largest-magnitude component of each column is positive, ties broken by the earliest row. v and -v are equally valid eigenvectors, so without a convention the direction you read back depends on the arithmetic that produced it and is not reproducible across builds or machines.

The sign is all that is pinned. A repeated eigenvalue leaves the basis within its eigenspace arbitrary — any rotation of those columns diagonalizes equally well — so read a degenerate block as a subspace, not column by column.

Python

solve_with_sens exposes the same capability from the cyipopt-compatible Python wrapper:

# pin_constraint_indices is required; pass deltas=..., compute_reduced_hessian=True,
# or both. Returns (x, info) — sensitivity outputs live in the info dict.
x, info = prob.solve_with_sens(x0, pin_constraint_indices=[2, 3],
                               deltas=[0.05, 0.0], sens_boundcheck=True)
# info["dx"], info["reduced_hessian"], info["reduced_hessian_eigenvalues"], ...

compute_reduced_hessian=True returns the reduced Hessian in info["reduced_hessian"]; rh_eigendecomp=True adds its eigendecomposition; sens_bound_eps=… tunes the bound projection. See python/notebooks/04_sensitivity.ipynb for a walkthrough.

Pyomo

pyomo_pounce wraps the same machinery in a declare-then-query interface: flag the parameters that matter while building the model (no perturbed values required), solve normally, then ask for derivatives. Parameters are declared with declare_sens_param (mutable Param or fixed Var, scalar or indexed); when declarations are present, SolverFactory("pounce").solve(m) runs in-process and keeps the converged KKT factorization, so every query afterwards is a single backsolve.

import pyomo.environ as pyo
import pyomo_pounce
from pyomo_pounce import declare_sens_param, gradient, estimate

m.p = pyo.Param(initialize=2.0, mutable=True)
declare_sens_param(m.p)                 # a flag, not a perturbation

pyo.SolverFactory("pounce").solve(m)    # ordinary solve

gradient(m.x, wrt=m.p)                  # dx*/dp (float)
gradient(m.con, wrt=m.p)                # d(multiplier of con)/dp
G = gradient(m.z, wrt=m.r)              # containers -> Gradient object
G[m.z[1], m.r[2]]; G.to_dataframe()     # element access / full Jacobian
estimate(m, [(m.p, 2.5)])               # first-order solution estimate at
                                        # new values, clamped to bounds

gradient returns exact first-order derivatives (unit-perturbation backsolves, no finite differencing); estimate combines the stored derivative columns for arbitrary perturbed values after the fact. Its perturbation is measured from the solve point, not the Param’s current value, so writing a measurement into the Param before asking (the receding-horizon pattern) does not change the answer. It also warns when the linear step leaves the variable bounds (a single-pass projection analogous to the CLI’s --sens-boundcheck) — with one exception, a bound written on a declared Param, covered in Declared Params in variable bounds below. Multiplier sensitivities are available for equality constraints. Models without declarations solve through the ordinary AMPL/CLI path, unchanged. See python/notebooks/25_pyomo_sensitivity.ipynb for a worked optimal-control example (initial conditions as parameters; the first-move gradient IS the NMPC feedback gain).

Declared Params in variable bounds

A limit is often most naturally written as a bound rather than a constraint:

m.u_max = pyo.Param(initialize=1.0, mutable=True)
declare_sens_param(m.u_max)
m.u = pyo.Var(m.t, bounds=(0, m.u_max))   # the cap, as a bound

pyomo.contrib.sensitivity_toolbox, which supplies the expression surgery underneath, substitutes declared Params in constraint expressions only. A Param left in a bound is written to the .nl file as a constant at its pre-perturbation value, so the bound never moves and gradient(m.u[t], wrt=m.u_max) reads exactly 0.0 — a wrong answer that is indistinguishable from a legitimate insensitivity.

POUNCE rewrites such a bound as a constraint over the substituted variable before the solve, so both spellings of the same limit give the same derivative. Expression bounds work too, e.g. bounds=(0, 2 * m.p + 1). Two kinds of variable are deliberately left alone: fixed Vars, whose bounds the solver never enforces, and Vars on deactivated Blocks.

This is a deliberate divergence from sensitivity_calculation, which still reports zero for the same model. Four things follow from it:

• The bound is dropped on the clone that is solved. m.x.ub reads None there and the NL row carries the reader’s no-bound sentinel 1e19 — finite, so an isinf() test will not catch it. The model you wrote is never modified.
• estimate() does not clamp against a rewritten bound, and raises no clamp warning for it. That is correct rather than an oversight: the bound now moves with the perturbation, so the linear step already respects it to first order.
• covariance()’s bound-active projection still fires. The value the bound held at the solve point is recorded and read back for the activity test, so a declare_fitted variable capped by a declared Param is still projected and still warns.
• It costs a row. A simple bound is handled directly in the barrier; a general inequality costs a slack and a Jacobian row. A model with many Param-dependent bounds trades roughly one row per bound. Only models that write a bound in terms of a declared Param pay this.

Solver options and warm starts

Solver options reach the in-process path the same two ways they reach an ordinary solve: factory-level (SolverFactory("pounce", options={...}) or solver.options[...]) and per-call (solve(m, options={...})), with the per-call mapping winning on conflict. Everything the CLI accepts works here: tolerances, max_iter, scaling, warm-start knobs.

With warm_start_init_point=yes (Python True works too) among the options, the initial multipliers come from the model’s suffixes, the same ones the ASL path uses: dual for equality multipliers, ipopt_zL_in / ipopt_zU_in for bound multipliers, matched by component name (a constraint rewritten by the declared-parameter surgery is reached through its internal alias). Sign conventions are handled: dual holds the AMPL marginal and ipopt_zU_in Ipopt’s negative-at-upper value, and both are translated to the solver’s internal conventions on the way in.

One deliberate improvement over the ASL path: entries you do not supply take the solver’s own default initialization rather than zero. Through a dense ASL array an absent entry is indistinguishable from a zero multiplier, and a zero bound multiplier on an active bound is a contradictory KKT certificate the solver must first recover from. A suffix knows which entries exist, so an explicit zero is honored (then floored at warm_start_mult_bound_push, exactly as a round-tripped inactive multiplier is) and absence means “initialize as you normally would”: the solver’s own bound_mult_init_val for bound multipliers, and for equality duals the warm path’s 0, which is not the cold path’s least-squares estimate. Seed everything from a prior solve and the two paths behave identically; seed partially and the in-process path degrades gracefully.

Watching the solve (tee=True)

SolverFactory("pounce").solve(m, tee=True) streams the solver’s full Ipopt-style log — banner, problem statistics, iteration table, and end-of-run summary — live to standard output, including inside a Jupyter notebook cell. The log is emitted by the engine itself (the same blocks the pounce CLI prints), so the in-process path just tails it: a long solve shows its iteration table as it runs rather than as one block at the end. Without tee=True the solve is silent, matching the Pyomo convention.

Parameter covariance and identifiability

For a parameter-estimation model whose objective is a plain sum of squared residuals, the factorization from ONE ordinary solve yields the asymptotic covariance of the fitted parameters. Declare the fitted variables (they stay free) and the residual container while building the model, solve, and ask:

from pyomo_pounce import covariance, declare_fitted, declare_residual

m.A = pyo.Var(); m.k = pyo.Var()        # the fitted parameters, free
declare_fitted(m.A, m.k)

m.r = pyo.Var(m.I)                      # residuals, one per data point
m.res = pyo.Constraint(m.I, rule=...)   # r[i] == y[i] - model(A, k, t[i])
declare_residual(m.r)

m.obj = pyo.Objective(expr=sum(m.r[i]**2 for i in m.I))
pyo.SolverFactory("pounce").solve(m)    # one solve

cov = covariance(m)                     # no further information needed

cov[m.A, m.k]               # covariance entry (either order)
cov.std_err[m.k]            # standard error of one parameter
cov.correlation[m.A, m.k]   # correlation matrix entry
cov.matrix                  # dense numpy array, ordered like cov.params
w, V = cov.eigen()          # eigendecomposition, for identifiability

The recipe: the parameter block of the inverse KKT matrix, one backsolve per parameter against the held factor, equals the inverse reduced Hessian of the eliminated problem, and for a sum-of-squares objective cov = 2 sigma^2 (K^-1)_pp. The factor 2 belongs to the unscaled sum of squares (a Gaussian negative log-likelihood objective, SSR / (2 sigma^2), would drop it). The scaling is pinned by test against the analytical linear-regression covariance sigma^2 inv(X^T X) (pyomo-pounce/tests/test_covariance.py).

The noise variance comes from, in order of precedence: sigma_sq= (known measurement variance); the declared residuals (estimated as SSR / (n - n_params), with both numbers derived from the container); or the n_data= fallback for models without explicit residuals, whose SSR is the objective value at the solve — like estimate()’s baseline, writing into the model afterwards (a measurement, a warm start for the next horizon) does not move the answer. The solve warns if the declared residuals do not reproduce the objective value (weights or regularization terms would silently corrupt the estimate).

Groups. declare_residual(m.r_conc, group="conc") partitions residuals into noise groups by arbitrary user strings: containers sharing a group (or all ungrouped containers) pool into one estimated variance; distinct groups get their own (cov.sigma_sq becomes a dict), and the covariance switches to the heteroscedastic sandwich form, whose per-group pieces come from the same backsolves. When groups genuinely differ, weighting the objective itself (dividing each group’s residuals by its sigma) is the statistically efficient fix; the sandwich is the truthful report on the unweighted fit.

cov.eigen() returns ascending eigenvalues and matching eigenvectors. An eigenvalue much larger than the rest flags a poorly identified problem: its eigenvector is the parameter combination the data cannot pin down, and the corresponding cov.correlation entries approach +/-1. The returned signs follow the project-wide eigenvector sign convention — the largest-magnitude component of each eigenvector is positive, ties broken by the earliest position in cov.params — so the direction reproduces across machines instead of coming back as whatever LAPACK’s build chose. information().eigen() is the same. covariance warns when the held factor carries inertia-correction perturbations (typically an exactly unidentifiable parameterization) and when the covariance diagonal comes out negative (not a least-squares minimum).

Bound and constraint activity is classified from the solve’s own barrier geometry, not a slack threshold. A STRONGLY ACTIVE bound pins its parameter: zero variance, correlations 0, conditional on the bound, warned. A WEAKLY ACTIVE bound (slack and multiplier vanish together) is KEPT at its full finite variance, corrected for the barrier weight the held factor carries; AMBIGUOUS (loosely converged) and UNIDENTIFIED (curvature below the model’s own noise scale) stay in the free block, each with a warning. A strongly active inequality CONSTRAINT over the fitted parameters pins a combination rather than a coordinate: the matrix is projected on the constraint’s null space, going singular by one per binding row, and the warning names the constraint, the pinned combination, and its conditional information. The same limit written as a bound or as a row returns the same matrix. A binding row that reaches the fitted parameters through free eliminated variables cannot be represented by a restricted normal and is kept unprojected with an explicit warning.

To classify honestly, the declaration-triggered solve sets bound_relax_factor = 0 (slacks must measure distance to your own bounds). This applies to every solve routed through the sensitivity session, not only ones that end in covariance(). If you need the relaxation, pass bound_relax_factor explicitly in options=: your value wins, and covariance() then refuses with a clear error rather than classifying against shifted slacks.

Relation to pounce.curve_fit. This uses the same scale-and-invert-the-reduced-Hessian recipe as pounce.curve_fit — both read a reduced-Hessian block from the held KKT factor and scale it by 2 sigma^2 with sigma^2 = SSR / (n - p) — but with one substantive difference for nonlinear models: curve_fit factors the Gauss-Newton Hessian (pcov = 2 sigma^2 (J^T J)^-1, the expected-information / scipy / pycse.nlinfit convention, always positive semidefinite), while covariance() here feeds the exact Lagrangian Hessian through the .nl bridge, so it reports the observed-information covariance — the full reduced Hessian including the residual-curvature term that Gauss-Newton drops. The two are identical for linear models and in the small-residual / large-n limit, and differ by O(residual x model curvature) otherwise (a few percent on a strongly-curved fit). Neither is uniquely “correct”: Gauss-Newton is the conventional, robust default (it cannot produce a negative variance); observed information is the honest local curvature of the objective you actually solved (Efron & Hinkley 1978) but can go indefinite — which is what the negative-variance warning above is telling you. covariance() offers both: the default hessian="lagrangian" inverts the exact reduced Hessian of the Lagrangian, and covariance(m, hessian="gauss-newton") rebuilds the expected-information form from the residual Jacobian, recovered from the same backsolves at no extra solve (declared residuals required). Reach for it when the numbers must match scipy/nls, when covariance() warns about a negative diagonal, or when the covariance must stay positive semidefinite by construction, e.g. feeding an arrival-cost update in moving horizon estimation. The other difference is the input surface. curve_fit(f, xdata, ydata, ...) is the batteries-included fitter for a callable model f(x, *params) and data arrays: it chooses a starting point, offers robust losses, per-point sigma weights, confidence intervals, prediction bands, dpopt/ddata, and out-of-core streaming, and it projects the covariance onto the active-constraint nullspace when a parameter sits on a bound. covariance() is the post-solve primitive for a model you have already written in Pyomo — residuals as constraints, arbitrary surrounding structure — where you want the covariance of the fit as posed without re-expressing it as f(x, *params). Use curve_fit when the fit is naturally a model-plus-data call; use covariance() to interrogate an existing Pyomo estimation model. Both project a bound-active fitted parameter onto the active-constraint nullspace: covariance() reports the covariance conditional on the active bound (zero variance in the pinned direction, computed by inverting the free block of the information matrix) and still warns, since boundary asymptotics are nonstandard. Only variable bounds on the fitted parameters themselves are detected; a parameter held at the same value by an active constraint row is treated as free (#362). A bound rewritten into a constraint by the rule in Declared Params in variable bounds is the one exception: the value it held at the solve point is recorded, so it is still detected and still projected.

Relation to pyomo.contrib.parmest. parmest is an estimation workflow: multi-experiment data management, bootstrap resampling, and likelihood-ratio confidence regions, at the price of restructuring the problem into its experiment framework, with covariance computed by finite differences or an ipopt re-solve. covariance() is a post-solve primitive: the model as written, one declaration per component, the asymptotic covariance and identifiability diagnostics from the factorization the solve already produced. Use parmest for multi-experiment campaigns and non-asymptotic intervals; use this to interrogate the fit you already have.

See python/notebooks/26_parameter_covariance.ipynb for a worked example with a Monte Carlo validated confidence ellipse and an identifiability diagnosis.

Activity classification

Which bounds and constraint rows are actually holding the solution is a question the converged iterate answers only ambiguously: at a weakly active bound the slack and its multiplier are both O(√μ), so no fixed threshold on either one alone separates “just touching” from “not binding”. Solver.classify_activity() keys on the ratio of barrier curvature to the model’s own curvature instead, which is O(μ), O(1) and O(1/μ) in the three regimes:

solver = pounce.Solver(problem)   # problem.add_option("bound_relax_factor", 0.0)
x, info = solver.solve(x0=x0)

rep = solver.classify_activity()
rep["var_status"]        # ["inactive", "unbounded", "fixed", "strongly_active"]
rep["row_status"]        # ["equality", "strongly_active"]
rep["var_ratio"]         # the ratio behind each call (NaN where nothing was classified)
rep["mu"]                # the barrier parameter the calls were made at

Statuses are inactive, weakly_active, strongly_active, ambiguous (the ratio fell in a gap where this μ cannot decide — re-solve tighter), and unidentified (the curvature is below noise scale, so the question does not arise). unbounded, fixed and equality mark entries with no barrier geometry to classify.

Both arrays are indexed in user space: var_* follows your n variables and row_* your m constraints, in your order. A variable that fixed_variable_treatment = make_parameter removed from the solve (lb == ub) reports fixed at its own index rather than shifting everything after it.

Two per-entry flags report on the assumptions rather than the geometry: off_central_path (s·z differs from μ by more than 10× on some side) and contaminated (classified inactive yet carrying barrier curvature well above the O(μ) an inactive bound should have — typically a bound that sits close enough to the optimum to bend it).

Inequality rows classify through the same rule, via the curvature along the constraint normal. That is the point of classifying rows at all: move a bound off a variable and onto a row and the activity disappears from the bound-multiplier view entirely, while any covariance or identifiability heuristic keyed on z alone silently stops seeing it (#362).

The call requires the solve to have run with bound_relax_factor=0 (the Ipopt default is 1e-8) and raises ValueError otherwise: relaxed bounds shift the very slacks the classifier reads. The guard tests the value that solve ran under, so setting the option after the fact does not change the answer — set it on the Problem and solve again.

The information matrix

information(model) is the un-inverted sibling of covariance(): the reduced Hessian over the declared fitted block, from the same single solve, in natural units with no sigma^2 anywhere. For a homoscedastic Lagrangian fit, covariance() equals 2*sigma^2*inv(information()) on the free block. hessian= selects the observed ("lagrangian", default) or expected ("gauss-newton") form exactly as in covariance().

The Lagrangian form is built by tangent recovery against the held factorization rather than by inverting the covariance back: the K-inverse columns’ x-blocks are T*M, so T = Zx*inv(M) exactly and R = T'HT with the exact Lagrangian Hessian. The barrier weight cancels multiplicatively, so equality and variable-bound activity carries machine precision at any barrier parameter, including on pinned parameters where a subtract-the-barrier route loses log10(Sigma/q) digits. A binding inequality row is the one exception: it couples through its slack barrier and leaves ~1e-6 relative residue at practical barrier parameters.

Membership and warnings follow covariance(). One disposition is opposite by design: a strongly active (pinned) parameter’s entry is S, the reduction onto the pinned set, NOT a zero row — zero information is the opposite of what a pinned parameter carries — conditional on the rest of the pinned set, with zero cross blocks to the free parameters. Binding constraint rows project the free block on both sides (the pseudo-inverse of the projected covariance). An indefinite Lagrangian block is returned as computed with a warning naming Gauss-Newton as the PSD alternative: refusing would withhold the finding that the point is not a minimum or the model is over-parameterized. eigen() reads identifiability directly: a near-zero eigenvalue is a direction the data does not inform; its eigenvector’s sign follows the project-wide convention.

Choosing the block: wrt=

Both accessors take wrt= to reduce onto any block of the solve’s variables off the held factor, post-solve; the declared fitted block is the default, so omitting it is exactly the prior behavior. Accepted forms: a Var (scalar or indexed, every member), an indexed slice (m.x[2, :]), a (Var, iterable) pair, data objects, or a list mixing these.

cov = covariance(m)                      # the fitted block, as before
cov_a = covariance(m, wrt=[m.a])         # one parameter's marginal
band = covariance(m, wrt=m.r)            # a predicted trajectory
info_a = information(m, wrt=[m.a])

Each call re-reduces onto its own argument, so one solve serves as many blocks as are asked about, and each block gets its MARGINAL: everything outside it is profiled out, not held fixed. Sigma estimation always divides by the fit’s own degrees of freedom (a property of the solve, not of the question being asked), so a sub-block’s numbers agree exactly with the corresponding entries of the default answer.

A rank-deficient block, one with more coordinates than the fit has degrees of freedom or with linearly dependent coordinates (a duplicated design point), is the trajectory-band case: covariance() returns its (rank-deficient) marginal, 2 sigma^2 M, the confidence band on the fitted trajectory (add the observation noise for a prediction band), with the membership handling bypassed, and information() raises an error pointing to covariance(), since such a block carries no information matrix. For information(), a block that parameterizes the constraint manifold (size equal to the degrees of freedom) gets the exact tangent construction; a sub-block of the fitted set gets its marginal as a Schur complement of the exact tangent R over the fitted block (never inverting a covariance, so a pinned member costs no digits); other blocks reduce off the held factor with the item-1 corrections, which is benign for free coordinates.

One exception is returned rather than hidden: a strongly active variable OUTSIDE the block is not deleted from the factor, so the block’s numbers are the values conditional on that bound, not the marginal over it. The result carries the list as .conditioned_on (empty when there is none); inside-block activity is membership, not conditioning, and is handled as before. The list is decided by the same classification the block members get, applied per candidate as a singleton block, so it is scale-invariant; only near-bound variables pay the extra backsolve.

Keeping and releasing the factor: retain_kkt(), release_kkt()

The solve factors the KKT matrix to solve the NLP; the only question is whether that factor is kept for post-solve queries. Any declaration keeps it. retain_kkt(model) keeps it with no declaration at all, which is what wrt= queries with nothing declared need: the MHE case, where the arrival state and the parameters are each queried by wrt= and neither is THE fitted set. It defaults off, so a solve with no sensitivity pays nothing.

retain_kkt(m)
SolverFactory("pounce").solve(m)
arrival = covariance(m, sigma_sq=s2, wrt=m.x[:, t0])
params = information(m, wrt=[m.k1, m.k2])
release_kkt(m)          # done asking: give the memory back now

The retention policy in one place: the factor is kept if anything is declared or retain_kkt() was called, and a Covariance or Information result whose lazy conditioned_on has not been read keeps the session alive through its pending computation until first access. release_kkt(model) is the exit: it drops the model’s hold on the factor immediately, freeing the memory, while declarations and the retain flag still apply to the next solve. Release drops the model’s hold, not a result’s: a Covariance or Information with a pending conditioned_on, and a Gradient (which reads the factor on every lookup), each hold their own reference, so they keep working across the release and keep the factor in memory until they are discarded. Noise is a separate question: retain_kkt() keeps the factor, not a noise model, and with nothing declared fitted the degrees of freedom for a noise ESTIMATE are unknown, so covariance() under retain-only needs sigma_sq=; the estimation routes (declared residuals, n_data=) raise an error saying so.

Like any declaration, retain_kkt() routes the solve through the in-process sensitivity path, whose solve() surface is not keyword-identical to the ordinary subprocess path (for example, load_solutions=False is not honored there). Adding it to an existing script changes how the solve runs, not just what is kept.

Units and NLP scaling

All sensitivity outputs are in natural (unscaled) units. The IPM holds its converged KKT factor in an internally scaled space whenever NLP scaling is active (the default nlp_scaling_method = "gradient-based" fires when an objective gradient or constraint row exceeds nlp_scaling_max_gradient = 100 at the starting point); pounce undoes that scaling in every held-factor back-solve, so dx, kkt_solve, and the reduced Hessian are independent of how the problem was scaled internally (#128).

That covers user scaling too, on all three of its axes. A per-variable scaling_factor is applied as a change of variables x̃ = d ⊙ x below the algorithm, so the held factor is the scaled problem’s; the factors are carried into the same translation, and every accessor answers in your units (#486). The factors a solve ran under are readable back from Solver.nlp_scaling["x_scaling"] (Python) / Solver::variable_scaling (Rust) — diagnostic rather than a correction to apply, since the outputs already carry it.

classify_activity() is scale-invariant for the same reason, and mostly by construction rather than by undoing anything: its ratios are formed so that rescaling a constraint row or the objective leaves them fixed. Writing a constraint as 1000·x ≥ 0 instead of x ≥ 0 does not move a status, and neither does the solver’s own per-row d_scale. A change of variables is the one case the ratios do not absorb on their own — the identification floor is a single number shared across entries, so a non-uniform d would move entries across it — and there the factors are divided out of the geometry before anything is classified, which keeps a status from depending on the conditioning you asked for. The values the report exports follow the natural-units contract like everything else: var_sigma and row_sigma are the barrier diagonals in the model’s own units, row_normal(j) is the constraint gradient with the solver’s per-row scale divided out, and hessian_vec(v) is the exact Lagrangian Hessian times a user-space vector with the objective scale divided out; classification happens on the scaled quantities internally, the report never shows them.

Variable indices are user-space, factor rows are not. Everything the sensitivity API reports or accepts — the .col file’s order, the activity report’s var_* arrays, row_normal(j)’s entries — indexes the variables you wrote. The converged factor does not: a variable whose bounds are equal is removed from the solve (fixed_variable_treatment = make_parameter, the default), so its column is absent and every later variable sits one row earlier. The two orders coincide exactly when the model has no fixed variable, which makes the difference easy to miss. Translate with Solver.primal_rows(indices) — None marks a removed variable — before indexing a kkt_solve or parametric_step_full result, just as multiplier_rows has always been required for the y_c block.

In particular, for a parameter-estimation NLP with the parameters pinned by equality constraints, -inv(info["reduced_hessian"]) is directly the parameter covariance — no per-problem scale factor, no need to set nlp_scaling_method = "none". (Sign convention: over pin constraint rows, B K⁻¹ Bᵀ equals the multiplier sensitivity ∂λ/∂p = −∂²f*/∂p², hence the minus in the covariance recipe.)

For callers that calibrated against the pre-#128 behavior, the solver-space value and the factors that relate the two are exposed:

• Python: info["reduced_hessian_scaled"], info["obj_scaling_factor"], info["pin_g_scaling"]; Solver.reduced_hessian(pins, scaled=True), Solver.kkt_solve(rhs, scaled=True), and the Solver.nlp_scaling dict ({"obj": df, "c_scale": …, "d_scale": …, "x_scaling": …}).
• Rust: SensResult::{reduced_hessian_scaled, obj_scaling_factor, pin_g_scaling}, Solver::{compute_reduced_hessian_scaled, kkt_solve_scaled, nlp_scaling, pin_g_scaling}, and PdSensBacksolver::solve_scaled_space.

The relation is H_scaled[i,j] = df / (dc_i·dc_j) · H[i,j], where df is the objective scaling factor and dc_i the pin rows’ constraint scaling factors.

One caveat: the IPM’s inertia-correction perturbations (δ_x, δ_s, δ_c, δ_d) are added to the factor in scaled space, so on a problem whose final factorization needed regularization (e.g. linearly dependent pin rows) the unscaling maps a slightly different perturbed system per scaling method. The perturbations are reported — info["kkt_perturbations"] / Solver.kkt_perturbations (Python), SensResult::kkt_perturbations / Solver::kkt_perturbations (Rust) — so a covariance workflow can assert they are all zero before trusting -inv(reduced_hessian); on well-posed estimation problems the final factor is unregularized and the invariance is exact.

Verification

All three entry points are verified against upstream sIPOPT 3.14.19’s parametric_cpp golden output to within roughly 6e-9 per component. The bound projection is a single-pass clamp; upstream’s iterative Schur refinement (re-factorize on each violation) is intentionally not ported.

Beyond one perturbation

Everything above answers “how does x* move for this \(\Delta\theta\)” — a first-order step off one converged factor. Repeat it and you are tracing a path, at which point the questions become where the linear prediction stops being good enough, when the active set changes under you, and what to do where \(\partial x^*/\partial\theta\) goes singular.

The Python frontend answers those with PathFollower, which turns the same held factor into a predictor–corrector continuation loop (and a pseudo-arclength mode that traces through folds), plus inverse_map_rhs for running the map backwards as an ODE. See Path Following & Inverse Mapping.

Sessions: Factor-Once / Solve-Many

POUNCE’s IPM converges to a KKT linear system that, once factored, answers a number of useful follow-up questions cheaply: parametric steps, reduced Hessians, custom back-solves. The session APIs let you hold that factor alive between operations, rather than rebuilding it on every call. The same machinery serves two workloads:

• Sensitivity / many-RHS. After one solve, issue many cheap operations against the converged factor — parametric steps for several parameter perturbations, reduced Hessians over several pinned-row sets, raw KKT back-solves.
• Factor-only. For non-IPM uses (shift-invert eigensolves, custom Newton iterations) the underlying [Factorization] handle in pounce-linsol exposes factor / refactor / back-solve directly, without the IPM in the loop.

Which layer do I want?

The session API does not rebuild the IPM. Each solve() call runs the full barrier method from scratch. What it reuses is the factor that exists at convergence: KKT back-solves and sensitivity operations skip the symbolic factor, AMD ordering, and numeric factorization.

Python

import pounce

problem = pounce.Problem(...)
solver = pounce.Solver(problem)

x, info = solver.solve(x0=x0)
assert solver.converged

# Parametric step ∂x*/∂p · Δp, with p pinned by g(x) row indices.
dx = solver.parametric_step([2, 3], [-0.5, 0.0])

# Reduced Hessian B K⁻¹ Bᵀ over the same pinned-row set.
hr = solver.reduced_hessian([2, 3])

# Raw KKT back-solve, useful for custom workflows.
dim = solver.kkt_dim
rhs = np.zeros(dim)
lhs = solver.kkt_solve(rhs)

# Which bounds and rows actually hold the solution, in user index order.
rep = solver.classify_activity()          # needs bound_relax_factor=0
rep["var_status"], rep["row_status"]

# Constraint-row gradient in user space, natural units.
a = solver.row_normal(j)

# Exact Lagrangian Hessian times a user-space vector, natural units.
hv = solver.hessian_vec(v)

The KKT compound vector is laid out as x || s || y_c || y_d || z_l || z_u || v_l || v_u. pin indices in parametric_step / reduced_hessian are 0-based row indices into g(x); they are mapped internally to the matching y_c rows (through the equality/inequality split, so inequalities may precede the pins).

All back-solves are in natural (unscaled) units — any NLP scaling the IPM applied internally is undone, so results are independent of nlp_scaling_method (#128). The solver-space values remain available via reduced_hessian(pins, scaled=True) / kkt_solve(rhs, scaled=True), and the factors via the Solver.nlp_scaling dict — see Sensitivity Analysis.

pounce.Problem.solve() and Problem.solve_with_sens() still work unchanged — each internally builds a fresh session — but new code that issues more than one sensitivity query per solve should prefer pounce.Solver to skip rebuilding the application.

C

IpoptProblem prob = CreateIpoptProblem(...);
AddIpoptStrOption(prob, "linear_solver", "feral");

/* Consumes prob — the IpoptSolver is now the sole owner.
   prob is NULLed; calling FreeIpoptProblem(prob) on the now-null
   pointer is harmless. */
IpoptSolver sol = IpoptCreateSolver(&prob);

double x[n], obj;
IpoptSolverSolve(sol, x, NULL, &obj, NULL, NULL, NULL, user_data);

Index dim = IpoptSolverGetKktDim(sol);     /* compound KKT dim     */
double rhs[dim], lhs[dim];                  /* memset rhs as needed */
IpoptSolverKktSolve(sol, rhs, lhs);

Index pins[2] = {2, 3};
double deltas[2] = {-0.5, 0.0};
double dx[n];
IpoptSolverParametricStep(sol, 2, pins, deltas, dx);

double hr[2 * 2];                           /* column-major dense   */
IpoptSolverReducedHessian(sol, 2, pins, 1.0, hr);

IpoptFreeSolver(sol);

The classic IpoptSolve API is unchanged and unaffected; the session handle lives alongside it.

Rust

Both session APIs come through the pounce-rs facade. Solver needs the sensitivity feature; the bare Factorization below needs convex or qp, whichever you are already using — either one enables pounce_rs::linsol.

#![allow(unused)]
fn main() {
use pounce_rs::sensitivity::Solver;

let mut solver = Solver::new(app, tnlp);
solver.solve();
assert!(solver.converged().is_some());

let dx = solver.parametric_step(&[2, 3], &[-0.5, 0.0])?;
let hr = solver.compute_reduced_hessian(&[2, 3], 1.0)?;

let mut lhs = vec![0.0; solver.kkt_dim().unwrap()];
solver.kkt_solve(&rhs, &mut lhs)?;
}

For purely linear-algebra uses with no IPM in the loop:

#![allow(unused)]
fn main() {
use pounce_rs::linsol::{Factorization, backend};

let mut fact = Factorization::new(dim, ia, ja, values, backend())?;
fact.solve(&mut rhs, 1)?;          // back-substitute in place
fact.refactor(&new_values)?;       // pattern preserved; numeric reuse
fact.solve_one(&mut another_rhs)?;
}

What’s preserved across operations

• Symbolic factor / AMD ordering. Owned by the linear-solver backend; reused on every back-solve and on refactor().
• Numeric factor. Reused on every back-solve until you refactor.
• The converged primal-dual state (x*, multipliers, g(x*), iteration stats).

What’s not preserved across solve() calls

The session is currently a factor-and-query value: one solve, many follow-up operations. A separate resolve() that re-runs the IPM while reusing the symbolic factor + AMD ordering across top-level solves (for MPC / B&B / warm-start workloads) is planned but not yet implemented. Each solve() call today runs a fresh IPM.

Verification

All session entry points are tested for numerical equivalence with the corresponding one-shot APIs:

• pounce.Solver.solve ≡ Problem.solve (1e-12).
• pounce.Solver.parametric_step ≡ Problem.solve_with_sens(deltas=…)['dx'] (1e-10).
• pounce.Solver.reduced_hessian ≡ Problem.solve_with_sens(compute_reduced_hessian=True)['reduced_hessian'] (1e-10).
• pounce_rs::sensitivity::Solver::parametric_step ≡ SensSolve::with_deltas (1e-10).

See python/tests/test_solver_session.py and crates/pounce-sensitivity/tests/solver_session.rs for the full test matrix.

Differentiable Solves & the DiffHandoff Contract

POUNCE solves are differentiable: a solve can sit inside a JAX or PyTorch model and pass gradients with respect to the problem parameters. This page documents the handoff contract — the well-defined bundle of post-convergence data every solve produces — so that any consumer (the built-in JAX/Torch layers, a downstream tool such as discopt, or your own autodiff code) can differentiate a POUNCE solve from one stable surface rather than from solver internals.

Design notes: dev-notes/diff-handoff-contract.md.

What a differentiable backward needs

The gradient of an optimal solution x*(p) with respect to a parameter p comes from the implicit-function theorem applied to the KKT conditions at the solution. To assemble it, a backward pass needs:

1. the primal solution x* and the constraint / bound multipliers;
2. the active set — which variable bounds bind and which constraint rows are active — so inactive directions drop out correctly;
3. (for performance) the converged KKT factorization, reused as a back-solve rather than rebuilt.

POUNCE produces all three. The first two ride out in the solve info dict; the third is reused automatically by JaxProblem (see Sessions).

The active-set masks (the DiffHandoff core)

Every NLP solve’s info dict carries a precomputed active set, derived once on the Rust side (pounce_sensitivity::DiffHandoff) so no consumer re-derives it under its own tolerance:

pinned_vars is the seam used for mixed-integer problems: a branch-and-bound leaf fixes integer variables at their optimal values, and those variables differentiate exactly like an active bound (dx/dp = 0). A producer of a fixed-integer leaf adds them to the mask (DiffHandoff::pin on the Rust side).

Multiplier conventions (canonical mapping)

The same dual quantity is named differently across POUNCE’s solver surfaces — deliberately, because each surface preserves an external contract. The canonical field is DiffHandoff.lambda (general constraint multipliers); this table maps every surface onto it so a consumer knows the correspondence:

Caution. The internal symbol lam is not a single quantity: in the NLP backward (jax/_diff.py) it is all constraint multipliers (= mult_g); in the QP backward (jax/_qp.py) it is the inequality-only duals (= z). Always map through the table above rather than assuming a shared name means a shared quantity.

These names are stable: the NLP keys are an external cyipopt contract and will not be renamed.

Consuming the contract

JAX / PyTorch (built in)

pounce.jax and pounce.torch already differentiate solves; you do not touch the masks directly. Use pounce.jax.solve / JaxProblem (or the torch equivalents) and call jax.grad / .backward() as usual. For batched and repeated solves, JaxProblem reuses the converged KKT factor in the backward (factor_reuse=True, default) — see Sessions.

Across a language / tool boundary (e.g. discopt)

A downstream tool that drives POUNCE as its NLP backend and composes its own autodiff reads the contract straight from the info dict returned by Problem.solve:

x, info = problem.solve(x0=...)
# primal + duals
lam   = info["mult_g"]      # general-constraint multipliers (the canonical λ)
z_L   = info["mult_x_L"]
z_U   = info["mult_x_U"]
# precomputed active set — do NOT re-derive |mult| > tol yourself
pinned = info["pinned_vars"]          # bool[n]: dx/dp = 0 on these
active = info["active_constraints"]   # bool[m]: rows in the KKT block
tol    = info["active_tol"]

Because the active set is computed once in the producer, every consumer sees the same masks under the same tolerance — which is what makes a gradient assembled on one side of the boundary agree with one assembled on the other.

Verification

The contract is exercised by the test suite:

• python/tests/test_problem.py::test_diff_handoff_masks_in_info asserts the masks against a problem with a known active set (HS071: one variable on its lower bound, a binding inequality, and an equality).
• python/tests/test_jax.py (85 finite-difference gradient checks) and python/tests/test_parity_jax_torch.py (JAX↔Torch gradient agreement) confirm the backward passes that rest on this data are correct and frontend-independent.

Interactive Solver Debugger

POUNCE ships an interactive debugger for the interior-point loop — a pdb for the IPM. You can pause the solve at well-defined points, inspect and mutate the live mathematical state (the iterate, multipliers, the barrier parameter μ), set breakpoints (by iteration, on a numeric condition, or on a solver event), step through an iteration’s internal phases, rewind to an earlier iterate, re-solve from a saved point with new options, and drop in automatically when a solve fails.

It has two front ends sharing one command engine:

• a human REPL (--debug) with history, Ctrl-R search, and Tab completion, and
• a newline-delimited JSON protocol (--debug-json) that an LLM agent, a script, or a visual debugger (e.g. a VS Code Debug Adapter) can drive programmatically.

No production NLP solver ships anything like this; if you have used ipopt you have had print_level and a log. This is a live debugger.

The same debugger spans every POUNCE solver — the NLP filter-IPM and the convex / conic interior-point solver share one command engine and one REPL. See Beyond the interior-point loop for the small set of commands whose availability is backend-conditional.

The debugger has zero effect on the solve when it is not attached. The checkpoint fire-sites short-circuit when no debugger is installed, so the standard regression suite is bit-for-bit identical with and without the feature compiled in.

Quick start

pounce problem.nl --debug          # human REPL, pauses at iteration 0
pounce problem.nl --debug-json     # JSON protocol on stdin/stdout
pounce problem.nl --debug-on-error # run freely; drop in only if it fails
pounce problem.nl --debug-on-interrupt   # run; Ctrl-C drops you in

A 30-second session (human REPL):

$ pounce --problem rosenbrock --debug

── pounce-dbg ── iter 0 @iter_start  mu=1.000e-1  obj=2.420000e1  inf_pr=0.00e0  inf_du=1.00e2
pounce-dbg> info
iter      = 0
mu        = 1.000000e-1
objective = 2.42000000e1
...
pounce-dbg> print x
x = [-1.200000e0, 1.000000e0]
pounce-dbg> break if inf_du<1e-6
conditional breakpoint: inf_du<1e-6
pounce-dbg> continue
... solver runs ...
── pounce-dbg ── iter 21 @iter_start  mu=...  inf_du=8.7e-7
   ↳ inf_du<1e-6
pounce-dbg> quit

The prompt is on stderr; the solver’s own iteration table stays on stdout, so a redirected log is unaffected.

The two front ends

On a non-TTY stdin (a pipe), the REPL falls back to a plain line reader (no history/Tab) but otherwise behaves identically — handy for scripted tests.

The JSON protocol is documented in full below.

Pausing and flow control

Checkpoints

The loop fires the debugger at these points (a pause reports which one via its checkpoint field):

By default the debugger only stops at iter_start (and terminated). The sub-iteration checkpoints fire every iteration but resume immediately unless you ask to stop at them.

Stepping into restoration. The same debugger drives the restoration inner IPM: when the solve enters restoration, the inner solve’s checkpoints fire too. A step/stepi that lands on an inner iteration pauses there with in_restoration: true (REPL banner shows [restoration]), and print x shows the restoration sub-NLP iterate. stop-at resto (pre_restoration_entry) is the easy way to catch the hand-off and then step inward.

Stepping

stop-at takes a checkpoint name or a friendly alias:

stop-at after_search_dir     # or:  stop-at kkt
stop-at pre_restoration_entry   # or:  stop-at resto
stop-at                      # list active stop-at checkpoints
stop-at clear

Aliases: mu → after_mu, kkt/search_dir → after_search_dir, step → after_step, resto → pre_restoration_entry, resto_exit → post_restoration_exit.

Breakpoints

Three kinds, all reported in break and surfaced as the pause reason.

By iteration

break 12            # pause at iteration 12   (alias: b 12)
tbreak 12           # one-shot: pause at 12, then delete itself (alias: tb)
break               # list all breakpoints
break del 12        # remove
break clear         # remove everything (iters + conditions + events)

Watchpoints (data breakpoints)

watchpoint x[3]        # pause when x[3] changes (alias: wp)
watchpoint x 1e-3      # pause when any x component moves by > 1e-3
watchpoint             # list; watchpoint del x[3]; watchpoint clear

Distinct from watch (which only displays): a watchpoint pauses the solve when the watched value changes by more than its threshold (default 0 = any change) between iterations. Useful for a component expected to stay put (e.g. a variable pinned at a bound).

Breakpoint command lists

Attach commands to a breakpoint that run automatically when it hits — semicolon-separated, ending with a flow command to auto-resume:

break 5
commands 5 print kkt ; set mu 0.1 ; continue   # at iter 5: inspect, tweak μ, go
commands 5 clear                               # remove
commands                                       # list all

When iteration 5 is reached, the debugger emits the pause, runs the attached commands (each result is reported), and if one of them resumes/stops, honors it without dropping to the prompt — otherwise it falls through to the interactive prompt as usual.

Conditional (with compound predicates)

break if inf_pr<1e-6
break if mu<1e-4 && inf_pr>1e-3
break if iter>10 && (inf_du>1e-2 || obj<0)
break clear cond

• Metrics: mu, inf_pr, inf_du, obj, err (overall NLP error), iter.
• Operators: <, <=, >, >=, == (== is float-tolerant).
• Compound: && and ||, evaluated strictly left-to-right with no precedence; parentheses are accepted but stripped (they don’t group). For real grouping, register several conditions — any one that holds fires.

Conditions are evaluated at iter_start.

On a solver event

break on regularized
break on resto_entered
break clear events

Events fire at whatever checkpoint makes them observable (e.g. regularized at after_search_dir), and pause with reason: "event: <name>".

Inspecting state

info                 # one-line summary: iter, mu, obj, inf_pr, inf_du, nlp_error, dims
print x              # a primal/dual block (alias: p x)
print dx             # a search-direction block (d + block name)
print mu             # a scalar: mu|obj|inf_pr|inf_du|err|compl|iter
print kkt            # KKT inertia + regularization (see below)
print rank           # SVD numerical rank of the equality Jacobian J_c (see below)
print active         # which bound categories are near-active (small slack)
watch mu             # auto-print a target at every pause (alias: display)
watch                # list watches; watch del mu; watch clear

watch <target> registers any print target (block, dx, scalar, kkt) to be shown automatically at every subsequent pause — the debugger’s equivalent of gdb’s display. In JSON mode the values arrive in the pause event’s watches array.

Blocks (the eight components of the primal-dual iterate):

Prefix any block with d (dx, dz_l, …) to print the corresponding block of the most recent Newton step.

Model names (.col / .row)

A solver-internal diagnostic that says “variable 132 in equation 3 looks singular” is far less actionable than one that says “T_reactor in energy_balance”. Lee et al. (2024) identify this gap — between detecting an issue numerically and tracing it back to a named equation in the modeling environment — as a central roadblock for debugging equation-oriented models.¹

AMPL .nl files carry no names, but AMPL emits two optional sibling files when the modeler sets option auxfiles rc;:

When these sit next to the .nl, pounce captures them (NlProblem::var_names / con_names) and exposes them through the ExpressionProvider::variable_name / constraint_name seam. Missing or malformed name files are non-fatal — names are a diagnostic aid, never load-blocking, so the debugger simply falls back to index labels.

print residuals uses these names directly. Residual values live in the solver’s split space (equalities and inequalities separated, fixed variables removed), so a name only labels the right row if it is carried through the same permutations. The TNLP publishes its .col/.row names under the conventional idx_names metadata key, and OrigIpoptNlp projects them into split space (x_not_fixed_map for variables, c_map for equalities, d_map for inequalities) — the debugger reads the result via DebugCtx::split_names. So a near-singular equality residual prints as

c[energy_balance] = +3.142e-04   |3.142e-04|

instead of c[3]. The same idx_names pool labels grad_x_L[...] (variable names) and grad_s_L[...] / d-s[...] (inequality names). The JSON payload keeps the numeric index and adds a name field.

Status. Capture, exposure, and print residuals labeling are live on the AMPL .nl path with names projected through the bound / c-d-split permutations. Presolve renumbers rows, so PresolveTnlp declines idx_names rather than risk mislabeling a permuted row — under presolve the debugger safely falls back to index labels. Carrying names through the presolve map and decorating print active are the next steps built on this foundation.

print equation — the algebra of a named constraint

Naming a culprit row is only half the story; the next question is always what does that equation actually say? Lee et al. (2024) make this the core of actionable equation-oriented diagnostics — a debugger should surface the named equation, not just a row index.¹ print equation closes that loop: once print residuals points at, say, c[energy_balance], you read the constraint’s source algebra directly.

(dbg) print equation energy_balance
energy_balance:  T_reactor*flow - 300*flow - Q = 0

(dbg) print equation 14          # by original .nl row index
c[14]:  x[3]^2 + x[7]^2 <= 1

A constraint is addressable by its model name (preferred, and robust to row reordering) or its original .nl row index. With no argument, print equation reports how many equations are available. The renderer works from the faithful Expr DAG the .nl parser built — not the lossy evaluation tape — so common-subexpressions, imported functions, and piecewise/conditional forms render as written. The affine part is printed with tidy signs (a - 2*b, not a + -2*b), zero-coefficient Jacobian placeholders are suppressed, and bounds render in their natural relation (= rhs, lo <= body <= hi, >= lo, <= hi). The JSON payload carries {index, name, equation}.

Equations are static model data in original .nl row order, so unlike residuals they need no split-space projection — print equation works regardless of presolve. It is available whenever a model was loaded from an .nl file; the JSON name field is present only when a .row auxfile supplied one.

print kkt — inertia and regularization

Available at/after after_search_dir (use stop-at kkt). This is the view a solver expert reaches for when a step looks wrong:

pounce-dbg> stop-at kkt
pounce-dbg> continue
── pounce-dbg ── iter 3 @after_search_dir ...
pounce-dbg> print kkt
dim       = 3
inertia   = n+=2 n-=1 (expected n-=1) → correct
delta_w   = 0.000000e0   (primal regularization)
delta_c   = 0.000000e0   (dual regularization)
status    = Success

The augmented (KKT) system has expected inertia (n₊ = n, n₋ = m, n₀ = 0) where m is the number of equality + inequality multipliers. A mismatch — or a nonzero delta_w/delta_c — is the classic signal that the step is being stabilized (the solver added regularization to fix the inertia).

For the matrix and factor themselves:

viz kkt     # the assembled augmented-system matrix (triplets) + inertia
viz L       # the LDLᵀ factor (strict-lower triplets + values)

viz kkt writes the KKT matrix as 1-based lower-triangle triplets (dim, irn, jcn, vals) alongside the inertia summary — point $POUNCE_DBG_VIEWER at a heatmap script. viz L writes the LDLᵀ factor (n, fill-reducing perm, strict-lower l_irn/l_jcn/l_vals in permuted coordinates), read out of the factor the solver actually computed.

Both are read-only and always show the most recent factorization: the current iteration’s system at an after_search_dir stop, or the previous iteration’s at the default iter_start pause (the step that produced where you’re standing). The matrix and factor are captured every iteration while the debugger is stepping; once you detach (run free) the capture is dropped — so on a large problem a free run doesn’t pay the O(nnz) assembly. If you viz kkt/viz L right after a free run, step once to re-capture.

print rank — numerical rank of the equality Jacobian

print kkt tells you that the dual system needed regularization (delta_c > 0) or that the inertia was wrong; the structural_singularity finding names equations that are dependent by sparsity pattern. print rank closes the last gap: a rank-revealing SVD of the equality Jacobian J_c at the current iterate. It factors the matrix the solver actually sees (constraint scaling already applied), so it localizes the dependency to specific equations — including dependencies that are numerical only (values that cancel over a full sparsity pattern), which the structural Dulmage–Mendelsohn pass cannot detect.

It doesn’t just name the culprit equations — it prints them. When a .nl model is loaded, each implicated row’s source algebra is rendered directly beneath it (the same DAG-faithful text print equation shows), so you read the dependency without a second command:

pounce-dbg> print rank
equality Jacobian J_c: 3 row(s) × 4 column(s)
numerical rank = 2 / 3  (deficiency 1)
σ_max = 3.162e0   σ_min = 0.000e0   cond = inf (σ_min = 0)   (rank tol τ = 1.40e-15)
singular values: [3.162e0, 1.414e0, 0.000e0]
rank-deficient: 1 equation(s) lie in the near-null space (linearly dependent / redundant) — the source of δ_c regularization:
  c[mass_balance]       (participation 0.50)
      x[0] + x[1] - 10 = 0
  c[mass_balance_dup]   (participation 0.50)
      x[0] + x[1] - 10 = 0

The two equations print identically — that is the redundancy, now visible on its face.

For the SVD J_c = U Σ Vᵀ, the left singular vectors u_k whose singular value σ_k ≈ 0 span the left null space — the row combinations u_kᵀ J_c ≈ 0 that vanish. Each row’s participation w_i = Σ_{k : σ_k ≤ τ} u[i,k]² ∈ [0, 1] localizes the dependency: a redundancy shared between two equations splits ≈ 0.5/0.5, while w_i = 1 means row i lies entirely in the null space. The numerical-rank threshold is the standard LAPACK/NumPy τ = σ_max · max(m, n) · ε; the implicated rows are resolved to model names through the same .row plumbing as print residuals / print equation.

The inline algebra is resolved by model name, so it appears for named rows. The rank report’s row index is the split equality position, not the original .nl row the equation source keys on, so an unnamed row can’t be mapped — there print rank falls back to a print equation <name> hint instead of guessing. When J_c has full row rank, that is reported as a positive signal (J_c has full row rank at this iterate.) with the σ_min/cond witnessing how far it is from degenerate — silence would be ambiguous. The command is available whenever the iterate has an equality block; a problem with no equality constraints returns a short explanatory error. The JSON payload is {iter, n_rows, n_cols, rank, deficiency, rank_deficient, sigma_max, sigma_min, cond, tol, singular_values, culprits: [{row, kind, index, name, label, weight, equation}]} (equation is the rendered source or null when unresolved; cond is null when σ_min = 0, since JSON has no infinity).

diagnose — a live, named health report

info, print residuals, and print kkt each expose one facet of the current iterate. diagnose (alias diag) runs a panel of heuristics over all of them at once and returns a ranked list of findings — and, crucially, names the culprit equation or variable behind each numerical symptom. That last step is the actionable-diagnostics path of Lee et al. (2024):¹ a report that says “mass_balance is the worst constraint residual” is worth far more than “row 13 is infeasible.”

pounce-dbg> diagnose
[  error] primal_infeasible: Primal infeasibility 1.70e+02; worst constraint
         residual is c[mass_balance] = +1.701e+02. Inspect this equation's
         feasibility and scaling (`print equation mass_balance`).
[warning] dual_infeasible: Dual infeasibility 9.84e-01; largest stationarity
         residual is grad_x_L[T_reactor] = -9.838e-01.
[warning] inertia_wrong: KKT inertia is wrong (n-=2 vs expected 1): the system
         was indefinite/singular and the step had to be stabilized.
[   info] bounds_pinned: 3 variable bound(s) are active (slack < 1e-6).

This is the live counterpart to the pounce-studio diagnose tool, which runs temporal heuristics over a finished solve report. The two share a {severity, code, message} shape so a client can treat them uniformly, but the live command sees what a saved report cannot: the current KKT inertia and regularization, and the named primal/dual residuals at this exact point. Findings are sorted error → warning → info; a clean iterate yields a single healthy finding. The checks:

KKT-derived findings (inertia_wrong, *_regularization) need a computed search direction, so they appear at/after after_search_dir. Names follow the same rule as print residuals: present on the .nl path with .col/.row files, index labels (c[13]) under presolve. The JSON payload is {iter, findings: [{severity, code, message}], n_findings}.

Structural rank: naming the dependent equations

inertia_wrong and dual_regularization detect a rank-deficient Jacobian, but only as a scalar — they tell you a redundancy exists, not which equations are redundant. structural_singularity closes that gap with a Dulmage–Mendelsohn decomposition of the equality Jacobian’s sparsity pattern (the same structural check at the heart of IDAES’s DiagnosticsToolbox). A maximum bipartite matching between equality rows and variables partitions the system; any over-determined block — more equations than the variables they jointly touch — forces at least one of those equations to be redundant or mutually inconsistent (LICQ fails). The finding lists those equations by model name, e.g.:

pounce-dbg> diagnose
[warning] structural_singularity: Constraint Jacobian is structurally singular
         (Dulmage–Mendelsohn): 2 equation(s) over-determine the 1 variable(s)
         they jointly touch (flow_rate), so ≥1 of them must be redundant or
         mutually inconsistent (LICQ fails on this block). Candidate dependent
         equations: mass_balance, mass_balance_dup. Inspect them with
         `print equation <name>`; this names the rows behind any δ_c
         dual-regularization / wrong-inertia signal.

This is the named-culprit payoff of Lee et al. (2024):¹ reporting “mass_balance and mass_balance_dup are linearly dependent” rather than “the Jacobian is singular.” The check is iterate-independent (it reads only the sparsity pattern), so unlike the KKT-derived findings it fires from iteration 0 — it can flag a structurally broken model before the solver ever stalls on it. It is suppressed for well-posed problems: an NLP with more variables than equality constraints is the normal case (the spare degrees of freedom are pinned by the objective, bounds, and inequalities), so only the over-determined side is reported, never the under-determined one. Available on the .nl path; names fall back to c[i]/x[i] when no .col/.row auxiliary files were emitted.

Numerical rank: the value-dependency the structure can’t see

structural_singularity reads only the sparsity pattern, so it is blind to a redundancy that lives in the values — three equations whose every entry is nonzero (a structurally full-rank pattern) but whose rows satisfy row₂ = row₀ + row₁ numerically. rank_deficient_jacobian is the numerical complement: it runs the same SVD as print rank over J_c at the current iterate and, when the numerical rank falls short, names the equations in the near-null space:

pounce-dbg> diagnose
[warning] rank_deficient_jacobian: Equality Jacobian J_c is numerically
         rank-deficient at this iterate: rank 2/3 (deficiency 1),
         σ_min=0.00e0, cond=inf (σ_min = 0). Linearly dependent or redundant
         equality constraints — the root cause behind δ_c regularization /
         wrong inertia. Implicated equations: c[mass_balance],
         c[mass_balance_dup].

Unlike the structural check, this one is iterate-dependent — it factors J_c at the current x, so it reflects the matrix the solver is actually regularizing and catches dependencies that only appear at certain points. The two checks are deliberately layered: structural_singularity fires from iteration 0 on the pattern alone; rank_deficient_jacobian confirms it numerically and, more importantly, surfaces the value-only dependencies the structural pass provably cannot. See print rank for the SVD math and the per-equation participation weights.

Mutating state

Mutations feed straight back into the solve.

set mu 0.5           # overwrite the barrier parameter
set x[2] 1.5         # overwrite one component of a block
set x 1.0,2.0,3.0    # overwrite a whole block (comma-separated)

Setting any block works (set z_l[0] 1e-3, …). Iterate edits rebuild the iterate with a fresh change-tag, so the cached derived quantities (curr_f, slacks, σ, …) invalidate correctly and the next step is computed from the new point — exactly as if the line search had produced it.

Staging a solver option (validated against the registry):

set opt mu_strategy adaptive
set opt linear_solver ma57

Staged options are not applied to the strategies already built for the running solve (they don’t re-read options mid-iteration). They take effect on a resolve or the next solve.

The read-side counterpart is get opt <name>, which reports an option’s current (or staged) value and its registry metadata — so you can confirm what a set opt actually staged before you resolve:

get opt mu_strategy        # → mu_strategy = adaptive  (staged)

Discovering options

opt                  # list every registered option
opt mu               # filter by name/category substring
complete pri         # completion candidates for a prefix

opt <exact-name> also prints the long description. In the REPL, Tab completes command verbs, block names, metric names (after break if), checkpoint names (after stop-at), event names (after break on), option names (after set opt / opt), and filesystem paths (after load / sweep / save / source — directories get a trailing /). The same contexts are available programmatically via the complete <prefix…> command (JSON complete), so an agent or GUI can offer the same completions.

Time travel

Rewind (goto / restart)

The debugger snapshots the primal-dual state (x, s, multipliers, μ, τ) every iteration. goto rewinds to a captured iteration and stays paused so you can re-tune before resuming:

goto 3               # rewind to the start of iteration 3
restart              # rewind to the earliest snapshot

Caveat — this is a soft rewind. Only the primal-dual state is restored; strategy history (the filter, the adaptive-μ oracle, the quasi-Newton memory) is not rolled back. So continuing from a rewound point is “resume from here,” not a bit-exact replay of the original run.

Re-solve from a saved point

resolve re-runs the solve from the current x with any set opt edits applied — a primal warm start with new options. Use it for “what if I change mu_strategy from here?”:

pounce-dbg> set opt mu_strategy adaptive
pounce-dbg> resolve
re-solving from current x with 1 staged option override(s)…
── pounce-dbg ── iter 0 @iter_start ...   # fresh solve, seeded from the captured x

Because each solve rebuilds its strategies from the options, the changes do take effect on the re-solve. The seed is dropped (falling back to the problem’s own start) if presolve / fixed-variable elimination changed the coordinate count.

Saving and visualizing artifacts

save                 # write the current iterate + residuals to a temp JSON
save /tmp/iter3.json # explicit path
viz x                # write a block and open it in an external viewer
viz dx               # a search-direction block
viz kkt              # the KKT inertia/regularization report

save writes every non-empty block, the search-direction blocks, and the residual scalars (iter, mu, objective, inf_pr, inf_du, nlp_error) — a self-contained artifact for external analysis.

load — the inverse of save

Typing a start point by hand is fine for a 2-variable toy and miserable for anything real. load reads a block straight into the live iterate, so you generate the point once (a prior solve, a surrogate, a sampler) and pull it in:

load /tmp/it0.json       # a `save` artifact: every block it contains is loaded
load start.csv           # a plain numeric file → x (comma/space/newline sep)
load start.csv s         # … into a named block instead of x

Two input shapes are accepted:

• A save artifact (JSON). Blocks are read from the top level or from an iterate object; every block present (x, s, multipliers, …) is written, each validated against the current dimension. So save→load round-trips a full point, and you can lift just the part that fits if dimensions changed.
• A plain numeric file — values separated by commas, whitespace, or newlines — written into the named block (default x). This is the many-variable escape hatch: numpy.savetxt("start.csv", x0) then load start.csv.

A loaded x becomes the seed for the next step (or for resolve — a warm start from an externally-computed point with no typing).

Interactive figures (pounce-dbg-viz)

viz writes a JSON artifact and hands it to a viewer. The Python package ships an interactive Plotly viewer that renders these properly — a spy/heatmap for viz kkt (the augmented matrix, colored by value, with the inertia/regularization in the title) and viz L (the LDLᵀ factor), and a bar chart for vector blocks (viz x, viz dx):

pip install 'pounce-solver[viz]'    # installs the `pounce-dbg-viz` script

When pounce-dbg-viz is on PATH, viz uses it automatically (opening an interactive figure in your browser). The launch order is:

1. $POUNCE_DBG_VIEWER — a command template ({} ← the artifact path), if set;
2. pounce-dbg-viz — the bundled Plotly viewer, if installed;
3. the OS opener (xdg-open / open) on the raw JSON.

So export POUNCE_DBG_VIEWER='python my_plot.py {}' overrides with your own plotter, and with nothing set + the viz extra installed it just works. The same pounce-dbg-viz <file.json> also renders a save artifact (the full iterate).

Multi-start and initialization sensitivity

Interior-point methods find a local solution, and which one depends on where you start. Two commands turn the debugger into an initialization-sensitivity probe: they run many full solves — each from a different start — and tabulate where each one ends up. Both build on the same re-solve machinery as resolve (so they need the restart cell the CLI wires by default; they error in contexts without it), and both leave you at a normal prompt on the final solve afterward.

sweep <file> — explicit starts

Run one solve per start point listed in a file (one start per line, comma/whitespace-separated; #/// comments and blank lines skipped):

pounce-dbg> sweep starts.txt
   sweep 1/4: Success                iters=21   obj=3.743990e-21 inf_pr=0.00e0
   sweep 2/4: Success                iters=15   obj=1.233088e-28 inf_pr=0.00e0
   sweep 3/4: Success                iters=14   obj=1.328861e-28 inf_pr=0.00e0
   sweep 4/4: Success                iters=29   obj=2.982346e-18 inf_pr=0.00e0

── sweep complete ── 4 solves, 4 succeeded, 1 distinct minima
     #  status                 iters       objective     inf_pr
     0  Success                   21    3.743990e-21     0.00e0
     1  Success                   15    1.233088e-28     0.00e0
     2  Success                   14    1.328861e-28     0.00e0
     3  Success                   29    2.982346e-18     0.00e0
   best: solve #2  obj=1.32886077e-28

Each start must have the same length as x (mismatches are reported with the line number). The summary clusters successful objectives to a relative 1e-6 to count distinct minima and flags the best (lowest-objective) solve. This is the “is this solve fragile to its start, and to which basins does it fall?” diagnostic — and unlike a black-box global search it leaves every solve’s trajectory observable: set a break on resto_entered or a stop-at kkt first and the sweep will pause inside whichever solve trips it.

multistart <N> [rel] — sampled restarts

When you don’t have a file of starts, multistart generates N of them:

pounce-dbg> multistart 8          # 8 starts
pounce-dbg> multistart 8 0.3      # wider jitter on any unbounded vars

Each variable that has a finite box [x_Lᵢ, x_Uᵢ] is sampled uniformly inside it — a genuine box multistart. Variables that are unbounded on either side fall back to a relative jitter ±rel·(|xᵢ|+1) around the current point (rel default 0.1, with a floor so components at zero still move). The command reports the split, e.g. multistart 8 (box 5/7 vars; 2 unbounded → jitter rel=0.1).

Start 0 is always the unperturbed current x (so the run includes where you already are), and the sampler is a fixed-seed PRNG, so a multistart run reproduces exactly.

The bounds are the ones the algorithm sees — full-length, post-scaling, after any bound_relax_factor — so every sampled start is a valid seed. For a problem with no finite bounds (a pure unconstrained NLP) multistart degrades to jitter around x; sweep an external sample if you want a specific spread there.

Driving a sweep from a file with load

The pieces compose. To seed a sweep from points computed elsewhere, write them with numpy.savetxt and sweep the file directly — or, for a single externally-computed warm start, load it and resolve:

import numpy as np
np.savetxt("starts.txt", sampler(n=32), delimiter=",")   # 32 starts, one per row

pounce-dbg> sweep starts.txt

sweep vs. find_minima

sweep/multistart are diagnostics: they show you how a handful of starts behave, with full visibility into each solve’s path. For an automated global search — Sobol sampling, deduplication, minimum certification (PSD Hessian), redundant-descent avoidance — reach for the Python pounce.find_minima, whose multistart and mlsl methods are the production tools. Rule of thumb: debugger sweep when you’re asking why a solve is start-sensitive; find_minima when you want the minima themselves.

Ask an LLM about the state

ask [question] packages the current paused state — checkpoint, residuals, step lengths, dimensions, and the KKT inertia/regularization — into a prompt and runs it through an LLM CLI (by default Claude Code, claude -p headless print mode), printing the reply inline. It’s AI-assisted debugging without leaving the loop:

pounce-dbg> stop-at kkt
pounce-dbg> continue
pounce-dbg> ask why is the dual infeasibility stalling?
# → the model's analysis of the state + suggested options to try

With no question it defaults to “explain the current state and suggest what to try next.”

Choosing the LLM

Set $POUNCE_DBG_LLM to pick the backend. It accepts either a bare provider keyword — which expands to that CLI’s correct non-interactive invocation — or a full command template:

export POUNCE_DBG_LLM=claude     # Claude Code      → claude -p   (default)
export POUNCE_DBG_LLM=codex      # OpenAI Codex CLI → codex exec <prompt>
export POUNCE_DBG_LLM=gemini     # Google Gemini    → gemini -p <prompt>
export POUNCE_DBG_LLM=llm        # simonw's llm     → llm <prompt>
# …or a full template:
export POUNCE_DBG_LLM='llm -m claude-opus'  # any prompt-on-stdin CLI
export POUNCE_DBG_LLM='mytool --ask {}'     # prompt substituted into {}

For a template, the prompt is fed on the tool’s stdin unless it contains a {} placeholder, in which case it is substituted as an argument. A bare word that isn’t a known provider is treated as a program name with the prompt on stdin.

Graceful when the CLI is absent. If the selected tool isn’t installed or on PATH, ask returns an error naming the tool and listing the provider keywords — the rest of the debugger (and the solve) is unaffected. ask is the only command that shells out; nothing else depends on an LLM being present.

In JSON mode the reply comes back in the result event’s data.reply.

Attaching to a run

You don’t have to single-step from iteration 0.

• Drop in on failure — --debug-on-error runs the solve freely and pauses at the terminated checkpoint only if the solve did not succeed, leaving you at the failing iterate for a post-mortem. (Plain --debug also pauses at terminated for a final-point inspect.)
• Attach with Ctrl-C — --debug-on-interrupt runs normally but installs a SIGINT handler; a first Ctrl-C drops you in at the next iteration (reason: "interrupt (Ctrl-C)"), a second Ctrl-C aborts. Ctrl-C also breaks into any other debug mode mid-continue.

Ctrl-C at the prompt. At a rustyline prompt Ctrl-C arrives as input, not a signal, so it has its own analogous double-tap: the first Ctrl-C cancels the current input line (readline convention), a second in a row stops the solve (a clean UserRequestedStop, same as quit). So whether you are running or sitting at the prompt, two Ctrl-Cs always get you out; quit/q and Ctrl-D (EOF, which detaches and finishes) remain the explicit exits.

Scripting

Run a sequence of debugger commands from a file — one per line, # and // comments and blank lines skipped:

# warmup.pdbg
break if inf_pr<1e-6
watch mu
stop-at after_search_dir
continue

pounce problem.nl --debug-script warmup.pdbg        # run at the first pause

pounce-dbg> source warmup.pdbg                       # or interactively

A script runs top-to-bottom and stops early if a command resumes or stops the solve (so ending with continue hands control back at the first breakpoint). --debug-script implies --debug when no --debug* mode is given, and runs once at the first pause (not on a resolve).

Example: a scripted initialization-sensitivity run

Because load, sweep, and set opt are ordinary commands, a whole diagnostic fits in a script file. This one watches each solve’s path and sweeps a set of externally-generated starts:

# sensitivity.pdbg — generate starts.txt first (e.g. numpy.savetxt)
break on resto_entered      # surface any start that falls into restoration
sweep starts.txt            # one solve per row; tabulated at the end

pounce model.nl --debug-script sensitivity.pdbg

Or compare a baseline against a what-if on the same starts by staging an option before the sweep:

# adaptive-vs-monotone.pdbg
set opt mu_strategy adaptive
multistart 16 0.2           # 16 sampled restarts, all under adaptive μ

Example: drive a multistart from a program (JSON protocol)

For many variables and many starts, hold the x0s as arrays in a driver program and let it assemble the commands — no point is ever typed. The --debug-json protocol emits a sweep_result per solve and a final sweep_summary:

import subprocess, json, numpy as np

p = subprocess.Popen(["pounce", "big.nl", "--debug-json"],
                     stdin=subprocess.PIPE, stdout=subprocess.PIPE, text=True)
send = lambda c, **k: (p.stdin.write(json.dumps({"cmd": c, **k}) + "\n"), p.stdin.flush())

recv = lambda: json.loads(p.stdout.readline())
recv()                                   # hello
recv()                                   # initial pause

# Option A — let the debugger sample: N restarts (uniform in finite boxes).
send("multistart", args=["32", "0.25"])

# Option B — supply your own starts via a file and sweep it:
# np.savetxt("starts.txt", my_sampler(n=32), delimiter=",")
# send("sweep", args=["starts.txt"])

results = []
for line in p.stdout:
    ev = json.loads(line)
    if ev.get("event") == "sweep_result":
        results.append((ev["status"], ev["objective"]))
    elif ev.get("event") == "sweep_summary":
        print(f"{ev['succeeded']}/{ev['solves']} ok, "
              f"{ev['distinct_minima']} distinct minima, "
              f"best obj {ev['best_objective']:.6e}")
        break

Each sweep_result carries index, status, iters, objective, inf_pr, and the seed it started from; the sweep_summary adds distinct_minima, best_index, and best_objective. A client can feature-detect support via hello.capabilities.sweep.

Exit model

Every non-kill path ends with a terminated event in JSON mode.

Command reference

The JSON protocol

--debug-json makes stdout a pure stream of newline-delimited JSON objects (the banner, problem stats, and final summary are routed to stderr, and print_level is forced to 0). A program reads one JSON object per line.

For an LLM agent: the whole contract

You do not need this page to drive the debugger — the protocol is self-describing. The contract is five lines:

1. Launch pounce <model> --debug-json (or --problem <name>), with the child’s stdin and stdout piped.
2. Read the first line — hello. It enumerates everything you can do: commands (the verbs), events (breakpoint triggers), checkpoints (where you can pause), metrics (the scalar field names), blocks (the inspectable vectors), and a capabilities map. Feature-detect off these lists, never off the version string.
3. Send commands, one JSON object (or bare string) per line, e.g. {"cmd":"break if inf_pr<1e-6","id":1} then {"cmd":"continue","id":2}. Set id to correlate the matching result.
4. Read events until you see the one you want. Every pause / progress / terminated event carries the same scalar metric fields, under the exact names listed in hello.metrics (objective, mu, inf_pr, inf_du, nlp_error, complementarity, iter) — so you can index them directly.
5. Finish with {"cmd":"continue"} to run to completion (then read terminated), or {"cmd":"quit"} to stop early.

A complete minimal transcript (→ sent, ← received), eliding long lines:

←  {"event":"hello","protocol":"pounce-dbg/1","commands":[…],"metrics":[…],…}
←  {"event":"pause","checkpoint":"iter_start","iter":0,"objective":24.2,…}
→  {"cmd":"break if inf_du<1e-6","id":1}
←  {"event":"result","request_id":1,"command":"break","ok":true,…}
→  {"cmd":"continue","id":2}
←  {"event":"progress","iter":1,"objective":4.7,"inf_du":2.1e1,…}
   … more progress events …
←  {"event":"pause","checkpoint":"iter_start","iter":21,"inf_du":8.7e-7,"reason":"inf_du<1e-6"}
→  {"cmd":"continue","id":3}
←  {"event":"terminated","status":"SolveSucceeded","iterations":21,…}

If you are wired in through the pounce-studio MCP server, you don’t even spawn the CLI yourself: call debug_start to open a live session and debug_command to step it (debug_state / debug_sessions / debug_close round it out) — the server owns the child process and the framing, and debug_start hands you the same hello handshake. Call debug_session_guide for the contract and a launch snippet if you’d rather drive --debug-json directly. The MCP analysis tools (diagnose, find_stalls, …) are post-mortem over a finished report; the debug_* tools and --debug-json are the live loop.

Session lifecycle

1. hello — emitted once, up front. The handshake.
2. pause — at each stop.
3. result — one per command, echoing the client’s request_id.
4. progress — one per iteration while running between pauses.
5. sweep_result / sweep_summary — during a sweep/multistart: one sweep_result per completed solve, then a sweep_summary at the end.
6. terminated — once, after the solve.

Commands

Write one per line to stdin, either a bare string or an object:

{"cmd": "print", "args": ["x"], "id": 7}
{"cmd": "break if inf_pr<1e-6", "id": 8}
"continue"

id (any JSON value) is echoed back as request_id on the matching result, for async correlation.

hello

{"event":"hello","protocol":"pounce-dbg/1","pounce_version":"0.4.0",
 "capabilities":{"inspect":true,"mutate_iterate":true,"mutate_mu":true,
   "conditional_breakpoints":"compound","request_ids":true,
   "viz":["block","delta","kkt","L"],"save":true,"load":true,"sweep":true,
   "kkt_inspect":true,"diagnose":true,"llm_assist":true,
   "pause_command":true,"equations":false,"structural_diagnose":false,
   "rewind":"primal_dual","resolve":true,"terminal_checkpoint":true,
   "interruptible":true,"progress_events":true,"async_pause":"checkpoint"},
 "checkpoints":["iter_start","after_mu","after_search_dir","after_step",
                "step_rejected","pre_restoration_entry",
                "post_restoration_exit","terminated"],
 "events":["resto_entered","resto_exited","regularized","tiny_step",
           "ls_rejected","mu_stalled","nan"],
 "commands":[…],"blocks":[…],"metrics":[…]}

A client should feature-detect off capabilities / checkpoints / events rather than the protocol string — those lists are additive as the debugger grows. A few capabilities are model-conditional: equations and structural_diagnose are true only when the solve came from an .nl file (which carries the source algebra and structural metadata) and false for a built-in problem, as shown above.

The handshake above is the NLP filter-IPM’s. Capabilities are answered for the backend that is actually running, so they agree with what the REPL will do: on the convex / conic IPM the backend-conditional entries (kkt_inspect, diagnose, mutate_mu, resolve, sweep, load, structural_diagnose) are false, viz is ["block","delta"], and blocks names that solver’s iterate blocks (x/s/y/z, plus tau/kappa on the HSDE drivers) — see the capability matrix below.

pause

{"event":"pause","checkpoint":"iter_start","status":null,
 "iter":3,"mu":2.0e-2,"objective":5.05,"inf_pr":0.0,"inf_du":2.7e-14,
 "nlp_error":0.0237,"complementarity":1.9e-2,"dims":{"x":2,"s":0,"y_c":0,
 "y_d":0,"z_l":2,"z_u":2,"v_l":0,"v_u":0},"breakpoints":[],"conditions":[],
 "reason":"mu<0.05"}

status is non-null only at the terminated checkpoint. reason carries the firing breakpoint / condition / event / interrupt.

result

{"event":"result","request_id":7,"command":"print x","ok":true,
 "output":["x = [-1.18e0, 1.38e0]"],"data":{"name":"x","values":[-1.18,1.38]}}

output is human-readable lines; data is the structured payload (present for inspection commands).

progress

{"event":"progress","iter":42,"mu":1.0e-5,"inf_pr":3.2e-7,"inf_du":1.1e-6,
 "objective":12.34,"nlp_error":1.1e-6,"complementarity":9.0e-6}

Emitted once per outer iteration during a continue, so a UI can show live progress instead of a hang. Carries the same scalar fields, under the same names, as pause — so hello.metrics names index directly off either event. Default on; toggle with the progress command.

terminated

{"event":"terminated","status":"SolveSucceeded",
 "status_message":"Optimal Solution Found.","iterations":6,
 "objective":4.9999999,"evals":{"obj":7,"obj_grad":7,"constr":1,
 "constr_jac":12,"hess":6}}

Async pause

A running continue can be interrupted two ways, both pausing at the next checkpoint with a reason:

• SIGINT — process.kill(pid, "SIGINT") (or Ctrl-C). This is what a Debug Adapter’s pause button maps to. Reason: "interrupt (Ctrl-C)".
• In-band command — send {"cmd":"pause"} on stdin while the solve is running (JSON mode). No signals, so it works on Windows. Reason: "pause (requested)".

hello.capabilities.async_pause is "checkpoint", and pause_command is true.

Tutorials

1. Why did this problem go to restoration?

$ pounce hard.nl --debug-json
{"cmd":"break on resto_entered"}
{"cmd":"continue"}
# → pause at checkpoint "pre_restoration_entry", reason "event: resto_entered"
{"cmd":"info"}                # how infeasible is the iterate?
{"cmd":"print kkt"}           # was the KKT singular / heavily regularized?
{"cmd":"print x"}

2. Catch a step that gets regularized

break on regularized
continue
# → pause at after_search_dir when delta_w > 0
print kkt        # inertia n- vs expected; delta_w / delta_c
print dx         # the (stabilized) Newton step

3. What-if: try a different μ strategy from here

break 5
continue                 # stop at iteration 5
set opt mu_strategy adaptive
resolve                  # re-solve from the iter-5 point with adaptive μ

4. Post-mortem on a failure

pounce maybe-infeasible.nl --debug-on-error

Runs unattended; if the solve returns anything but success you land at the final iterate:

── pounce-dbg ── TERMINATED (LocalInfeasibility)  iter 11  obj=1.13e0  inf_pr=5.0e-1  inf_du=1.2e-8
pounce-dbg> print x
pounce-dbg> print kkt

5. Drive it from a program / agent

import subprocess, json
p = subprocess.Popen(["pounce", "hs071.nl", "--debug-json"],
                     stdin=subprocess.PIPE, stdout=subprocess.PIPE, text=True)

def send(cmd, **kw): p.stdin.write(json.dumps({"cmd": cmd, **kw}) + "\n"); p.stdin.flush()
def recv():          return json.loads(p.stdout.readline())

hello = recv()                       # capabilities / vocabulary
print(recv())                        # initial pause
send("break if inf_du<1e-6", id=1)
print(recv())                        # result, request_id=1
send("continue")
for line in p.stdout:                # progress … pause … terminated
    ev = json.loads(line)
    if ev["event"] == "terminated": break

6. Is this solve sensitive to its start?

break on resto_entered       # flag any start that falls into restoration
multistart 16                # 16 restarts (uniform in each finite box)
# → per-solve lines, then a table: succeeded / distinct minima / best

Swap multistart 16 for sweep starts.txt to run your own start points (numpy.savetxt("starts.txt", X0, delimiter=",")). See Multi-start and initialization sensitivity.

Beyond the interior-point loop

Everything above is the NLP filter-IPM. The same debugger — same command engine, same REPL — drives the other solvers too.

Convex and conic solves

The convex LP/QP interior-point solver and the HSDE conic drivers (SOCP, the exponential / power cones, and small PSD cones) expose the same checkpoints and commands as the NLP loop. The iterate blocks follow the QP standard form — x (variables), s (cone slacks), y (equality multipliers), z (inequality / cone multipliers) — and the HSDE drivers additionally expose the homogenizing scalars tau / kappa as 1-element blocks (print tau). set <block> and goto work as on the NLP path; set mu is rejected, because the convex μ is derived from ⟨s, z⟩ (edit s/z to move it).

pounce model.nl --debug                 # LP / convex-QP (auto-routed) — IPM REPL
pounce_cblib model.cbf --debug          # SOCP / exp / power / PSD (conic) — IPM REPL
pounce_cblib model.cbf --debug-script s.pdbg

Capability matrix

The flow-control core — checkpoints, stepping, breakpoints, watchpoints, block/scalar inspection, diff, goto/restart, save, ask, and the JSON protocol — works identically on every backend. The table below is just the commands whose availability is backend- or model-conditional; anything not listed is universal. A command that isn’t available on the current backend returns an explicit “not available for this solver” error (it never silently no-ops), and a JSON client should feature-detect off hello.capabilities rather than this table.

✅ available · ⚠️ model-conditional · ➖ reports “not available for this solver” · ❌ explicitly rejected with an explanation

The streamed scalar metric vocabulary (iter, mu, objective, inf_pr, inf_du, nlp_error, complementarity) is the same on every backend — see hello.metrics. Each backend maps its native quantities onto these NLP-centric names; the convex IPM, for instance, reports nlp_error = max(pinf, dinf, μ). A backend that has no value for a metric reports it as JSON null (never a dropped field), and a test pins the emitted set to that single advertised vocabulary so it can’t drift.

A third backend — an interactive branch-and-bound tree debugger for a spatial global solver — is not part of this release.

Limitations

• Soft rewind only. goto/restart restore the primal-dual state, not strategy history (see the caveat above).
• set opt is staged, not hot-applied to a running solve; it takes effect on resolve / the next solve.

1. A. Lee, R. B. Parker, S. Poon, D. Gunter, A. W. Dowling, and B. Nicholson, “Model Diagnostics for Equation-Oriented Models: Roadblocks and the Path Forward,” Systems and Control Transactions 3:966–974 (2024). https://doi.org/10.69997/sct.147875 ↩ ↩2 ↩3 ↩4

Pyomo

Because POUNCE speaks the AMPL NL/SOL protocol, it drops into Pyomo through the AMPL Solver Library interface — exactly how Pyomo drives Ipopt.

The pyomo-pounce package registers pounce as a Pyomo SolverFactory solver:

import pyomo_pounce  # registers 'pounce'
from pyomo.environ import ConcreteModel, Var, Objective, SolverFactory

model = ConcreteModel()
model.x = Var(bounds=(-10, 10))
model.obj = Objective(expr=(model.x - 3) ** 2)

solver = SolverFactory('pounce')
solver.solve(model)

Options pass through the usual Pyomo mechanism:

solver.solve(model, options={'tol': 1e-10, 'max_iter': 500})

Under the hood, Pyomo writes the model to an AMPL .nl file, invokes pounce problem.nl -AMPL, and reads the result back from the .sol file. See Running Solves for the -AMPL solver mode.

Which pounce binary runs

import pyomo_pounce is required before SolverFactory('pounce'). Without it Pyomo does not know the solver and raises a clear UnknownSolver / “plugin not registered” error — it does not silently run some other pounce. With it imported, the plugin runs the binary bundled in the pounce-solver wheel, independent of PATH; only a source/dev install lacking that wheel falls back to a pounce on PATH (and the plugin warns when it does).

Because two builds can report the same version string (X.Y.Z) while behaving differently — a binary from before and after a fix does — a stale pounce on PATH is otherwise hard to notice. To see exactly which executable will run, its build (the git commit from pounce --about), and whether a different pounce earlier on PATH would shadow it:

import pyomo_pounce
pyomo_pounce.check_binary()   # prints a report; returns a dict

Which interface runs — and why it matters for timing

Pyomo has more than one way to drive an NL/SOL solver, and they are genuinely different code paths, not aliases. All of these reach POUNCE (verified against Pyomo 6.10.1):

The first three are pyomo-pounce’s own registrations and the supported routes; they are the only ones that bring the rest of this page with them: the scaling_factor suffix handling, the sensitivity path, the preflight/repair helpers, the guard against handing a model with live integer variables to a continuous solver, and the bundled-binary resolution above. The generic routes run the same solver and return the same answer, but silently do without all of that.

ipopt and asl are Pyomo’s legacy solver interface; ipopt_v2 is the newer pyomo.contrib.solver one. Driving POUNCE through ipopt_v2 needs a build carrying the two ASL-compatibility fixes noted in the CHANGELOG under “Pyomo’s modern solver interface could not drive POUNCE at all” — before them it failed on every model, because Pyomo v2 passes options as key="value" in a single argv entry (quotes and all, since no shell is involved) and because POUNCE’s .sol wrote an Options count of 0, which the v2 .sol reader rejects.

Choosing between the legacy and v2 interfaces

Both of pyomo-pounce’s interfaces carry the same extras and return the same numbers — a test in pyomo-pounce/tests/test_v2.py solves one model through both and compares primals, objective, duals and reduced costs, so this is checked on every CI run rather than asserted here.

import pyomo_pounce
from pyomo.environ import SolverFactory
from pyomo.contrib.solver.common.factory import SolverFactory as SolverFactoryV2

solver = SolverFactory('pounce')       # legacy interface
solver = SolverFactoryV2('pounce')     # v2 interface (a Results object)
solver = SolverFactory('pounce_v2')    # v2 engine, legacy-style API

The v2 route needs Pyomo ≥ 6.10.1 (where the SolutionLoader / get_vars API it builds on landed — pyomo.contrib.solver.common exists from 6.9.2, but 6.9.2–6.10.0 ship the older SolutionLoaderBase / get_primals) and pounce-solver > 0.9.0 (Pyomo’s asl_sol_reader is strict where the legacy reader is lenient and needs the per-model .sol Options echo added after 0.9.0). pip install pyomo-pounce[pyomo-v2] asks for both. Neither applies to SolverFactory('pounce'): on an older Pyomo the legacy plugin works exactly as before and pyomo_pounce.HAVE_V2_INTERFACE reports False.

They differ in API and in per-solve overhead. The v2 interface returns a Results object and hands the solution back through a solution loader (so load_solutions=False gives you the values without touching the model); the legacy one returns a SolverResults and loads into the model as a side effect. Options are solver_options={...} on v2 against options={...} on the legacy route.

The v2 route can be materially faster outside the solve, and how much depends on the model’s shape. Same POUNCE binary, wall clock around solve() minus POUNCE’s own reported time:

On the plain model the two are indistinguishable; on the IDAES-shaped ones the legacy interface adds roughly 0.25 s per solve (~1.8×). If your models are of that kind and you are solving many of them, the v2 route is worth taking. (Figures from the #552 measurements; the first row was measured on Linux, the other two on Windows, so read down the columns rather than across the rows.)

If you are benchmarking, put both solvers on the same interface. solver.solve(model) is not only the solve: it is Pyomo writing the .nl, launching the process, reading the .sol back and loading it into the model. Timing around that call and subtracting the solver’s own reported time leaves a remainder that is mostly Pyomo’s work, and — as the table above shows — it is not the same work on every interface. On the 3,010-variable collocation model that remainder breaks down as ~0.082 s Pyomo writing the .nl, ~0.020 s process spawn plus POUNCE’s own .nl read and setup, and ~0.008 s Pyomo reading the .sol and loading it. So comparing SolverFactory('pounce') against SolverFactory('ipopt_v2', …) compares two Pyomo interfaces as much as two solvers, and attributing the remainder to either solver’s file handling will mislead you. Use the same interface on both sides, or compare the solvers’ own reported times.

User scaling with the scaling_factor Suffix

A badly conditioned model converges poorly, and often you know its natural units better than the solver can infer from gradients at x0. The standard Pyomo channel for saying so is the scaling_factor Suffix, read exactly as Ipopt reads it:

model.scaling_factor = Suffix(direction=Suffix.EXPORT)
model.scaling_factor[model.obj] = 1e-3           # objective in MW, not W
model.scaling_factor[model.mass_balance] = 1e2   # one constraint
model.scaling_factor[model.energy_balance] = 1e2 # or a whole container

solver.solve(model, options={'nlp_scaling_method': 'user-scaling'})

Both halves are required: without nlp_scaling_method=user-scaling the Suffix is inert (a scaling_factor Suffix also drives Pyomo’s own core.scale_model transformation, which never involves the solver), and without the Suffix the option has nothing to apply — pyomo-pounce warns in that case rather than leaving you to wonder.

Rules, matching AMPL/Ipopt:

• Only an export-enabled Suffix counts (Suffix.EXPORT or Suffix.IMPORT_EXPORT).
• Components you do not list are unscaled, as are components listed with a factor of 0.
• An entry on a container applies to every member.
• Entries on inactive constraints/objectives and on fixed variables are skipped — none is a row or column of the problem the solver is handed.
• Scaling changes conditioning, never the answer: solutions, duals, and everything the sensitivity accessors report come back in your model’s units.

Variables can be scaled, and a factor on a Var is applied as a change of variables inside the solver: the algorithm works in the scaled coordinates and the solution, the duals, and the bound multipliers come back in your model’s own units. No clone of the model is made and no propagate_solution step is needed, which is what distinguishes this from Pyomo’s core.scale_model transformation. Factors must be positive and finite. A negative factor would reverse a variable’s direction and swap its bounds, so it raises rather than being applied.

This works on both solve paths: the ordinary ASL/subprocess solve and the in-process path taken when the model carries sensitivity declarations — including the accessors themselves. covariance(), information(), gradient() and estimate() read the solver’s KKT factorization directly rather than through the scaling layer, so they carry the factors through their own natural-units translation and answer in your model’s units on a variable-scaled solve (issue #486).

Preflight and initialization

A Var whose .value was never set is written as 0 into the .nl file, so an uninitialized model actually starts at the origin (see Initialization and Warm Starts). The package ships a preflight check plus an initialization pipeline for exactly this:

import pyomo_pounce

report = pyomo_pounce.preflight(model)   # what will POUNCE see at x0?
print(report)                            # unset vars, bound/constraint
if report.fatal:                         # violations, NaN/inf evaluations
    ...

# fill -> repair -> block-solve, with the decisions held constant:
rep = pyomo_pounce.initialize(model, decisions=[m.feed, m.reflux])
if not rep.block.square:
    print(rep)          # names of what you forgot to specify

preflight evaluates every active constraint and the objective at the current values with unset values treated as 0 (exactly what the NL writer sends), restores the model untouched, and reports what iteration 0 will see; report.fatal means the solve would abort with Invalid_Number_Detected.

initialize follows the workflow you would run by hand on, say, a distillation column: set the decisions (feed, reflux, boilup), solve for a physical profile with them held constant, then let the optimizer move them. Its three stages are also available individually:

pyomo_pounce.initialize_missing_values(model)   # bounds-aware fill
                                                # (midpoint / one unit
                                                # inside / zero)

pyomo_pounce.project_to_feasible(model)         # min-norm repair: move the
                                                # current point onto the
                                                # model's own constraints
                                                # (one POUNCE solve)

rep = pyomo_pounce.block_initialize(            # solve the equality
    model, decisions=[m.feed, m.reflux])        # system's square blocks
                                                # in calculation order

initialize_missing_values fills each variable independently, so the fill can be internally inconsistent (mole fractions that do not sum to one); project_to_feasible repairs that by minimizing sum((v - v0)**2) subject to the model’s active constraints and bounds — the full nonlinear projection, solved with POUNCE, with the original objective restored afterwards.

Both stages guarantee that a failed solve leaves variable values exactly as they were: a diverged projection restores the pre-projection point, and a failed block solve restores that block’s seeds and stops, so initialization can never make your starting point worse than it found it.

block_initialize is IDAES-flavored initialization without hand-written routines. decisions= holds the listed variables at their current values for the solve and releases them afterwards (each must have a value). The active equality constraints are decomposed (Dulmage-Mendelsohn, via pyomo.contrib.incidence_analysis); the square part is solved block by block in topological order by Pyomo’s solve_strongly_connected_components (1x1 blocks by Newton, larger blocks by POUNCE), filling Var.value along the way. When the system is not square, report.square is False and the offending variables and constraints are reported by name — underconstrained_variables is the list of things you forgot to specify or flag as decisions, overconstrained_constraints the redundant or conflicting specifications. Permanently-known inputs can simply be fix()ed instead of listed as decisions.

The analysis half is also available on its own:

rep = pyomo_pounce.block_analyze(               # the DM partition only:
    model, decisions=[m.feed, m.reflux])        # nothing seeded or solved
rep.underconstrained_variables                  # VarData objects, uncapped
rep.n_extra_degrees_of_freedom                  # how many specs are missing
rep.variable_blocks                             # the calculation order

block_analyze runs the same decision handling and the same Dulmage-Mendelsohn decomposition, but touches nothing: no values are read or written (so, unlike block_initialize, the decisions do not need values), and no solve happens. Where the initialization reports cap their name lists for display, block_analyze returns the full partition as the component objects themselves: the underconstrained and overconstrained subsystems, the square part, and its block-triangular calculation order. Use it to diagnose a large model’s specification, or as the structural front end for tooling that decides what to specify before calling initialize / block_initialize to do the work.

Repairing a bad specification

Some specifications are structurally wrong, not just badly started. On a distillation column at steady state, holding all the flow controls leaves the drum levels undetermined while the holdup balances become redundant — square by count, singular in structure, and no starting point fixes that. block_repair_plan plans a valid specification instead of failing on the broken one:

plan = pyomo_pounce.block_repair_plan(
    model,
    decision_candidates=[m.LT, m.VB, m.D, m.B])  # what you would like held
plan.decisions   # candidates a square system can hold
plan.pruned      # candidates the equalities claim: solved for instead
plan.pinned      # what nothing determines: hold at values you choose

The candidates are pruned to the subset a valid specification can hold: matching prefers plain variables over candidates, which provably minimizes the number pruned, and among candidates earlier-listed ones are preferentially kept, so the listing order acts as an implicit priority when a pruning tie could go either way. The pins need no user input: a variable is pinned when every one of its edges is provably unusable — the key case being an equation 0 == f/g, which cannot determine a variable appearing only in the denominator g, since its sensitivity there vanishes at every solution. That is exactly the shape substituting d/dt = 0 into a dynamic balance produces, which is how loose integrators (drum levels with no weir feedback) hide in steady-state models. Like block_analyze it is a plan, not an action: nothing is fixed, read, or written, and no values are needed. loose_variables (undetermined, not repairable) and redundant_constraints (satisfiable by no specification) are genuine model defects.

initialize and block_initialize run the same check on their decisions automatically (repair="auto", the default). A square specification is used exactly as given, the shipped behavior. A broken one is repaired: the decisions become the candidate pool, conflicting ones are pruned (they need no values), pins are seeded bounds-aware and never at zero when valueless (a pin lives in denominators, so zero is the one forbidden seed), and report.repair records the plan (None when nothing was needed). Pass repair="off" for the strict path: decisions held exactly as given, and a non-square specification is reported (report.square, the name lists) instead of repaired. The repair is call-scoped exactly like the decisions themselves (fixed flags restored, values only), so it never changes your model’s own specification. To apply a plan to a model you intend to solve — a square simulation, say — fix plan.decisions and plan.pinned and leave plan.pruned free; which variables to fix is a modeling decision, so the plan leaves it to you.

GAMS

POUNCE plugs into GAMS as an NLP solver, so a model can hand its problem to POUNCE with:

option nlp = pounce;
solve mymodel using nlp minimizing obj;

There are two ways to make POUNCE available to GAMS. Pick one:

Both register POUNCE under the same name (pounce) for NLP, DNLP, and RMINLP models — POUNCE is a continuous local NLP solver, so mixed-integer and conic model types are not offered here.

The pip route

1. Install

pip install pounce-solver[gams]

The [gams] extra pulls in gamsapi[core] — GAMS’s own expert-level GMO/GEV Python bindings — and PyYAML. The bindings dlopen the GAMS C libraries from your local install, so gamsapi must match your GAMS version. POUNCE itself redistributes nothing GAMS-owned. If your GAMS and gamsapi versions disagree, install the matching one from your GAMS system (GAMS ships a gamsapi wheel under apifiles/Python/), or:

pip install 'gamsapi[core]==<your GAMS X.Y.Z>'

2. Check the install

pounce-gams status

reports whether gamsapi imports, the config directory POUNCE will register into, and whether POUNCE is already registered:

gamsapi:       available
               gamsapi 53.2.0 importable
config dir:    /Users/you/Library/Preferences/GAMS
gamsconfig:    /Users/you/Library/Preferences/GAMS/gamsconfig.yaml (missing)
POUNCE solver: not registered

3. Register

pounce-gams register

This writes a tiny launcher script and a solverConfig entry into your GAMS per-user gamsconfig.yaml. It merges — any other solvers already in that file (CONOPT overrides, discopt, …) are preserved — and is idempotent (re-running just updates POUNCE in place). The per-user config directory GAMS searches is OS-specific:

Override the target with --config-dir <path> (e.g. to register into the GAMS system directory instead). To undo, pounce-gams unregister.

No sudo is needed and nothing is written into the GAMS system directory, so a GAMS upgrade does not wipe the registration.

4. Solve

option nlp = pounce;
solve mymodel using nlp minimizing obj;

GAMS invokes the launcher with a control file; the launcher runs the Python link, which reads the model through GMO/GEV, solves it with POUNCE, and writes the primal/dual solution and GAMS model/solve status back.

Marginals

Equation marginals (.M) follow the usual GAMS convention: they are stated against the objective as you wrote it, for both minimizing and maximizing models. POUNCE always minimizes internally, so for a maximizing model it solves min(−f) and converts its multipliers back via pi = −obj_sign · λ (see gams_pi in python/pounce/gams/link.py and the matching block in gams/gams_pounce.c).

Variable marginals (.M on variables, i.e. reduced costs) are z_L − z_U, carrying the same obj_sign factor.

Before v0.9.1 the conversion omitted the obj_sign factor, so equation marginals on maximizing models came back with the wrong sign (#272). minimizing models, objective values, variable marginals and status mapping were never affected. Both the pip link and the native C link were affected identically, so the install method made no difference.

Option files

If a model sets mymodel.optfile = 1, POUNCE reads pounce.opt (.op2, .op3, … for higher optfile values). Each line is a keyword value pair using POUNCE’s option names; lines starting with * or # are comments. The GAMS iterlim and reslim are honored as max_iter and max_wall_time.

* pounce.opt
tol        1e-10
max_iter   500

Machine-readable solve report

Set json_output in pounce.opt to also emit a structured pounce.solve-report/v1 JSON report (identical to the CLI’s --json-output, consumable by pounce-studio):

json_output  my_solve.json
json_detail  full        * "summary" or "full"; default is "full"

See the JSON Solve Report schema for the format.

Notes & limitations

• Version match. The single most common failure is a gamsapi ↔ GAMS version mismatch; pounce-gams status diagnoses it.
• Warm starts. The active-set-SQP working-set / state-file warm-start features (algorithm active-set-sqp, sqp_state_file) are currently only in the native C link, where each solve reuses an in-process state. The pip link runs each solve as a fresh process; full warm-start parity is a planned follow-up.

Python API

POUNCE ships a Python wrapper that is intentionally cyipopt-compatible: code written for cyipopt typically runs against POUNCE by changing only the import.

Install

cd python
pip install maturin
maturin develop --release    # builds the native extension into your venv

Optional extras:

pip install -e .[jax]        # JAX integration
pip install -e .[torch]      # PyTorch integration
pip install -e .[dev]        # tests + jax + torch + scipy

cyipopt-style interface

import numpy as np
import pounce

class HS071:
    def objective(self, x):
        return x[0]*x[3]*(x[0]+x[1]+x[2]) + x[2]
    def gradient(self, x):
        return np.array([
            x[0]*x[3] + x[3]*(x[0]+x[1]+x[2]),
            x[0]*x[3],
            x[0]*x[3] + 1.0,
            x[0]*(x[0]+x[1]+x[2]),
        ])
    def constraints(self, x):
        return np.array([np.prod(x), np.dot(x, x)])
    def jacobianstructure(self):
        return (np.repeat([0, 1], 4), np.tile([0, 1, 2, 3], 2))
    def jacobian(self, x):
        return np.array([
            x[1]*x[2]*x[3], x[0]*x[2]*x[3], x[0]*x[1]*x[3], x[0]*x[1]*x[2],
            2*x[0], 2*x[1], 2*x[2], 2*x[3],
        ])

prob = pounce.Problem(
    n=4, m=2,
    problem_obj=HS071(),
    lb=[1]*4, ub=[5]*4,
    cl=[25, 40], cu=[2e19, 40],
)
prob.add_option('tol', 1e-8)
x, info = prob.solve(x0=np.array([1.0, 5.0, 5.0, 1.0]))
print(info['status_msg'], info['obj_val'], x)

Verifying convergence / trustworthy duals

info carries the final KKT residuals so a consumer can independently check how converged a returned point is — useful when the duals (info["mult_g"], info["mult_x_L"], info["mult_x_U"]) feed a downstream certificate (e.g. dual bound tightening). Two flavors:

• final_kkt_error / final_dual_inf / final_constr_viol / final_compl — the residuals the convergence test saw, in the internally-scaled NLP space (the nlp_scaling_method factors).
• final_unscaled_kkt_error / final_unscaled_dual_inf / final_unscaled_constr_viol / final_unscaled_compl — the same residuals with the scaling divided back out, i.e. in your original problem units. Equal to the scaled values when no scaling activates.

On ill-conditioned problems nlp_scaling can deflate the scaled residual enough that the default test reports Solve_Succeeded while the unscaled duals have drifted. The robust guard is to read the residual yourself rather than trust the status enum alone — info["status"] is a coarse signal, and some callers treat Solve_Succeeded (0) and Solved_To_Acceptable_Level (1) identically:

x, info = prob.solve(x0=...)
converged = info['final_unscaled_kkt_error'] <= 1e-6   # your own threshold

If you’re feeding the duals into a downstream certificate (e.g. a dual bound), prefer building a safe bound from the multipliers — valid for any dual-feasible point, so it doesn’t hinge on the solver’s exact convergence.

Two convenience knobs back this up:

• Tighten the (unscaled) component tolerances — dual_inf_tol, constr_viol_tol, compl_inf_tol gate on the unscaled residuals, so the solver keeps iterating until it actually meets them (or exits non-success).
• kkt_fidelity_tol (default 0 = off) — a defensive post-solve relabel: a Solve_Succeeded whose final_unscaled_kkt_error exceeds it is demoted to Solved_To_Acceptable_Level. Note this only helps a caller that distinguishes those two statuses; if yours doesn’t, gate on the residual directly as above.

Where the time went (info["timing"])

Every Problem.solve attaches a per-subsystem wall-clock breakdown so you can attribute a solve’s runtime without patching or rebuilding the solver. info["wall_time"] is the overall-algorithm total (seconds); info["timing"] is a dict of the same total plus its components:

x, info = prob.solve(x0=...)
t = info["timing"]
print(t["overall_alg"])                    # total solve wall time
print(t["linear_system_factorization"],    # KKT factorization vs …
      t["linear_system_back_solve"],        # … back-solve, and their
      t["linear_system_total"])             # sum (total linear algebra)
print(t["eval_objective"], t["eval_gradient"],
      t["eval_constraints"], t["eval_constraint_jacobian"],
      t["eval_lagrangian_hessian"])         # per-callback eval time

The scipy-style pounce.minimize facade mirrors these onto the result as res.wall_time and res.timing (also in res.info). The callback split is what lets you see, for example, that a reduced-space / variable-aggregation solve converges in few iterations but spends most of its time in a densified Lagrangian-Hessian evaluation — the func/Jacobian/Hessian story becomes a direct measurement rather than an inference. All values are wall-clock seconds; unused subsystems read 0.0.

Caller-supplied KKT ordering (set_ordering)

A structure-aware presolve can hand pounce a fill-reducing permutation for the KKT linear solver that the built-in AMD/METIS pass cannot derive — a block-triangular / Schur ordering (Parker, Garcia & Bent, arXiv:2602.17968) or a tearing ordering from equation-oriented decomposition. Install it on the low-level Problem before solving:

prob = pounce.Problem(n, m, problem_obj=...)
prob.set_ordering(perm)        # 0-based new-to-old permutation (list / int array)
x, info = prob.solve(x0=...)
# prob.get_ordering()  -> the installed permutation, or None
# prob.clear_ordering() -> restore the feral_ordering default

perm[k] is the original index that becomes index k. Its length must equal the augmented KKT system dimension (variables + slacks + constraint duals), not the problem’s n; for an unconstrained problem that is n, but with constraints it is larger. The ordering is validated inside FERAL as a bijection — a wrong length or a duplicate fails the factorization and the solve returns a non-success status (e.g. Error_In_Step_Computation) rather than crashing or returning a wrong answer, since a permutation only affects fill and pivot order, never the computed solution. set_ordering is persistent config (it applies to every subsequent solve() until clear_ordering()) and is honored only by the default FERAL backend. This maps to FERAL’s OrderingMethod::External (feral#107).

Block-triangular / Schur KKT solve (set_kkt_schur_block)

If a presolve can identify a reducible block of the KKT system — e.g. the nonsingular block-triangular submatrix a reduced-space / variable-aggregation analysis exposes (Parker, Garcia & Bent, arXiv:2602.17968) — it can hand that block to pounce, which Schur-complements it out and factorizes only the two diagonal blocks, recovering the full-system inertia a priori via Sylvester’s law:

prob = pounce.Problem(n, m, problem_obj=...)   # needs an exact Hessian
prob.set_kkt_schur_block(indices)              # KKT-space indices of the Schur block
x, info = prob.solve(x0=...)
# prob.get_kkt_schur_block()  -> installed indices, or None
# prob.clear_kkt_schur_block()

indices are KKT-space indices into 0..dim where dim = n + n_slack + n_eq + n_ineq, in the solver’s internal x, slack, eq-dual, ineq-dual block order (e.g. for an all-equality problem the constraint-dual block is range(n, n + n_eq), and the primal block is the positive-definite eliminated block — the classic range/null-space split). The method wins only when the Schur block is much smaller than the eliminated block (the dense Schur complement is O(n_schur²) to store and O(n_schur³) to factor). When the partition is unsuitable — too large a fraction of the system, malformed, or a diagonal block turns out singular — the solver falls back to the standard full-space path transparently, so the hook can never break a solve; it only changes how the identical system is factored, never the solution. Honored on the default feral + exact-Hessian path.

Building a model in memory (NlExpr / build_nl_problem)

pounce.read_nl("model.nl") gives you pounce’s native reverse-mode-AD evaluators for an AMPL .nl file on disk. Two sibling entry points reach the same machinery without a file:

• pounce.parse_nl_text(text, var_names=None, con_names=None) — the same parser, fed a string. For a frontend that already generates .nl, this drops the temp file and its cleanup. There are no sibling .col / .row files to read, so names are passed explicitly.
• pounce.build_nl_problem(...) — skip .nl entirely and hand over expression trees built from pounce.NlExpr.

Both return the same NlProblem class read_nl does, with the same surface: objective, gradient, constraints, jacobian / jacobian_structure, hessian / hessian_structure, hessian_vector_product, and variant. They also feed solve_nlp_batch.

import pounce

x = pounce.NlExpr.vars(2)                      # [Var(0), Var(1)]
rosen = (1 - x[0]) ** 2 + 100 * (x[1] - x[0] ** 2) ** 2

p = pounce.build_nl_problem(
    n=2,
    objective=rosen,
    constraints=[x[0] ** 2 + x[1] ** 2],
    g_l=[0.0], g_u=[2.0],
    x0=[-1.2, 1.0],
)

p.objective(p.x0)          # float
p.gradient(p.x0)           # ndarray[n]
(x_star, info), = pounce.solve_nlp_batch([p])

Bounds default to unbounded (±1e19, the .nl sentinel) and x0 to zeros. minimize=False maximizes; as with a parsed maximize model, the returned objective/gradient/Hessian are negated so that minimizing them solves the model, and p.minimize records the original sense.

NlExpr supports the Python arithmetic operators (+ - * / ** -, abs, with plain numbers accepted on either side) plus method-form transcendentals: sqrt exp log log10 sin cos tan asin acos atan sinh cosh tanh asinh acosh atanh erf. Multi-argument and control-flow nodes are static methods: NlExpr.sum(iterable), NlExpr.atan2(y, x), NlExpr.min(*args), NlExpr.max(*args), NlExpr.compare(op, a, b), NlExpr.select(cond, then_, else_), and NlExpr.logical_and / logical_or / logical_not.

Comparison is spelled NlExpr.compare("<", a, b) rather than a < b: overloading Python’s comparison operators would break every ordinary use of an expression in a container. The result is piecewise constant (zero derivative), and pairs with NlExpr.select.

Why not just write .nl? Because the round trip is lossy. .nl writers commonly refuse atan2 (no two-argument funcall path) and min/max (they force a DNLP model type), and AMPL has no erf opcode at all — yet pounce’s tape differentiates all three natively. Built here, they survive:

x = pounce.NlExpr.vars(2)
p = pounce.build_nl_problem(n=2, objective=pounce.NlExpr.sum([
    pounce.NlExpr.atan2(x[0], x[1]),
    pounce.NlExpr.min(x[0], x[1]),
    x[0].erf(),
]))

Operands are shared, not copied. a * b references its operands rather than deep-copying them, so building an expression costs the same whether the pieces are two variables or two half-million-node models. Two consequences worth knowing:

• Accumulating in a Python loop is linear in the number of terms. It still nests one level deeper per term, though, and nesting is capped (below) — so a many-term sum still belongs in one flat NlExpr.sum(terms) node, which tapes better and is one level whatever its length. The same goes for min / max, flat in their argument count too.
• Reusing a Python name reuses the subexpression. t = x[0] * x[1] used in ten places is one shared body on the tape, evaluated once per sweep, with its adjoint summing the ten contributions — the same value and the same derivatives as writing it out ten times, off a tape a tenth the size. Expressions that are only tractable as a DAG work too: for _ in range(40): e = e * e is 40 nodes describing x ** 2**40, and it builds, tapes, and differentiates in under a millisecond.

e = pounce.NlExpr.const_(0.0)
for t in terms:                    # linear, but len(terms) levels deep
    e = e + t

e = pounce.NlExpr.sum(terms)       # linear and one level — prefer this

Nesting is capped at NlExpr.max_depth (10 000), and exceeding it raises ValueError. Every consumer of an expression — the tape builder, the problem assembler, freeing it, and the .nl parser that produces one — recurses once per level, so a deep enough tree overflows the stack, which is a hard crash rather than an exception. Two things keep that unreachable: those walks run on a worker thread with a 64 MB stack, so what is survivable does not depend on the calling thread (8 MB on a macOS/Linux main thread, 1 MB on Windows, less on a threading.Thread), and the cap then keeps the depth well inside it.

The same limit applies to read_nl and parse_nl_text, which enforce it on what they parsed rather than as it is built — a model that arrives already built cannot be capped during construction. A .nl file that spells a long sum as an o0 (binary +) chain rather than o54 (n-ary sum) is the way to hit it. The cap bounds nesting, not size: NlExpr.sum and o54 are one level regardless of term count, so wide models are unaffected. Each expression’s .depth is readable.

For checking a subexpression before wiring it into a model, NlExpr has .eval(x), .gradient(x), and .variables(), which build a one-off tape for that expression alone.

Two things NlExpr does not do: it cannot carry AMPL imported (external) functions — build_nl_problem has nowhere to put the F-segment declarations that bind them, so use read_nl / parse_nl_text for a model that needs them — and it cannot be pickled. copy.copy and copy.deepcopy do work.

Hessian-vector products

NlProblem.hessian_vector_product(x, v, lam=None, obj_factor=1.0) returns (obj_factor·∇²f + Σᵢ lamᵢ·∇²gᵢ) · v without ever forming the Hessian — one forward-over-reverse AD pass per tape, seeded with v directly. hessian(...) instead runs one such pass per Hessian color and decodes the compressed columns into the sparse lower triangle, so on a large model the matrix-free call is cheaper by roughly the chromatic number of the coloring. It is the operator a Newton–Krylov / truncated-CG step wants.

Hv = p.hessian_vector_product(x, v)                 # objective block only
Hv = p.hessian_vector_product(x, v, lam, 1.0)       # full Lagrangian

Available on every NlProblem, however it was built — read_nl, parse_nl_text, build_nl_problem, or variant.

Dense and sparse directions. v may be any of:

The shape rule is the same dense or sparse: (n,) or (n, k). A (1, n) row vector raises rather than being guessed at — for a square-ish block it is indistinguishable from k directions of the wrong length. Watch for this with SciPy sparse matrices, which shape a 1-D input as a row: csr_matrix(v) on a length-n v is (1, n) and will be refused. Pass v[:, None], or use the 1-D sparse array API — coo_array(v) is genuinely (n,), on SciPy >= 1.14.

import scipy.sparse as sp

p.hessian_vector_product(x, sp.csc_matrix(V))   # sparse block of directions
p.hessian_vector_product(x, np.eye(n))          # densify: H, in one call

A sparse v is densified on the way in, and an all-zero direction is skipped, so a mostly-empty block costs only the columns that carry signal. The sparsity that actually pays here is the model’s, not v’s: every pass is O(tape ops), never O(n²), whichever way v arrives. On a model with a tridiagonal Hessian — the usual IPM shape — that is the whole game.

The block form is not just a loop: the forward sweep depends only on x, so k directions share one sweep per tape where k separate calls would repeat it. Only the forward-tangent and reverse-over-tangent passes are per-direction.

The result is always dense, including for sparse input. ∇²L · v is dense in general even when both ∇²L and v are sparse, so a sparse return type would advertise an economy the product does not have. When you want the sparse Hessian itself, hessian_structure() + hessian(x) give it directly as a COO lower triangle:

hr, hc = p.hessian_structure()
lower = sp.coo_matrix((p.hessian(x), (hr, hc)), shape=(p.n, p.n)).tocsr()
H = lower + lower.T - sp.diags(lower.diagonal())    # full symmetric matrix

NaN and Inf do not spread through structural zeros. AD never multiplies by an entry that is not in the tape, so a NaN in one component of v stays confined to the variables actually coupled to it. A dense H @ v computes 0 * nan and smears NaN across every row. On a block-diagonal Hessian with v = [nan, 0, 1, 0], the dense product is [nan nan nan nan] while the HVP is [nan nan 2.42 3.08]. Arguably the better semantics, but it does mean the HVP is not bit-equivalent to a dense product on non-finite input.

Sharing one NlProblem across threads

An NlProblem may be built on one thread and evaluated — or garbage collected — on any other. Threaded hosts (a branch-and-bound worker pool, a ThreadPoolExecutor) can hold one shared evaluator rather than one tape per worker:

p = pounce.read_nl("model.nl")

with ThreadPoolExecutor(max_workers=8) as pool:
    values = list(pool.map(p.objective, points))   # one tape, N workers

The evaluators do not release the GIL, so concurrent calls serialize rather than overlap — the win is memory (one copy of the tapes) and the absence of thread-affinity ceremony, not parallel throughput. For actual parallelism across instances, use solve_nlp_batch, which releases the GIL and runs the whole batch on a Rayon pool.

DenseLU and SparseLU carry the same guarantee: factor on one thread, back-solve on another.

Solver, QpFactorization and QpSensitivity are the exceptions — their held Ipopt / KKT factorizations are genuinely thread-affine, so each must be used and released on the thread that created it. Keep them in a threading.local (not a dict keyed by thread id: CPython clears a threading.local on the owning thread as it exits, so the object is both built and dropped where it belongs), which is what pounce.jax’s JaxProblem does internally. Using one from another thread raises a PanicException — note that this derives from BaseException, so an except Exception will not catch it.

Batched NLP solving (solve_nlp_batch)

pounce.solve_nlp_batch solves N independent NLPs and returns one (x, info) pair per input, in input order — for parametric sweeps, multi-start, MPC chains, or branch-and-bound node relaxations where each sibling differs only in tightened bounds.

import numpy as np
import pounce

base = pounce.read_nl("model.nl")          # native-Rust evaluators

# One parsed structure, many variations (cheap clones of the AD tapes):
rng = np.random.default_rng(0)
batch = [base.variant(x0=np.asarray(base.x0) + rng.normal(0, 0.01, base.n))
         for _ in range(24)]

results = pounce.solve_nlp_batch(batch, options={"tol": 1e-8})
for x, info in results:
    print(info["status_msg"], info["obj_val"])

NlProblem.variant(x0=, x_l=, x_u=, g_l=, g_u=) builds a sibling instance with per-instance starting point / bounds; everything structural (expression DAG, AD tapes, sparsity, coloring) is shared work that is not redone.

Native vs. callback inputs — the GIL caveat. Both kinds solve in parallel, with different ceilings:

• NlProblem inputs (from read_nl / variant) are native-Rust reverse-mode-AD evaluators. The batch runs on a Rayon thread pool with the GIL fully released; each worker solves its instance end-to-end with an inner-serial factorization (outer-parallel / inner-serial, the same model as solve_qp_batch).
• Callback-based Problem inputs (pass x0s=, one starting point per instance) also run one instance per worker, but every objective / gradient / constraints / jacobian / hessian call re-acquires the GIL. The Python share of the work is therefore serialized: the speedup scales with the Rust/Python work ratio — medium and large problems whose factorizations dominate parallelize well (~4x on 4 cores for an n=800 banded NLP with vectorized NumPy callbacks); tiny problems whose callbacks dominate won’t. Each Problem’s own add_option settings are honored per instance, with options= as a batch-level overlay.

With parallel=False either path solves one instance at a time, letting each factorization parallelize internally — better for a few large instances. For the batch, print_level defaults to 0 (N workers interleaving iteration tables is noise); pass an explicit print_level to override.

Warm-start chaining (MPC / B&B). Feed one batch’s results into the next solve of a nearby batch:

results = pounce.solve_nlp_batch(batch_t)              # cold
results = pounce.solve_nlp_batch(batch_t1, warms=results)  # warm

Each instance is seeded with the previous x and duals, the converged barrier parameter (info["mu"]) is threaded into mu_init, and warm_start_init_point=yes is forced. A warm start changes iteration counts, never solutions (re-solving the 24-instance gaslib sweep warm drops 482 total iterations to 120). A dimension-mismatched warm entry falls back to that instance’s cold start.

Partial multiplier seeds (Problem.solve / Solver.solve). The lagrange=, zl=, zu= arguments take the solver’s internal conventions (+λ with L = f + λᵀg, non-negative bound multipliers). Under warm_start_init_point=yes, a NaN entry means “unseeded”: the warm-start initializer substitutes its own resolved default (bound_mult_init_val for bound multipliers, the warm path’s 0 for equality duals) before its clamps, so a partial seed never turns into a zero bound multiplier on an active bound, which is a contradictory KKT certificate. This contract belongs to the warm-start initializer only: the batch warms= hand-off above and the SQP working_set arrays do not route through it and must not carry NaN.

Identical-sparsity batches (share_structure=True). When every instance shares its KKT sparsity (parametric sweeps, multi-start, B&B siblings), this opt-in keeps each worker’s factorization backend alive across instances so the symbolic analysis (fill-reducing ordering, supernode structure) runs once per worker rather than once per instance. Always correct — a pattern change just triggers a fresh analysis — but pooled solver state means results are within solver tolerance of, not bit-identical to, the default fresh-backend solves. The win scales with how expensive ordering is for your model (small models: negligible; large sparse models: worth measuring).

scipy.optimize-style

import numpy as np
from pounce import minimize

res = minimize(lambda x: (x - 1) @ (x - 1) + 1, x0=np.zeros(5))
print(res.fun, res.x)

minimize is a thin facade over pounce.Problem shaped after scipy.optimize.minimize, so SciPy code ports with few changes — including as a method= callable handed to scipy.optimize.minimize itself. It returns a genuine scipy.optimize.OptimizeResult (res.x, res.fun, res.success, res.status, res.message, res.nit, and the res.nfev / res.njev / res.nhev evaluation counters), with pounce-specific extras under res.info and a back-compat shim so a key absent at the top level falls back to res.info.

Compatibility with scipy.optimize.minimize

minimize(fun, x0, args=(), jac=None, hess=None, bounds=None,
         constraints=None, callback=None, **options)

Conventions that match SciPy (so constraints port directly):

• Inequalities use the SciPy sign convention g(x) ≥ 0; equalities are g(x) = 0. A LinearConstraint(A, lb, ub) becomes lb ≤ A x ≤ ub.
• The result object is a genuine scipy.optimize.OptimizeResult (subset of fields + an info map).

Gaps worth knowing:

• NonlinearConstraint objects are not accepted — pass nonlinear constraints as the dict form {"type": …, "fun": …, "jac": …}. (Bounds and LinearConstraint objects are accepted.)
• A constraint dict’s Jacobian is dense; for large sparse Jacobians use the Problem class directly (a LinearConstraint may carry a sparse A, which is honored).
• options={"maxiter": 100} now works (scipy synonyms are mapped), but the underlying pounce option is still max_iter; an unrecognized key is forwarded verbatim to the backend.

Solver routing in minimize

By default minimize uses the general NLP filter line-search interior-point method and does no structure probing — an expensive fun pays nothing. Opt in with solver_selection="auto" (the same key the CLI uses) and minimize probes the callables: a problem that is provably a linear program or a convex quadratic program is dispatched to the specialized convex interior-point solver (pounce.solve_qp, the HSDE driver), and a provably convex QCQP (convex-quadratic objective and/or constraints) is reformulated to a second-order cone program and dispatched to the conic solver (pounce.solve_socp). Both reach a global optimum in materially fewer iterations; everything else falls through to the NLP solver.

The catch is that minimize only sees opaque callables — it cannot read a .nl expression tree the way the CLI can. So instead of reading the structure it probes it: it evaluates fun/jac/hess at several points, fits a linear/quadratic model, and then validates that model against the true callables at held-out points before trusting it. The two misclassification directions are not symmetric, and the validation gates the dangerous one:

• A convex LP/QP/QCQP mistakenly sent to the NLP solver is merely slower — the filter-IPM still solves it correctly.
• A genuinely nonlinear or nonconvex problem sent to the convex solver would return a silently wrong answer.

So any probe that raises, any model mismatch beyond route_tol, a non-constant Hessian/Jacobian, an indefinite objective Hessian (a nonconvex QP), a quadratic equality, or a quadratic inequality whose feasible set is nonconvex (a non-PSD constraint Hessian) all fall back to the NLP solver. You never get a wrong “optimum” from a misclassification.

Forcing the solver

The solver_selection option (passed in options=) overrides the automatic choice — mirroring the CLI option of the same name: