pounce-algorithm 0.15.0

//! SQP outer-loop options. Mirrors the spirit of
//! `crate::alg_builder::ConvCheckOptions` but only what the SQP
//! driver itself needs; QP-subproblem-specific options pass
//! through `pounce_qp::QpOptions` (which `SqpAlgorithm` owns).

use pounce_common::Number;
use pounce_common::timing::Deadline;
use pounce_qp::HessianInertia;

/// Choice of SQP globalization strategy. Defaults to filter
/// (Fletcher-Leyffer 2002) per the design note §4.1; l1-elastic
/// merit (Han-Powell with adaptive penalty) is the alternative.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum SqpGlobalization {
    #[default]
    Filter,
    L1Elastic,
}

/// Hessian source for the SQP subproblem.
///
/// - `Exact`: use `nlp.eval_h(x, 1.0, λ_g, 0)` directly. Indefinite
///   on nonconvex problems; the QP subproblem solver handles
///   indefinite reduced Hessian via inertia control (§4.5).
/// - `DampedBfgs`: dense `n×n` Powell-damped BFGS, guaranteed PSD
///   (Powell 1978). See [`crate::sqp::bfgs::DampedBfgs`].
/// - `Lbfgs`: limited-memory BFGS with a circular `(s, y)` buffer
///   of `lbfgs_max_history` pairs (Nocedal-Wright §7.2; Byrd-
///   Nocedal-Schnabel 1994). See [`crate::sqp::lbfgs::LBfgs`].
/// - `QuasiNewtonAuto`: what `hessian_approximation=limited-memory`
///   selects — `DampedBfgs` up to [`SqpOptions::dense_hessian_max_n`]
///   variables, `Lbfgs` above it, decided per solve from the problem's
///   `n` (gh#1054). The dense update is the more robust of the two on
///   small models (issue #358) but needs `n²` storage and makes every QP
///   factorization dense, which a 15 000-variable model cannot afford.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum SqpHessianSource {
    #[default]
    Exact,
    DampedBfgs,
    Lbfgs,
    QuasiNewtonAuto,
}

impl SqpHessianSource {
    /// The concrete source a solve of `n` variables uses: resolves
    /// `QuasiNewtonAuto`, passes every other source through.
    pub fn resolve(self, n: usize, dense_hessian_max_n: usize) -> Self {
        match self {
            Self::QuasiNewtonAuto if n <= dense_hessian_max_n => Self::DampedBfgs,
            Self::QuasiNewtonAuto => Self::Lbfgs,
            other => other,
        }
    }

    /// Whether this source approximates the Hessian (and is therefore
    /// positive definite by construction) rather than evaluating it.
    pub fn is_quasi_newton(self) -> bool {
        !matches!(self, Self::Exact)
    }
}

/// Default for [`SqpOptions::qp_schur_min_dim`].
pub const DEFAULT_QP_SCHUR_MIN_DIM: usize = 500;

/// Default for [`SqpOptions::qp_schur_indefinite_max_dim`]. The rule applies
/// to warm step QPs only, so it does not depend on whether cold ones are
/// seeded; the value is the size `algorithm=active-set-sqp` already treats as
/// large (the interior-point handover and cold-start seed start here too),
/// and it matches `pounce_qp::SOFT_DROP_MIN_DIM`, so a warm indefinite QP
/// moved to the refactor path gets the soft drops as well. The reporter's QP is
/// `n + m = 28 950`; nothing measured places the crossover more precisely.
pub const DEFAULT_QP_SCHUR_INDEFINITE_MAX_DIM: usize = 5000;

/// Default for [`SqpOptions::qp_ipm_cold_start_min_dim`]: the size at which
/// `algorithm=active-set-sqp` already starts from the interior point and
/// hands over (`IpoptApplication::sqp_auto_handover`), for the same reason
/// one level down.
pub const DEFAULT_QP_IPM_COLD_START_MIN_DIM: usize = 5000;

#[derive(Debug, Clone)]
pub struct SqpOptions {
    pub globalization: SqpGlobalization,
    pub hessian: SqpHessianSource,

    /// KKT tolerance (max-norm) on the stationarity residual.
    pub tol: Number,
    /// Tolerance on constraint violation (max-norm).
    pub constr_viol_tol: Number,
    /// Tolerance on stationarity residual (max-norm).
    pub dual_inf_tol: Number,
    /// Outer-iteration cap.
    pub max_iter: u32,

    /// l1-merit penalty parameter ν. Used only when
    /// `globalization = L1Elastic`. Filter globalization ignores
    /// this. Default is a moderate starting value; the Han-Powell
    /// ν-adaptation in `l1_merit_line_search` grows ν to dominate
    /// `‖λ_qp‖_∞ + l1_penalty_safety` on every iteration.
    pub l1_penalty: Number,

    /// Additive safety margin in the Han-Powell ν update:
    /// `ν_new = max(ν, ‖λ_qp‖_∞ + l1_penalty_safety)`. Default 0.1
    /// per Nocedal-Wright §18.4. Only consulted when
    /// `globalization = L1Elastic`.
    pub l1_penalty_safety: Number,

    /// Upper clamp on ν. Prevents catastrophic Armijo failure on
    /// pathological problems where `‖λ_qp‖` momentarily spikes.
    /// Default 1e10 — large enough not to interfere with normal
    /// runs, small enough to keep the merit numerically stable.
    /// Only consulted when `globalization = L1Elastic`.
    pub l1_penalty_max: Number,

    /// Backtracking line-search reduction factor.
    pub bt_reduction: Number,
    /// Minimum step before declaring line-search failure.
    pub bt_min_alpha: Number,

    /// `print_level`: 0 = silent, 1 = per-iteration summary,
    /// 2+ = trace (planned).
    pub print_level: u8,

    /// Number of `(s, y)` history pairs retained when
    /// `hessian = Lbfgs`. Mirrors the upstream
    /// `limited_memory_max_history` default of 6 (Nocedal-Wright
    /// recommends 3–20). Ignored for `Exact` and `DampedBfgs`.
    pub lbfgs_max_history: u32,

    /// Largest `n` at which a quasi-Newton Hessian is handed to the QP as a
    /// dense matrix (`sqp_dense_hessian_max_n`). Above it, `Lbfgs` poses the
    /// QP in its lifted low-rank form
    /// ([`crate::sqp::SqpQpData::build_low_rank`]) and `QuasiNewtonAuto`
    /// resolves to `Lbfgs` rather than `DampedBfgs`. At or below it nothing
    /// changes from before gh#1054. Default 1000.
    pub dense_hessian_max_n: usize,

    /// Size from which the step QP uses Schur-complement updates
    /// (gh#1054): when `Some(d)`, each solve sets
    /// `QpOptions::use_schur_updates` to `n + m ≥ d`; `None` leaves the QP
    /// options as given. The driver reads `sqp_qp_use_schur_updates` only
    /// when it is set explicitly, and then passes `None` here.
    ///
    /// Why size: the refactor path pays a full sparse factorization per
    /// working-set change, which is what made the cold first QP of a
    /// 15 000-variable model run for 30 minutes, while on small models the
    /// two paths cost the same and the Schur path is the less robust one
    /// (measured on the CLI fixture corpus under `algorithm=active-set-sqp`).
    /// See `docs/src/active-set-sqp.md` for the Maros-Mészáros comparison
    /// behind the default.
    pub qp_schur_min_dim: Option<usize>,

    /// Size from which a *warm-started* step QP whose Hessian is not known
    /// PSD goes back to the refactor path, overriding
    /// [`Self::qp_schur_min_dim`] (gh#1083, gh#1091). `None` never does; cold
    /// solves and PSD step QPs are not affected.
    ///
    /// Such a QP carries the inertia hint `Indefinite`, and on it every
    /// *drop* refactors anyway (a rank-2 update cannot re-run inertia control
    /// on the enlarged null space), so only adds take the update. A warm step
    /// QP is drop-dominated: the reporter's took 187 Schur updates against
    /// 4181 working-set changes over a 3600 s run. It then pays for the larger
    /// `K_max` — dearer factorizations, back-solves about 5× the small KKT's —
    /// and gets almost nothing back; on its handover QP the Schur path took
    /// 126 s, the refactor path 72 s.
    ///
    /// Two step QPs carry the hint: an exact Hessian's, and the lifted
    /// L-BFGS one above `sqp_dense_hessian_max_n`
    /// ([`crate::sqp::SqpQpData::build_low_rank`]), whose Hessian
    /// `diag(σI, ±1)` is indefinite although the `B` it lifts is not. The rule
    /// reads the hint, which is what the QP solver branches on, not the
    /// Hessian source; until gh#1091 it read the source and missed the lift.
    ///
    /// Cold solves keep the size rule because they are the add-dominated
    /// ones the Schur path wins on, and at this size they are not all seeded
    /// from the interior point: when the hint is not `Psd` the seed certifies
    /// the Hessian PSD numerically (`ipm_seeded_step_qp`) and declines if it
    /// is not. The lifted QP's `-1` entries never pass, and a nonconvex
    /// model's exact Hessian need not, so such a cold QP still builds its
    /// working set one add at a time. (A convex model's exact Hessian does
    /// pass, and its cold step QPs are seeded.)
    pub qp_schur_indefinite_max_dim: Option<usize>,

    /// NLP size `n + m` from which a step QP solved without a working set,
    /// or with one that pins to an infeasible point, is first solved by the
    /// convex interior point, when its Hessian is PSD, and the active-set
    /// solver warm-started from the working set that solution identifies
    /// (gh#1061). `None` never does. The warm case: the reporter's second
    /// step QP went to the ℓ1-elastic phase 1 from its infeasible warm start
    /// and reached a 300 s cap after 564 working-set changes; seeded, 12.5 s.
    ///
    /// A cold active-set start builds its working set one pivot at a time,
    /// and on a large model that is the whole cost: gh#1061's first step QP
    /// (`n = 14 952`, 10 818 equalities, `H = I`) went through the elastic
    /// phase 1 for more than 30 minutes on either factorization path without
    /// finishing. The interior point solved it in 16 iterations, 4.3 s, and
    /// the active-set solve from its working set took zero working-set
    /// changes and one factorization, 0.4 s, to the same objective. The
    /// active-set solve still has the last word, so the step and its
    /// multipliers are what the active-set method certifies, and anything
    /// short of an interior-point `Optimal` falls back to the cold start.
    ///
    /// At the same size a warm start that pins feasibly is given a budget of
    /// four active-set iterations per iteration the last seed took, and is
    /// seeded once it runs over. gh#1061's warm step QPs took 68–350
    /// working-set changes and 51–263 s each at about one factorization per
    /// change; on the handover QP the budget and the seed together take 34 s
    /// against 72 s warm (126 s with Schur updates).
    pub qp_ipm_cold_start_min_dim: Option<usize>,

    /// Solve-wide time budget (`max_wall_time` / `max_cpu_time`), the
    /// same [`Deadline`] the interior-point arm checks. `None` means
    /// unlimited and is the default for a directly constructed driver.
    ///
    /// A shared handle rather than two durations, because the SQP is not
    /// always the first thing a solve runs: crossover hands it the
    /// interior-point phase's own deadline, and a duration measured from
    /// `optimize` would grant it a fresh budget (gh#1054).
    pub deadline: Option<Deadline>,
}

impl Default for SqpOptions {
    fn default() -> Self {
        Self {
            globalization: SqpGlobalization::default(),
            hessian: SqpHessianSource::default(),
            tol: 1e-8,
            constr_viol_tol: 1e-6,
            dual_inf_tol: 1e-4,
            max_iter: 200,
            l1_penalty: 1.0,
            l1_penalty_safety: 0.1,
            l1_penalty_max: 1e10,
            bt_reduction: 0.5,
            bt_min_alpha: 1e-12,
            print_level: 0,
            lbfgs_max_history: 6,
            dense_hessian_max_n: 1000,
            qp_schur_min_dim: Some(DEFAULT_QP_SCHUR_MIN_DIM),
            qp_schur_indefinite_max_dim: Some(DEFAULT_QP_SCHUR_INDEFINITE_MAX_DIM),
            qp_ipm_cold_start_min_dim: Some(DEFAULT_QP_IPM_COLD_START_MIN_DIM),
            deadline: None,
        }
    }
}

impl SqpOptions {
    /// Whether a step QP of size `n + m = dim` uses Schur updates, by
    /// [`Self::qp_schur_min_dim`] and [`Self::qp_schur_indefinite_max_dim`].
    /// `warm` is the inertia hint of a QP warm-started from a working set,
    /// `None` for a cold one, whose hint the rule does not read. `None` back
    /// when the caller chose (`qp_schur_min_dim` unset), which leaves the QP
    /// options as given.
    pub fn step_qp_schur_updates(&self, dim: usize, warm: Option<HessianInertia>) -> Option<bool> {
        let min = self.qp_schur_min_dim?;
        // `!= Psd`, as in the QP solver's own drop branch: `Unknown` resets on
        // a drop just as `Indefinite` does.
        let refactor = warm.is_some_and(|h| h != HessianInertia::Psd)
            && self.qp_schur_indefinite_max_dim.is_some_and(|d| dim >= d);
        Some(dim >= min && !refactor)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    /// gh#1083 / gh#1091: from the ceiling on, a warm-started step QP not
    /// known PSD — an exact Hessian's, or the lifted L-BFGS one — is back on
    /// the refactor path. A cold one keeps its updates (it is the
    /// add-dominated case, and not seeded when not PSD), as does a PSD one,
    /// and an explicit choice (`qp_schur_min_dim = None`) is left alone.
    #[test]
    fn warm_indefinite_step_qps_leave_schur_updates_off_at_size() {
        use HessianInertia::{Indefinite, Psd, Unknown};
        let o = SqpOptions::default();
        let (lo, hi) = (
            DEFAULT_QP_SCHUR_MIN_DIM,
            DEFAULT_QP_SCHUR_INDEFINITE_MAX_DIM,
        );
        assert!(lo < hi);
        for warm in [None, Some(Indefinite), Some(Psd), Some(Unknown)] {
            let at = |d| o.step_qp_schur_updates(d, warm);
            assert_eq!(at(lo - 1), Some(false), "{warm:?}");
            assert_eq!(at(lo), Some(true), "{warm:?}");
            assert_eq!(at(hi - 1), Some(true), "{warm:?}");
        }
        for h in [Indefinite, Unknown] {
            assert_eq!(o.step_qp_schur_updates(hi, Some(h)), Some(false), "{h:?}");
        }
        assert_eq!(o.step_qp_schur_updates(hi, None), Some(true), "cold");
        assert_eq!(o.step_qp_schur_updates(hi, Some(Psd)), Some(true));

        let unset = SqpOptions {
            qp_schur_indefinite_max_dim: None,
            ..SqpOptions::default()
        };
        assert_eq!(
            unset.step_qp_schur_updates(hi, Some(Indefinite)),
            Some(true)
        );
        let chosen = SqpOptions {
            qp_schur_min_dim: None,
            ..SqpOptions::default()
        };
        assert_eq!(chosen.step_qp_schur_updates(hi, Some(Indefinite)), None);
    }
}