HerculesABQP

HerculesABQP is a solver for convex box-constrained quadratic programs of the form

minimize    1/2 x^T Q x + c^T x
subject to  l <= x <= u

It is especially useful when the same quadratic model is solved many times with changing bounds and warm starts.

How It Works

The solver does two things:

1. It takes accelerated projected gradient steps on the box-constrained QP, with adaptive restart and optional diagonal scaling.
2. Once the first-order phase is close, (optionally) it applies a small number of active-set polishing steps in the style of ProxQP to clean up the final solution. This is based on the main iteration idea from proxqp.

Features

• accelerated projected gradient descent with adaptive gradient restart
• Gershgorin or auto Lipschitz estimation
• optional Hessian-diagonal scaling
• active-set polishing with reduced solves
• dense and sparse matrix support
• warm starts and a PreparedSolver path for repeated solves
• explicit and implicit operator interfaces
• benchmarks and validation tests

Quick Start

Run the Criterion benchmarks:

cargo bench --bench solver_bench

Build the Python extension in a virtual environment:

maturin develop --manifest-path python/Cargo.toml

Install Python test dependencies and run the binding tests:

python -m pip install "./python[test]"
python -m pytest tests/test_python_bindings.py -q

Example Use

Solve one box-constrained QP directly:

use herculesabqp::matrix::QuadraticMatrix;
use herculesabqp::solver::{solve_box_qp, SolverOptions};
use ndarray::arr2;

fn main() {
    let q = QuadraticMatrix::dense(arr2(&[
        [4.0, 1.0],
        [1.0, 2.0],
    ]));
    let c = vec![-1.0, -0.5];
    let lb = vec![0.0, 0.0];
    let ub = vec![1.0, 1.0];
    let options = SolverOptions::default();

    let result = solve_box_qp(&q, &c, &lb, &ub, &options)
        .expect("failed to solve box-constrained QP");
    println!("x = {:?}", result.x);
    println!("objective = {}", result.objective);
}

Prepare a solver once and reuse it for repeated solves with changing bounds:

use herculesabqp::matrix::QuadraticMatrix;
use herculesabqp::solver::{PreparedSolver, SolverOptions};
use ndarray::arr2;

fn main() {
    let q = QuadraticMatrix::dense(arr2(&[
        [4.0, 1.0],
        [1.0, 2.0],
    ]));
    let c = vec![-1.0, -0.5];
    let options = SolverOptions::default();

    let prepared = PreparedSolver::new(&q, &c, &options)
        .expect("failed to prepare reusable solver state");

    let root = prepared
        .solve_subproblem(&[0.0, 0.0], &[1.0, 1.0], None)
        .expect("failed to solve root QP");

    let child = prepared
        .solve_subproblem(&[1.0, 0.0], &[1.0, 1.0], Some(&root.x))
        .expect("failed to solve child QP");

    println!("root x = {:?}", root.x);
    println!("child x = {:?}", child.x);
}

PreparedSolver is the right entry point when Q and c stay fixed across many solves. It prepares the matrix-side work once, then lets later solves vary only in the bounds and optional warm start.

Python Interface

The optional Python bindings support:

• dense numpy.ndarray inputs for Q
• sparse SciPy matrices or arrays for Q through their tocsr() interface
• dense NumPy vectors for c, lb, ub, and optional x0
• Python-defined implicit operators with n and matvec(x)

The main Python entry points are:

• herculesabqp.solve_box_qp(...) for one-off solves
• herculesabqp.PreparedSolver(...) for repeated solves with shared Q and c
• herculesabqp.solve_box_qp_implicit(...) for Python-defined implicit operators
• herculesabqp.PreparedImplicitSolver(...) for repeated implicit solves

The Python extension lives in the thin wrapper crate under python/, so normal Rust dependencies only build the Rust library API and do not trigger a cdylib link step.

Notes

• The solver assumes a convex QP.
• By default the solver symmetrizes the supplied matrix. Set options.assume_symmetric = true if Q is already symmetric and you want to skip that preprocessing step.
• For repeated solves, prefer the prepared solver APIs.
• For implicit operators, diagonal scaling is only available when the operator exposes diagonal(). If scaling="hessian_diag" is requested without a diagonal, the implicit path falls back to an unscaled solve.
• The implicit path is first-order only and currently skips the final polishing phase.
• Additional examples live in examples/.
• Core references are collected in references.bib.