herculesabqp 0.1.2

# HerculesABQP

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

```text
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:

```bash
cargo bench --bench solver_bench
```

Build the Python extension in a virtual environment:

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

Install Python test dependencies and run the binding tests:

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

## Example Use

Solve one box-constrained QP directly:

```rust
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:

```rust
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/](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/](examples/).
- Core references are collected in [references.bib](references.bib).