herculesabqp 0.1.2

A convex box-constrained quadratic programming solver with warm starts and active-set polishing.
Documentation
/// Options controlling first-order iterations, diagnostics, scaling, and polishing.
#[derive(Clone, Debug)]
pub struct SolverOptions {
    /// Optional variable scaling applied before solving.
    pub scaling: ScalingOptions,
    /// Optional warm start in the original variable coordinates.
    pub x0: Option<Vec<f64>>,
    /// Trust that the supplied quadratic matrix is already symmetric.
    ///
    /// When false, the solver defensively replaces `Q` by `0.5 * (Q + Q^T)`
    /// before solving. Set this to true only when the caller already knows the
    /// matrix is symmetric and wants to skip that preprocessing step.
    pub assume_symmetric: bool,
    /// Strategy used to choose the gradient Lipschitz constant.
    pub lipschitz: LipschitzOptions,
    /// Iteration limits and convergence tolerances.
    pub stopping: StoppingOptions,
    /// Optional active-set polish settings.
    pub polish: PolishOptions,
    /// Verbose logging configuration.
    pub logging: LoggingOptions,
}

impl Default for SolverOptions {
    fn default() -> Self {
        Self {
            scaling: ScalingOptions::default(),
            x0: None,
            assume_symmetric: false,
            lipschitz: LipschitzOptions::default(),
            stopping: StoppingOptions::default(),
            polish: PolishOptions::default(),
            logging: LoggingOptions::default(),
        }
    }
}

/// Optional variable scaling applied before solving.
#[derive(Clone, Debug)]
pub struct ScalingOptions {
    /// Solve in the original coordinates or a diagonally scaled variable space.
    pub mode: ScalingMode,
    /// Small positive floor used when building diagonal scaling.
    pub eps: f64,
    /// Lower clamp for each diagonal scaling entry.
    pub min: f64,
    /// Upper clamp for each diagonal scaling entry.
    pub max: f64,
}

impl Default for ScalingOptions {
    fn default() -> Self {
        Self {
            mode: ScalingMode::HessianDiag,
            eps: 1e-12,
            min: 1e-6,
            max: 1e6,
        }
    }
}

/// Optional variable scaling applied before solving.
#[derive(Clone, Debug, PartialEq, Eq)]
pub enum ScalingMode {
    /// Solve in the original coordinates.
    None,
    /// Scale variables using the square root of the Hessian diagonal.
    HessianDiag,
}

/// Strategy for estimating the Lipschitz constant of the gradient.
#[derive(Clone, Debug)]
pub struct LipschitzOptions {
    /// User-supplied Lipschitz constant. When absent, the solver estimates one.
    pub value: Option<f64>,
    /// Strategy used when estimating the Lipschitz constant.
    pub method: LipschitzMethod,
}

impl Default for LipschitzOptions {
    fn default() -> Self {
        Self {
            value: None,
            method: LipschitzMethod::Auto,
        }
    }
}

/// Strategy for estimating the Lipschitz constant of the gradient.
#[derive(Clone, Debug)]
pub enum LipschitzMethod {
    /// Use the Gershgorin row-sum bound directly.
    Gershgorin,
    /// Use the smaller of Gershgorin and a short safeguarded power estimate.
    Auto,
}

impl LipschitzMethod {
    /// Human-readable identifier used in logs and benchmark output.
    pub fn name(&self) -> &'static str {
        match self {
            Self::Gershgorin => "gershgorin",
            Self::Auto => "auto",
        }
    }
}

/// Iteration limits and convergence tolerances for the first-order solve.
#[derive(Clone, Debug)]
pub struct StoppingOptions {
    /// Maximum number of accelerated projected-gradient iterations.
    pub max_iter: usize,
    /// Main solver tolerance for relative gap and KKT stopping tests.
    pub tol: f64,
    /// Compute the affine-minorant dual certificate and relative gap.
    ///
    /// When disabled, the solver skips gap computation and uses the KKT-style
    /// box residual alone for iteration stopping. This is useful on very large
    /// problems where certification is substantially more expensive than the
    /// projected-gradient steps themselves.
    pub dual_certification: bool,
    /// Frequency for recomputing objective and diagnostics.
    pub check_every: usize,
    /// Tolerance used when deciding whether a variable is on a bound.
    pub bound_tol: f64,
}

impl Default for StoppingOptions {
    fn default() -> Self {
        Self {
            max_iter: 5_000,
            tol: 1e-6,
            dual_certification: true,
            check_every: 100,
            bound_tol: 1e-10,
        }
    }
}

/// Optional active-set polish settings applied after the first-order loop.
#[derive(Clone, Debug)]
pub struct PolishOptions {
    /// Enable the active-set polishing phase after first-order iterations.
    pub enabled: bool,
    /// Bound tolerance used during active-set detection in polish.
    pub bound_tol: f64,
    /// Gradient tolerance used during active-set detection in polish.
    pub grad_tol: f64,
    /// Proximal shifts tried during reduced polish solves.
    pub rho_sequence: Vec<f64>,
    /// Maximum number of outer polish passes.
    pub num_passes: usize,
    /// Prefer a dense direct solve for small dense reduced systems.
    pub prefer_direct: bool,
    /// Maximum reduced dimension for the dense direct solve path.
    pub direct_max_n: usize,
    /// Relative tolerance for CG in reduced sparse/dense iterative solves.
    pub cg_rtol: f64,
    /// Maximum number of CG iterations in reduced solves.
    pub cg_maxiter: usize,
    /// Enable diagonal preconditioning in reduced CG solves.
    pub use_preconditioner: bool,
    /// Maximum number of backtracking steps during polish line search.
    pub max_backtracks: usize,
}

impl Default for PolishOptions {
    fn default() -> Self {
        Self {
            enabled: true,
            bound_tol: 1e-7,
            grad_tol: 1e-8,
            rho_sequence: vec![0.0, 1e-10, 1e-8, 1e-6, 1e-4],
            num_passes: 2,
            prefer_direct: true,
            direct_max_n: 1500,
            cg_rtol: 1e-10,
            cg_maxiter: 10_000,
            use_preconditioner: true,
            max_backtracks: 20,
        }
    }
}

/// Verbose logging configuration for iterative solves.
#[derive(Clone, Debug)]
pub struct LoggingOptions {
    /// Print iteration progress to stdout.
    pub verbose: bool,
    /// Print cadence when `verbose` is enabled.
    pub print_every: usize,
}

impl Default for LoggingOptions {
    fn default() -> Self {
        Self {
            verbose: false,
            print_every: 500,
        }
    }
}

/// Summary of the solver run and final diagnostics.
#[derive(Clone, Debug)]
pub struct SolverResult {
    /// Final primal solution in the original variable coordinates.
    pub x: Vec<f64>,
    /// Final objective value `0.5 x^T Q x + c^T x`.
    pub objective: f64,
    /// Number of first-order iterations executed.
    pub iterations: usize,
    /// Number of accelerated restarts triggered.
    pub num_restarts: usize,
    /// Compact convergence summary for the returned point.
    pub quality: SolverQuality,
    /// Time spent in the first-order and polish phases.
    pub timing: SolverTiming,
    /// Lipschitz constant used internally for the working problem.
    pub lipschitz: f64,
    /// Step size used internally for the working problem.
    pub step_size: f64,
    /// Scaling information for the solve.
    pub scaling: ScalingSummary,
}

/// Core quality diagnostics for a solver run.
#[derive(Clone, Debug)]
pub struct SolverQuality {
    /// Absolute primal-minus-certificate gap based on the affine minorant box bound.
    ///
    /// This is `NaN` when dual certification is disabled.
    pub gap: f64,
    /// Relative primal-dual gap based on the affine minorant box bound.
    ///
    /// This is `NaN` when dual certification is disabled.
    pub rel_gap: f64,
    /// Certified lower bound obtained by minimizing the affine minorant over the box.
    ///
    /// This is `NaN` when dual certification is disabled.
    pub certified_lower_bound: f64,
    /// Infinity-norm box KKT residual.
    pub kkt_inf: f64,
}

/// Timing measurements for the first-order loop and optional polishing.
#[derive(Clone, Debug)]
pub struct SolverTiming {
    /// Time spent in the accelerated projected-gradient loop.
    pub apgd_time_sec: f64,
    /// Time spent in the optional polish phase.
    pub polish_time_sec: f64,
    /// Total solve time across all phases.
    pub total_time_sec: f64,
}

/// Scaling and step-size information for the solve.
#[derive(Clone, Debug)]
pub struct ScalingSummary {
    /// Whether variable scaling was applied.
    pub applied: bool,
    /// Human-readable scaling mode name.
    pub name: &'static str,
    /// Smallest scaling factor used.
    pub scale_min: f64,
    /// Largest scaling factor used.
    pub scale_max: f64,
}

/// Internal diagnostics kept richer than the public result so the solver can
/// make stopping and polish decisions without exposing every research metric.
#[derive(Clone, Debug)]
pub(crate) struct Diagnostics {
    pub kkt_inf: f64,
    pub rel_kkt: f64,
    pub gap: f64,
    pub rel_gap: f64,
    pub certified_lower_bound: f64,
}