fix(qp): align fast/cvxopt objective convention, reject dropped equalities, tune PDIPM

src/cbfkit/optimization/quadratic_program/qp_solver_cvxopt.py

@@ -104,8 +104,12 @@ def solve_with_details(
 
     ``init_params`` is accepted for interface compatibility but ignored
     (CVXOPT does not support warm-starting).
+
+    Note: the registry convention is ``min x'Hx + f'x`` while CVXOPT's native
+    form is ``min 1/2 x'Px + q'x``; the ``P = 2H`` conversion happens here so
+    all ``get_solver()`` backends agree (see ``solver_registry`` docstring).
     """
     from cbfkit.optimization.quadratic_program.solver_registry import QpSolution
 
-    primal, success = solve(h_mat, f_vec, g_mat, h_vec, a_mat, b_vec)
+    primal, success = solve(2.0 * h_mat, f_vec, g_mat, h_vec, a_mat, b_vec)
     return QpSolution(primal=primal, status=1 if success else 0, params=None)

src/cbfkit/optimization/quadratic_program/qp_solver_fast.py

@@ -13,6 +13,7 @@
 from jax import Array
 
 from cbfkit.optimization.quadratic_program.qp_solver_pdipm import (
+    DEFAULT_MAX_ITER,
     PdipmState,
     solve_qp_pdipm,
 )
@@ -24,14 +25,18 @@ def solve_qp_fast(
     G: Array,
     h: Array,
     warm_start: Optional[Array] = None,
-    max_iter: int = 25,
+    max_iter: Optional[int] = None,
     tol: float = 1e-6,
 ) -> Tuple[Array, int, Array]:
     """Backward-compatible alias for solve_qp_pdipm.
 
     Returns ``(x, status, dual)``. For full state including primal slacks
     (needed for richer warm-starting), call ``solve_qp_pdipm`` directly.
 
+    ``max_iter=None`` defers to ``qp_solver_pdipm.DEFAULT_MAX_ITER`` rather
+    than pinning its own copy of the budget, so tuning the solver default is
+    not silently masked here.
+
     Note: the legacy interface accepted an Array as ``warm_start`` (the dual
     only). To preserve that signature, we wrap it in a minimal ``PdipmState``
     with slacks defaulting to ones, which the PDIPM init will then clamp
@@ -55,7 +60,7 @@ def solve_qp_fast(
         G,
         h,
         warm_start=warm,
-        max_iter=max_iter,
+        max_iter=DEFAULT_MAX_ITER if max_iter is None else max_iter,
         tol=tol,
     )
     return x, status, state.dual

src/cbfkit/optimization/quadratic_program/qp_solver_pdipm.py

@@ -46,42 +46,79 @@ def _step_to_boundary(s: Array, ds: Array) -> Array:
 _EPS_PD = 1e-10  # PD regularization on H before Cholesky
 
 
-def _solve_newton_reduced(P: Array, G: Array, s: Array, lam: Array, rhs: Array) -> Array:
-    """Solve (P + G^T diag(lam/s) G + eps*I) dx = rhs via Cholesky.
+def _factor_newton_reduced(P: Array, G: Array, s: Array, lam: Array) -> Array:
+    """Cholesky factor of (P + G^T diag(lam/s) G + eps*I).
 
     The eps*I term guards against rounding-induced non-PD near the optimum.
+
+    The predictor and corrector of a Mehrotra iteration share this matrix and
+    differ only in their right-hand side, so it is factored once per iteration
+    and reused via ``_solve_with_factor``.
     """
     n = P.shape[0]
     D = lam / s
     H = P + G.T @ (D[:, None] * G) + _EPS_PD * jnp.eye(n)
-    L = jnp.linalg.cholesky(H)
-    # dx = H^-1 rhs via two triangular solves
+    return jnp.linalg.cholesky(H)
+
+
+def _solve_with_factor(L: Array, rhs: Array) -> Array:
+    """Solve H dx = rhs given the Cholesky factor L of H, via two triangular solves."""
     y = jax.scipy.linalg.solve_triangular(L, rhs, lower=True)
-    dx = jax.scipy.linalg.solve_triangular(L.T, y, lower=False)
-    return dx
+    return jax.scipy.linalg.solve_triangular(L.T, y, lower=False)
 
 
 _STEP_SAFETY = 0.995  # fraction of step-to-boundary; keeps strict interior
 
-
-def _pdipm_iteration(
+# Default iteration budget. Because the fori_loop always runs to completion,
+# this is the per-solve cost, so it is chosen as the smallest value that is
+# still exact rather than a loose upper bound.
+#
+# Calibrated on two 250-QP suites (n=5, m=24), comparing against max_iter=60:
+#   - benign QPs (well-conditioned P, slack interior): converged by iteration 8.
+#   - hard QPs (cond(P) ~ 1e6 from slack penalty weighting, optimum on the
+#     barrier boundary with many near-active constraints) needed far more:
+#       max_iter=10 -> 123/250 QPs off by >1e-8;  max_iter=12 -> 30/250;
+#       max_iter=14 -> 1/250;   max_iter=16 -> 0/250, all status==1.
+# 16 is the first budget that is exact on both suites. Benign problems are
+# unaffected in accuracy since they freeze on convergence.
+DEFAULT_MAX_ITER = 16
+
+
+def _residuals(
     P: Array,
     q: Array,
     G: Array,
     h: Array,
     x: Array,
     s: Array,
     lam: Array,
+) -> Tuple[Array, Array]:
+    """Dual and primal residuals (r_d, r_p).
+
+    Computed once per iteration and shared by the freeze-on-converge test and
+    the Newton step, which would otherwise each rebuild them.
+    """
+    r_d = P @ x + G.T @ lam + q
+    r_p = G @ x + s - h
+    return r_d, r_p
+
+
+def _pdipm_iteration(
+    P: Array,
+    G: Array,
+    x: Array,
+    s: Array,
+    lam: Array,
+    r_d: Array,
+    r_p: Array,
 ) -> Tuple[Array, Array, Array]:
-    """One Mehrotra predictor-corrector PDIPM step.
+    """One Mehrotra predictor-corrector PDIPM step, given precomputed residuals.
 
     Returns updated (x, s, lam) all strictly positive in s, lam.
     """
     m = G.shape[0]
 
-    # --- Residuals (predictor: r_c with mu = 0) ---
-    r_d = P @ x + G.T @ lam + q
-    r_p = G @ x + s - h
+    # Predictor complementarity residual (r_c with mu = 0).
     r_c_aff = s * lam
 
     # --- Predictor (affine) step ---
@@ -90,7 +127,9 @@ def _pdipm_iteration(
     #                            = -(r_d + G^T ((lam*r_p - s*lam) / s))
     # which simplifies via r_c_aff = s*lam to -(r_d + G^T ((lam*r_p - r_c_aff)/s)).
     rhs_aff = -(r_d + G.T @ ((lam * r_p - r_c_aff) / s))
-    dx_aff = _solve_newton_reduced(P, G, s, lam, rhs_aff)
+    # Predictor and corrector share this factorization; see _factor_newton_reduced.
+    L = _factor_newton_reduced(P, G, s, lam)
+    dx_aff = _solve_with_factor(L, rhs_aff)
     ds_aff = -r_p - G @ dx_aff
     dlam_aff = -(r_c_aff + lam * ds_aff) / s
 
@@ -110,7 +149,7 @@ def _pdipm_iteration(
     # --- Corrector ---
     r_c = s * lam + ds_aff * dlam_aff - sigma * mu
     rhs = -(r_d + G.T @ ((lam * r_p - r_c) / s))
-    dx = _solve_newton_reduced(P, G, s, lam, rhs)
+    dx = _solve_with_factor(L, rhs)
     ds = -r_p - G @ dx
     dlam = -(r_c + lam * ds) / s
 
@@ -125,21 +164,11 @@ def _pdipm_iteration(
     return x_new, s_new, lam_new
 
 
-def _combined_residual(
-    P: Array,
-    q: Array,
-    G: Array,
-    h: Array,
-    x: Array,
-    s: Array,
-    lam: Array,
-) -> Array:
+def _combined_residual(r_d: Array, r_p: Array, s: Array, lam: Array) -> Array:
     """Scalar residual used both for in-loop freeze-on-converge and for the
     post-loop status check: ||r_d||_inf + ||r_p||_inf + mu.
     """
-    r_d = P @ x + G.T @ lam + q
-    r_p = G @ x + s - h
-    mu = jnp.sum(s * lam) / G.shape[0]
+    mu = jnp.sum(s * lam) / s.shape[0]
     return jnp.max(jnp.abs(r_d)) + jnp.max(jnp.abs(r_p)) + mu
 
 
@@ -150,7 +179,7 @@ def solve_qp_pdipm(
     G: Array,
     h: Array,
     warm_start: Optional[PdipmState] = None,
-    max_iter: int = 25,
+    max_iter: int = DEFAULT_MAX_ITER,
     tol: float = 1e-6,
 ) -> Tuple[Array, Array, PdipmState]:
     """Mehrotra predictor-corrector PDIPM for min 0.5 x^T P x + q^T x s.t. G x <= h.
@@ -161,7 +190,14 @@ def solve_qp_pdipm(
     fixed for the remaining iterations rather than stepping further. This
     prevents a degenerate (s ≈ 0, lam ≈ 0) post-convergence state from
     propagating NaN through ``lam/s`` divisions in subsequent iterations.
-    The total flop count is still ``max_iter`` Newton solves per call.
+    The total flop count is still ``max_iter`` Newton solves per call, which
+    is why ``max_iter`` defaults to the tight ``DEFAULT_MAX_ITER`` rather than
+    a loose upper bound; raise it for QPs harder than the CBF-QP shapes it was
+    calibrated on (see that constant for the measurements).
+
+    ``fori_loop`` is deliberate and load-bearing: ``lax.while_loop`` has no
+    reverse-mode differentiation rule, and downstream differentiable-CBF-QP
+    work takes ``jax.grad`` through this solver.
 
     Returns:
         (x, status, state) where status is:
@@ -193,8 +229,10 @@ def solve_qp_pdipm(
     # --- Outer loop (fixed iterations, freeze-on-converge) ---
     def body(_, carry):
         x_, s_, lam_ = carry
-        already_converged = _combined_residual(P, q, G, h, x_, s_, lam_) < tol
-        x_new, s_new, lam_new = _pdipm_iteration(P, q, G, h, x_, s_, lam_)
+        # One residual evaluation feeds both the convergence test and the step.
+        r_d, r_p = _residuals(P, q, G, h, x_, s_, lam_)
+        already_converged = _combined_residual(r_d, r_p, s_, lam_) < tol
+        x_new, s_new, lam_new = _pdipm_iteration(P, G, x_, s_, lam_, r_d, r_p)
         # If already converged, keep current (good) state rather than stepping.
         # Both branches of jnp.where are traced; the step is cheap relative to
         # the cost of NaN propagation if we let lam/s explode after convergence.
@@ -206,7 +244,8 @@ def body(_, carry):
     x_final, s_final, lam_final = lax.fori_loop(0, max_iter, body, (x0, s0, lam0))
 
     # --- Status: solved if final residual norm below tol ---
-    res_norm = _combined_residual(P, q, G, h, x_final, s_final, lam_final)
+    r_d_final, r_p_final = _residuals(P, q, G, h, x_final, s_final, lam_final)
+    res_norm = _combined_residual(r_d_final, r_p_final, s_final, lam_final)
     status = jnp.where(res_norm < tol, jnp.int32(1), jnp.int32(2))
 
     state = PdipmState(x=x_final, s=s_final, dual=lam_final, iter_num=max_iter)

src/cbfkit/optimization/quadratic_program/solver_registry.py

@@ -3,6 +3,19 @@
 Provides a common ``QpSolution`` return type and factory functions that wrap
 each backend (jaxopt, cvxopt, casadi) behind a single callable signature.
 
+**Objective convention.** Every solver returned by :func:`get_solver` solves
+
+.. math::
+
+    \\min_x \\; x^T H x + f^T x \\quad \\text{s.t.} \\quad G x \\le h \\;(, A x = b)
+
+This matches the CBF-CLF-QP generator, which passes ``f = -2 H u_nom`` so the
+unconstrained minimizer is ``u_nom``.  Backends whose native form is
+``min 1/2 x^T P x + q^T x`` (OSQP, CVXOPT, PDIPM) are adapted at this boundary
+(``P = 2 H`` or equivalently ``q = f/2``); the raw modules
+(``qp_solver_pdipm.solve_qp_pdipm``, ``qp_solver_cvxopt.solve``) keep their
+native halved convention.
+
 Usage::
 
     from cbfkit.optimization.quadratic_program import get_solver
@@ -168,7 +181,7 @@ def casadi_solver() -> QpSolverCallable:
 # ---------------------------------------------------------------------------
 
 
-def fast_solver(max_iter: int = 25, tol: float = 1e-6) -> QpSolverCallable:
+def fast_solver(max_iter: Optional[int] = None, tol: float = 1e-6) -> QpSolverCallable:
     """Fast PDIPM solver for small CBF-CLF problems.
 
     Mehrotra predictor-corrector primal-dual interior-point method. Designed
@@ -181,15 +194,24 @@ def fast_solver(max_iter: int = 25, tol: float = 1e-6) -> QpSolverCallable:
     (see ``benchmarks/qp_solver_comparison.py``). JIT-compatible and
     warm-startable across consecutive control steps.
 
+    Inequality constraints only: passing ``a_mat``/``b_vec`` raises
+    ``NotImplementedError`` rather than dropping them.
+
     Args:
-        max_iter: Maximum PDIPM iterations (default 25; ~10-15 typically suffice).
+        max_iter: Maximum PDIPM iterations. ``None`` (default) defers to
+            ``qp_solver_pdipm.DEFAULT_MAX_ITER`` so the solver's calibrated
+            budget is not silently pinned here.
         tol: Combined primal/dual/complementarity residual tolerance.
     """
     from cbfkit.optimization.quadratic_program.qp_solver_pdipm import (
+        DEFAULT_MAX_ITER,
         PdipmState,
         solve_qp_pdipm,
     )
 
+    if max_iter is None:
+        max_iter = DEFAULT_MAX_ITER
+
     def solve_with_details(
         h_mat: Array,
         f_vec: Array,
@@ -199,8 +221,19 @@ def solve_with_details(
         b_vec: Optional[Array] = None,
         init_params: Any = None,
     ) -> QpSolution:
+        if a_mat is not None or b_vec is not None:
+            raise NotImplementedError(
+                "The 'fast' (PDIPM) backend solves inequality-constrained QPs only "
+                "(min x'Hx + f'x s.t. Gx <= h), but was given equality constraints "
+                "via a_mat/b_vec. Silently dropping them would return a solution "
+                "that violates Ax = b — e.g. an MPC trajectory ignoring its own "
+                "dynamics. Use get_solver('jaxopt') or get_solver('casadi'), which "
+                "support equality constraints."
+            )
+
         if g_mat is None or h_vec is None:
-            x = jnp.linalg.solve(h_mat, -f_vec)
+            # Registry convention min x'Hx + f'x  =>  2H x* = -f.
+            x = jnp.linalg.solve(2.0 * h_mat, -f_vec)
             return QpSolution(primal=x, status=1, params=None)
 
         # Extract warm-start state from previous QpSolution.params
@@ -213,8 +246,11 @@ def solve_with_details(
             elif isinstance(init_params, PdipmState):
                 warm = init_params
 
+        # Convention adapter: solve_qp_pdipm natively solves
+        # min 1/2 x'Px + q'x; the registry convention is min x'Hx + f'x,
+        # so pass P = 2H (the jaxopt wrapper adapts via q = f/2 instead).
         sol, status, state = solve_qp_pdipm(
-            h_mat,
+            2.0 * h_mat,
             f_vec,
             g_mat,
             h_vec,

tests/test_optimization/test_fast_qp.py

@@ -1,7 +1,5 @@
 """Tests for the fast small-QP solver and its integration with CBF-CLF-QP."""
 
-import pytest
-import jax
 import jax.numpy as jnp
 from jax import random
 
@@ -110,15 +108,21 @@ def test_get_fast_solver(self):
         assert solver.solver_name == "fast"
 
     def test_fast_solver_via_registry(self):
+        """Registry convention: min x'Hx + f'x (matches jaxopt/casadi backends).
+
+        With H=I, f=-2*[1,2] the minimizer is [1,2] — NOT [2,4], which is
+        what the pre-fix wrapper returned by feeding (H, f) straight into
+        the native ``min 1/2 x'Px + q'x`` PDIPM form.
+        """
         solver = get_solver("fast")
         P = jnp.eye(2)
         q = jnp.array([-2.0, -4.0])
         G = jnp.array([[1, 0], [-1, 0], [0, 1], [0, -1.0]])
         h = jnp.array(
             [3.0, 3.0, 5.0, 5.0]
-        )  # box [-3,3] x [-5,5] — unconstrained opt [2,4] is feasible
+        )  # box [-3,3] x [-5,5] — unconstrained opt [1,2] is feasible
         sol = solver(P, q, G, h)
-        assert jnp.allclose(sol.primal, jnp.array([2.0, 4.0]), atol=1e-4)
+        assert jnp.allclose(sol.primal, jnp.array([1.0, 2.0]), atol=1e-4)
         assert int(sol.status) == 1
 
     def test_fast_solver_warm_start_via_registry(self):