QP Solver Backend: HiGHS

QP Solver Backend: HiGHS#

This notebook demonstrates the HiGHS-based QP/MIQP solver backend in discopt. When HiGHS (highspy) is available, discopt automatically dispatches quadratic programs to HiGHS instead of the general-purpose NLP interior-point method, yielding significant performance improvements.

HiGHS [Huangfu and Hall, 2018] is a high-performance open-source solver for linear programming, mixed-integer programming, and quadratic programming. The discopt integration extracts the quadratic objective (the Hessian matrix \(Q\)) and linear constraints from the expression DAG and passes them directly to HiGHS via its Python API.

By the end of this notebook you will understand how discopt dispatches QP and MIQP problems to HiGHS, see the performance difference compared to the NLP path, and formulate a mixed-integer quadratic program solved natively by HiGHS.

import os

os.environ["JAX_PLATFORMS"] = "cpu"
os.environ["JAX_ENABLE_X64"] = "1"

import time

import discopt.modeling as dm
import numpy as np

1. QP Formulation and Solve via HiGHS#

We begin with a simple convex QP: minimize a quadratic objective subject to a linear constraint. When discopt detects that the objective is quadratic and all constraints are linear (with all continuous variables), it classifies the problem as a QP and dispatches to HiGHS 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.0291s

# Verify the problem is classified as QP and dispatched to HiGHS
from discopt._jax.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: HiGHS vs NLP Path#

To demonstrate the performance advantage of the HiGHS QP backend, we solve a portfolio optimization QP using both backends and compare wall-clock times. The HiGHS backend uses a specialized active-set QP solver, while the NLP path uses a general-purpose interior-point method with JAX autodiff for gradient and Hessian evaluation.

from discopt.solver import _solve_qp_jax


def build_portfolio_qp(n_assets, target_return=0.08, seed=42):
    """Build a Markowitz portfolio QP with n_assets."""
    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)

    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),
        )
    )
    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


# Solve via HiGHS (default dispatch for QP)
m_highs = build_portfolio_qp(20)
t0 = time.perf_counter()
result_highs = m_highs.solve()
time_highs = time.perf_counter() - t0

# Solve via JAX IPM (bypass HiGHS by using the internal JAX path)
m_jax = build_portfolio_qp(20)
t0 = time.perf_counter()
result_jax = _solve_qp_jax(m_jax, time.perf_counter())
time_jax = time.perf_counter() - t0

header = f"{'Backend':<12s} {'Time (s)':>10s} {'Objective':>12s} {'Status':>10s}"
print(header)
print("-" * 48)
print(
    f"{'HiGHS':<12s} {time_highs:10.4f} {result_highs.objective:12.6f} {result_highs.status:>10s}"
)
print(f"{'JAX IPM':<12s} {time_jax:10.4f} {result_jax.objective:12.6f} {result_jax.status:>10s}")

if time_jax > 0:
    speedup = time_jax / time_highs
    print(f"\nHiGHS is {speedup:.1f}x faster than JAX IPM.")

Backend        Time (s)    Objective     Status
------------------------------------------------
HiGHS            0.0300     0.000987    optimal
JAX IPM          0.1756     0.000987    optimal

HiGHS is 5.8x faster than JAX IPM.

The HiGHS backend avoids JAX tracing overhead and uses a compiled C++ solver optimized for QP structure. For small to medium QPs, this overhead reduction is the dominant factor. For larger problems, the algorithmic advantage of HiGHS’s active-set method (which exploits sparsity in the KKT system) becomes increasingly important.

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 an MIQP, it dispatches to HiGHS which solves it via branch-and-bound with QP relaxations at each node.

# 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: 1.5932s

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 solver dispatch logic in discopt follows this sequence for QP/MIQP problems.

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

QP (continuous, quadratic objective, linear constraints): HiGHS QP solver via passHessian(). Falls back to the JAX QP IPM if HiGHS is unavailable.
MIQP (integer/binary variables, quadratic objective, linear constraints): HiGHS MIQP solver (branch-and-bound with QP relaxations). Falls back to discopt’s own B&B with JAX QP relaxations if HiGHS is unavailable.

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#

The HiGHS QP backend provides a significant speedup for quadratic and mixed-integer quadratic programs by leveraging a compiled, specialized solver instead of the general NLP interior-point path. The dispatch is transparent: users simply call model.solve() and discopt selects the best available backend automatically. When HiGHS is not installed, the solver gracefully falls back to the JAX-based QP IPM, so all code remains portable.