pounce_qp/problem.rs
1//! QP problem definition and the solution / warm-start types it
2//! pairs with. Storage borrows from `pounce-linalg` types directly;
3//! the QP solver never copies the Hessian or Jacobian.
4
5use crate::error::{QpError, QpStatus};
6use crate::working_set::WorkingSet;
7use pounce_common::Number;
8use pounce_linalg::triplet::{GenTMatrix, SymTMatrix};
9use std::time::Duration;
10
11/// Caller-supplied hint about the inertia of `H`. Lets a strictly-
12/// convex problem skip the inertia-correction probe of §4.5.
13/// `Unknown` is always safe (the solver detects indefiniteness from
14/// the LDLᵀ factor of the KKT block) and is the default.
15#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
16pub enum HessianInertia {
17 /// `H` is symmetric positive semi-definite.
18 Psd,
19 /// `H` has (potentially) negative eigenvalues; the inertia-
20 /// control path is required.
21 Indefinite,
22 /// Caller offers no claim; solver probes via factor inertia.
23 #[default]
24 Unknown,
25}
26
27/// A convex-or-nonconvex sparse QP:
28/// ```text
29/// min ½ xᵀ H x + gᵀ x
30/// s.t. bl ≤ A x ≤ bu
31/// xl ≤ x ≤ xu
32/// ```
33/// Two-sided general bounds (equality is `bl = bu`); two-sided
34/// variable bounds (fixed is `xl = xu`; free is `±NLP_*_BOUND_INF`).
35/// `H` is symmetric — only the upper triangle is stored.
36///
37/// The lifetime parameter ties the borrowed problem data to the
38/// caller's storage; the solver never copies any of these fields.
39pub struct QpProblem<'a> {
40 pub n: usize,
41 pub m: usize,
42 pub h: &'a SymTMatrix,
43 pub g: &'a [Number],
44 pub a: &'a GenTMatrix,
45 pub bl: &'a [Number],
46 pub bu: &'a [Number],
47 pub xl: &'a [Number],
48 pub xu: &'a [Number],
49 pub hessian_inertia: HessianInertia,
50}
51
52impl<'a> QpProblem<'a> {
53 /// Validate every dimension and bound-ordering invariant the
54 /// solver relies on. Called once at the top of `solve` before
55 /// any work happens.
56 pub fn validate(&self) -> Result<(), QpError> {
57 if self.h.space().dim() as usize != self.n {
58 return Err(QpError::DimensionMismatch(format!(
59 "H is {}×{} but n = {}",
60 self.h.space().dim(),
61 self.h.space().dim(),
62 self.n
63 )));
64 }
65 if self.g.len() != self.n {
66 return Err(QpError::DimensionMismatch(format!(
67 "g.len() = {} but n = {}",
68 self.g.len(),
69 self.n
70 )));
71 }
72 if self.a.space().n_rows() as usize != self.m || self.a.space().n_cols() as usize != self.n
73 {
74 return Err(QpError::DimensionMismatch(format!(
75 "A is {}×{} but expected {}×{}",
76 self.a.space().n_rows(),
77 self.a.space().n_cols(),
78 self.m,
79 self.n
80 )));
81 }
82 if self.bl.len() != self.m || self.bu.len() != self.m {
83 return Err(QpError::DimensionMismatch(format!(
84 "bl.len() = {}, bu.len() = {}, but m = {}",
85 self.bl.len(),
86 self.bu.len(),
87 self.m
88 )));
89 }
90 if self.xl.len() != self.n || self.xu.len() != self.n {
91 return Err(QpError::DimensionMismatch(format!(
92 "xl.len() = {}, xu.len() = {}, but n = {}",
93 self.xl.len(),
94 self.xu.len(),
95 self.n
96 )));
97 }
98 for (i, (&l, &u)) in self.bl.iter().zip(self.bu.iter()).enumerate() {
99 if l > u {
100 return Err(QpError::InvertedBounds(format!(
101 "constraint row {i}: bl = {l} > bu = {u}"
102 )));
103 }
104 }
105 for (i, (&l, &u)) in self.xl.iter().zip(self.xu.iter()).enumerate() {
106 if l > u {
107 return Err(QpError::InvertedBounds(format!(
108 "variable {i}: xl = {l} > xu = {u}"
109 )));
110 }
111 }
112 Ok(())
113 }
114}
115
116/// Warm-start seed: previous primal-dual iterate plus working set.
117/// Passed to [`crate::QpSolver::solve`] as `Some(ws)`; `None` is the
118/// cold-start path through phase-1 elastic mode (§4.3).
119#[derive(Debug, Clone)]
120pub struct QpWarmStart {
121 pub x: Vec<Number>,
122 /// Lagrange multipliers for the general constraints, length `m`.
123 pub lambda_g: Vec<Number>,
124 /// Bound multipliers, length `n`, packed signed
125 /// (`z_l − z_u`). Positive ⇒ lower-bound active, negative ⇒
126 /// upper-bound active.
127 pub lambda_x: Vec<Number>,
128 pub working: WorkingSet,
129}
130
131/// Solver output. `working` is the new working set, suitable for
132/// passing as the next solve's warm start.
133#[derive(Debug, Clone)]
134pub struct QpSolution {
135 pub x: Vec<Number>,
136 pub lambda_g: Vec<Number>,
137 pub lambda_x: Vec<Number>,
138 pub working: WorkingSet,
139 pub obj: Number,
140 pub status: QpStatus,
141 pub stats: QpStats,
142 /// The certified recession direction `d` behind a
143 /// [`QpStatus::Unbounded`] verdict: `Hd ≈ 0`, `d` feasible for every
144 /// step length, `∇q(x)ᵀd < 0`. `None` for every other status.
145 ///
146 /// Callers that embed this QP as a *subproblem* need the witness, not
147 /// just the verdict: an unbounded step QP is only evidence about the
148 /// local model, so an SQP driver has to re-test the ray against the
149 /// true NLP before it can report the NLP unbounded (gh #388). The
150 /// direction is not normalized.
151 pub unbounded_ray: Option<Vec<Number>>,
152}
153
154/// Per-solve counters and timers reported alongside the solution.
155/// Phase 5a uses these for the §8.2 scaling-sweep plots and the
156/// §8.5 warm-start sweep.
157#[derive(Debug, Clone, Default)]
158pub struct QpStats {
159 /// Total active-set changes (adds + drops) across the solve.
160 pub n_working_set_changes: u32,
161 /// Number of times the cached factorization was discarded and
162 /// the base KKT was refactored (the §4.2 reset cycle).
163 pub n_refactor: u32,
164 /// Number of Schur-complement rank-1 updates applied (the
165 /// bounded-cost path between refactors).
166 pub n_schur_updates: u32,
167 /// Whether the solve passed through phase-1 elastic mode.
168 pub used_phase1: bool,
169 /// Wall-clock time spent inside `solve`.
170 pub time: Duration,
171 /// What a [`QpSolver::solve_parametric`] call actually reused from the
172 /// pair it was handed. `None` on every other entry point.
173 ///
174 /// `solve_parametric` has three internal outcomes and only one of them is
175 /// the homotopy: the eligibility guards decline a changed `H` or a changed
176 /// equality/fixed topology, and the tracer itself can return "path could
177 /// not be started". Both fall through to the working-set hint, and a
178 /// previous solve that did not reach `Optimal` falls through again to a
179 /// cold solve. All three return `Ok`, all three are conclusive, and
180 /// nothing in the returned solution distinguished them — so a caller
181 /// counting parametric reuse counted declines as successes, and a
182 /// benchmark asking "is the warm path engaging?" could not be answered
183 /// (gh #769, found in review by @GermanHeim).
184 ///
185 /// [`QpSolver::solve_parametric`]: crate::QpSolver::solve_parametric
186 pub parametric_source: Option<ParametricSource>,
187 /// What the **second-order** test found at the returned point (gh #848).
188 ///
189 /// A summary of [`crate::negcurv::SecondOrder`] without its witness — the
190 /// direction belongs to the escape that consumed it, and when it survives
191 /// as a certificate it survives in [`QpSolution::unbounded_ray`].
192 ///
193 /// Callers that re-derive a status from the returned point need this: a
194 /// saddle point of an indefinite QP satisfies the first-order KKT
195 /// conditions exactly, so *any* verifier built on the KKT residual will
196 /// promote it to `Optimal` no matter what status the engine assigned.
197 /// [`pounce-convex`'s `verify_status`] did precisely that. The field is
198 /// the only channel through which a second-order finding reaches such a
199 /// verifier.
200 ///
201 /// [`pounce-convex`'s `verify_status`]: https://github.com/jkitchin/pounce
202 pub second_order: SecondOrderVerdict,
203}
204
205/// [`crate::negcurv::SecondOrder`] reduced to a verdict, for
206/// [`QpStats::second_order`].
207///
208/// `Default` is [`NotChecked`](SecondOrderVerdict::NotChecked) — the honest
209/// value for every solve that did not run the test, and the reason adding this
210/// field did not have to touch the ~135 `QpSolution` literals in the
211/// workspace: they all build their stats with `..Default::default()`.
212#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
213pub enum SecondOrderVerdict {
214 /// The test did not run, or ran without reaching a conclusion. Never a
215 /// claim that the point is a minimum.
216 #[default]
217 NotChecked,
218 /// The reduced Hessian on the working set's null space is positive
219 /// definite: the point is a local minimum.
220 Certified,
221 /// A feasible direction of strictly negative curvature exists at the
222 /// returned point. It is not a local minimum, at any tolerance.
223 NegativeCurvature,
224}
225
226/// What a [`QpSolver::solve_parametric`](crate::QpSolver::solve_parametric)
227/// call reused from the `(qp_prev, sol_prev)` pair it was handed.
228///
229/// Ordered by how much of the previous solve survived: [`Homotopy`] traces the
230/// path from the previous *solution*, [`WorkingSet`] keeps only the discrete
231/// state, [`Cold`] keeps nothing.
232///
233/// [`Homotopy`]: ParametricSource::Homotopy
234/// [`WorkingSet`]: ParametricSource::WorkingSet
235/// [`Cold`]: ParametricSource::Cold
236#[derive(Debug, Clone, Copy, PartialEq, Eq)]
237pub enum ParametricSource {
238 /// The homotopy path was traced from the previous solution — the
239 /// parametric route proper.
240 Homotopy,
241 /// The homotopy was declined by a guard, or the tracer could not start
242 /// the path. The previous working set was reconciled onto the new problem
243 /// and used as a hint instead, which is worth having even when it is
244 /// stale: what it encodes is *which constraints bind*, and that survives a
245 /// perturbation of `H` far better than the iterate does (gh #602).
246 WorkingSet,
247 /// Nothing from the previous solve was usable — the previous status was
248 /// not `Optimal`, or its working set was dimensionally invalid here — so
249 /// the call solved cold. Also reported when the call was cancelled by the
250 /// deadline before it could choose a route: in both, nothing was reused.
251 Cold,
252}