QP Solver Backend: POUNCE

QP Solver Backend: POUNCE#

This notebook demonstrates the specialized QP/MIQP dispatch in discopt. When the problem classifier recognizes a quadratic objective over linear constraints, discopt routes the model to a dedicated quadratic engine instead of the general-purpose spatial branch-and-bound path, yielding a large speedup.

The engine is POUNCE, a pure-Rust port of the Ipopt [Wächter and Biegler, 2006] primal-dual interior-point algorithm. This is a deliberate design choice: the QP path has no HiGHS [Huangfu and Hall, 2018] backend (issue #359) – HiGHS solves only pure LP and MILP models, with discopt-verified certificates – in favor of a self-contained Rust core, and the older pure-JAX QP interior-point method was retired along with it. A continuous QP is therefore solved by POUNCE and nothing else — there is no second engine and no fallback.

That last point is a soundness decision, not an oversight. _solve_qp re-checks the point POUNCE returns for primal feasibility and for a stationary KKT residual, and refuses a point that fails either. The old JAX rescue path re-solved the same QP and reported status="optimal", gap=0 on nothing but its own internal convergence flag — so the one situation that could reach it was exactly the situation where a verified engine’s answer had just been thrown out. Under the rule that a global solver’s product is its certificate, an unverified rescue is worse than no rescue.

By the end of this notebook you will understand how discopt classifies and dispatches QP and MIQP problems, see the speedup relative to the general nonlinear path, and formulate a mixed-integer quadratic program.

import time

import discopt.modeling as dm
import numpy as np

1. QP Formulation and Solve#

We begin with a simple convex QP: minimize a quadratic objective subject to a linear constraint. When discopt detects that the objective is at most quadratic and all constraints are linear (with all continuous variables), it classifies the problem as a QP and dispatches to POUNCE automatically.

# Simple QP: min 0.5*x^2 + y^2  s.t. x + y >= 1, x,y >= 0
m = dm.Model("simple_qp")
x = m.continuous("x", lb=0, ub=10)
y = m.continuous("y", lb=0, ub=10)
m.minimize(0.5 * x**2 + y**2)
m.subject_to(x + y >= 1)

result = m.solve()
print(f"Status:    {result.status}")
print(f"Objective: {result.objective:.6f}")
print(f"x = {float(result.value(x)):.6f}")
print(f"y = {float(result.value(y)):.6f}")
print(f"Wall time: {result.wall_time:.4f}s")

Status:    optimal
Objective: 0.333333
x = 0.666667
y = 0.333333
Wall time: 0.0303s

# Verify the problem is classified as QP and dispatched to the quadratic path
from discopt._relax.problem_classifier import ProblemClass, classify_problem

pc = classify_problem(m)
print(f"Problem class: {pc}")
assert pc == ProblemClass.QP, "Expected QP classification"

Problem class: ProblemClass.QP

2. Performance Comparison: QP Path vs General Nonlinear Path#

To show what the specialized dispatch buys, we solve the same Markowitz portfolio problem two ways. The first is the ordinary solve(), which the classifier sends down the QP path. The second adds a 0 * exp(w[0]) term to the objective — numerically a no-op, but enough to push the model out of the quadratic class, so it goes through the general nonlinear path with its DAG evaluation and relaxation machinery.

Both should reach the same objective; the difference is entirely in how much work is done to get there.

The head-to-head runs at a deliberately small size, five assets. Spatial branch-and-bound on a dense covariance objective degrades quickly with dimension, and past roughly ten assets it stops certifying optimality within any reasonable budget — it returns a good incumbent it cannot prove. Comparing the two paths where both still close the bound makes the speedup a like-for-like measurement rather than a race against a time limit. The cell ends by running the QP path alone at twenty assets, where the contrast is starkest: still a root solve.

def build_portfolio_qp(n_assets, target_return=0.08, seed=42, force_nonlinear=False):
    """Build a Markowitz portfolio QP with n_assets.

    ``force_nonlinear`` adds a numerically-zero ``exp`` term, which keeps the
    optimum identical but takes the model out of the quadratic class so the
    classifier routes it down the general nonlinear path.
    """
    rng = np.random.RandomState(seed)
    mu = rng.uniform(0.02, 0.15, n_assets)
    L = rng.randn(n_assets, n_assets) * 0.05
    Sigma = L @ L.T + 0.01 * np.eye(n_assets)

    m = dm.Model(f"portfolio_{n_assets}")
    w = m.continuous("w", shape=(n_assets,), lb=0, ub=1)

    obj = dm.sum(
        lambda i: dm.sum(
            lambda j: Sigma[i, j] * w[i] * w[j],
            over=range(n_assets),
        ),
        over=range(n_assets),
    )
    if force_nonlinear:
        obj = obj + 0 * dm.exp(w[0])
    m.minimize(obj)
    m.subject_to(dm.sum(w) == 1, name="budget")
    m.subject_to(
        dm.sum(lambda i: mu[i] * w[i], over=range(n_assets)) >= target_return,
        name="min_return",
    )
    return m


# Compare at a size the *general* path can still prove optimal on. Spatial
# branch-and-bound does not scale on a dense covariance objective: measured on
# this branch, n=8 needs ~32 s and 661 nodes, n=10 uses the full 120 s budget,
# and n=20 only reaches `feasible`. The QP path solves every one of them at the
# root. Comparing where both paths *certify* optimality keeps the speedup below
# an honest one rather than a race against a time limit.
N_CMP = 5

# Specialized QP path (POUNCE), selected automatically by the classifier
m_qp = build_portfolio_qp(N_CMP)
print(f"classified as {classify_problem(m_qp)}")
t0 = time.perf_counter()
result_qp = m_qp.solve()
time_qp = time.perf_counter() - t0

# General nonlinear path. The time limit is a guard, not the expected exit:
# at this size the path converges well inside it.
m_nlp = build_portfolio_qp(N_CMP, force_nonlinear=True)
print(f"classified as {classify_problem(m_nlp)}")
t0 = time.perf_counter()
result_nlp = m_nlp.solve(time_limit=120)
time_nlp = time.perf_counter() - t0

header = f"{'Path':<16s} {'Time (s)':>10s} {'Objective':>12s} {'Status':>10s} {'Nodes':>8s}"
print()
print(header)
print("-" * 62)
print(
    f"{'QP (POUNCE)':<16s} {time_qp:10.4f} {result_qp.objective:12.6f} "
    f"{result_qp.status:>10s} {result_qp.node_count:8d}"
)
print(
    f"{'general NLP':<16s} {time_nlp:10.4f} {result_nlp.objective:12.6f} "
    f"{result_nlp.status:>10s} {result_nlp.node_count:8d}"
)

assert result_qp.status == "optimal", "QP path did not certify optimality"
assert result_nlp.status == "optimal", "NLP path did not certify optimality"
assert abs(result_qp.objective - result_nlp.objective) < 1e-6, "paths disagree on the optimum"
if time_qp > 0:
    print(f"\nThe QP path is {time_nlp / time_qp:.1f}x faster, at the same objective.")

# The gap widens with size. At n=20 the QP path is still a root solve, while
# the general path above could not close the bound within its budget at all.
m_qp20 = build_portfolio_qp(20)
t0 = time.perf_counter()
result_qp20 = m_qp20.solve()
print(
    f"n=20 on the QP path: {result_qp20.status}, {result_qp20.node_count} nodes, "
    f"{time.perf_counter() - t0:.3f} s, objective {result_qp20.objective:.6f}"
)

classified as ProblemClass.QP
classified as ProblemClass.NLP

Path               Time (s)    Objective     Status    Nodes
--------------------------------------------------------------
QP (POUNCE)          0.0124     0.004388    optimal        0
general NLP          1.4439     0.004388    optimal       83

The QP path is 116.3x faster, at the same objective.
n=20 on the QP path: optimal, 0 nodes, 0.095 s, objective 0.000987

The QP path wins because it never builds the machinery the general path needs. extract_qp_data reads the \(Q\) matrix, linear cost \(c\), and constraint matrix \(A\) straight off the expression DAG — no automatic differentiation, no relaxation, no tree — and hands a matrix problem to POUNCE. The general path must compile the DAG, evaluate derivatives, and construct convex relaxations, all to solve a problem whose structure was already known in closed form.

3. MIQP Example: Cardinality-Constrained Portfolio#

Mixed-integer quadratic programs (MIQPs) extend QPs with integer or binary variables. A common application is portfolio optimization with a cardinality constraint: limit the number of assets held to at most \(k\). This requires binary indicator variables \(z_i \in \{0, 1\}\) where \(z_i = 1\) means asset \(i\) is included in the portfolio.

When discopt detects a convex MIQP it dispatches to its own branch-and-bound over QP relaxations [Land and Doig, 1960], with POUNCE solving each node. Convexity is checked first and the check is load-bearing: _solve_miqp_bb assumes each node relaxation is solved to global optimality, so on an indefinite objective it would return a local stationary point and certify it as global. An indefinite or otherwise non-convex MIQP therefore falls through to the sound spatial McCormick branch-and-bound instead.

Note

The cell below prints a line from the outer-approximation route:

OA: nonlinear objective is not convex in the optimization sense; disabling master lower-bound updates and skipping objective OA cuts

The objective is convex — Sigma is built as L @ L.T + 0.005 I, positive definite by construction. OA’s own convexity classifier false-negatives on a quadratic written as a nested dm.sum, which is #1352. It is a bound-quality defect, not a soundness one: declaring the objective non-convex makes OA decline to update its master lower bound and skip its objective cuts, which is the conservative direction. The answer below is correct and certified; the route just pays for it. The line is left in the output rather than filtered, so the cost is on the record.

# Cardinality-constrained portfolio: hold at most k=3 of 8 assets
n_assets = 8
k_max = 3
np.random.seed(42)

mu = np.array([0.12, 0.10, 0.07, 0.03, 0.15, 0.09, 0.11, 0.06])
L = np.random.randn(n_assets, n_assets) * 0.04
Sigma = L @ L.T + 0.005 * np.eye(n_assets)
target_return = 0.08

m = dm.Model("cardinality_portfolio")
w = m.continuous("w", shape=(n_assets,), lb=0, ub=1)
z = m.binary("z", shape=(n_assets,))  # z[i]=1 if asset i is held

# Minimize portfolio variance
m.minimize(
    dm.sum(
        lambda i: dm.sum(
            lambda j: Sigma[i, j] * w[i] * w[j],
            over=range(n_assets),
        ),
        over=range(n_assets),
    )
)

# Budget constraint
m.subject_to(dm.sum(w) == 1, name="budget")

# Minimum return
m.subject_to(
    dm.sum(lambda i: mu[i] * w[i], over=range(n_assets)) >= target_return,
    name="min_return",
)

# Linking: w[i] <= z[i] (can only invest if asset is selected)
for i in range(n_assets):
    m.subject_to(w[i] <= z[i], name=f"link_{i}")

# Cardinality: at most k assets
m.subject_to(dm.sum(z) <= k_max, name="cardinality")

# Verify MIQP classification
pc = classify_problem(m)
print(f"Problem class: {pc}")

result = m.solve()
print(f"\nStatus:    {result.status}")
print(f"Variance:  {result.objective:.6f}")
print(f"Wall time: {result.wall_time:.4f}s")

w_vals = result.value(w)
z_vals = result.value(z)
print(f"\nSelected assets (z=1): {np.where(z_vals > 0.5)[0]}")
print(f"Weights: {np.round(w_vals, 4)}")
print(f"Return:  {np.dot(mu, w_vals):.4f}")
print(f"Assets held: {int(np.sum(z_vals > 0.5))} (max {k_max})")

Problem class: ProblemClass.MIQP

Status:    optimal
Variance:  0.002238
Wall time: 0.0628s

Selected assets (z=1): [0 2 5]
Weights: [0.3107 0.     0.3579 0.     0.     0.3315 0.     0.    ]
Return:  0.0922
Assets held: 3 (max 3)

4. How the Dispatch Works#

The problem classifier walks the Rust expression IR to check whether the objective is at most quadratic and all constraints are linear. Based on that classification, the presence of integer variables, and a convexity check, the solver selects a path:

QP (continuous, quadratic objective, linear constraints): solved by POUNCE via _solve_qp, with the returned point re-verified for primal feasibility and KKT stationarity. There is no fallback engine — a point that fails either check is refused rather than replaced with an unverified one. A QP whose objective is not known convex is forced onto the spatial path instead.
MIQP (integer/binary variables, quadratic objective, linear constraints): if the objective is known convex, discopt’s own branch-and-bound with QP relaxations at each node; otherwise the spatial McCormick branch-and-bound.

The algebraic coefficient extraction (extract_qp_data) walks the expression DAG to directly read off the \(Q\) matrix, linear cost \(c\), and constraint matrix \(A\) without any automatic differentiation. This makes the data extraction essentially free compared to the solve time.

5. Summary#

Recognizing quadratic structure and dispatching to a specialized engine gives a large speedup over the general nonlinear path, because the \(Q\), \(c\), and \(A\) data can be read algebraically off the DAG rather than rediscovered by differentiation and relaxation. The dispatch is transparent: users call model.solve() and discopt selects the path.

Two things worth remembering from this notebook. First, the quadratic path is guarded by convexity, not just by shape — a quadratic objective that is not known convex is deliberately routed to the slower spatial solver, because the fast path’s certificate would not be valid. Second, there is no fallback engine behind POUNCE, by design: a rescue path that answers without checking feasibility and stationarity produces a certificate the solver cannot stand behind, which is worse than an honest failure.