diff --git a/CHANGELOG.md b/CHANGELOG.md index f65d99cf..7e9b3975 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1372,38 +1372,48 @@ Both are detailed under *Behaviour changes to plan for*. than its interior, which a plain fixed-floor search reliably stalls short of. - **Honest gap, not a silent one:** the hardest classical textbook examples — - Motzkin's polynomial (`x⁴y²+x²y⁴−3x²y²+1`, PSD but not SOS; classically SOS - after multiplying by `x²+y²`) and Robinson's form — are *not* yet - certified by this search, and the tests say so directly - (`motzkin_reports_undecided_rather_than_a_false_certificate`, - `psd_search_does_not_yet_reach_homogeneous_motzkin_times_sum_of_squares`) - rather than asserting a false positive. This was diagnosed, not merely - observed: a diagnostic trajectory - (`psd::diag::diag_step1_step2_trajectory_and_family_sanity`) shows the - annealed search converging *monotonically* toward Motzkin's boundary - certificate (minimum eigenvalue running from roughly `−1.6` to roughly - `−0.0018` as the floor anneals to `0`) without fully closing the last, - asymptotically slow stretch to exactly `0` — the textbook signature of - alternating projection at a tangential (non-transversal) set intersection, - a known hard case for this class of method, not a bug. That the mechanism - itself is sound was checked independently two ways: an exact sanity check - that the affine Gram-matrix family constructed for Motzkin really does - reproduce the target polynomial at an arbitrary rational point in the - family, and a synthetic planted rank-deficient PSD example of the same - nullspace dimension (`psd::diag::diag_step3_planted_singular_example`), - which *is* found and exactly re-verified. Recording `undecided` on Motzkin - and Robinson is the correct behaviour for now, not a workaround; escaping - a tangential intersection reliably (e.g. Douglas–Rachford with - over-relaxation, or a facial-reduction preprocessing step) is future work. + **Motzkin and Robinson now both certify — the gap that motivated this + feature is closed, though not unconditionally.** The first version of this + search (annealed alternating projection with several random restarts) + reliably fell short of both classical textbook examples: a diagnostic + trajectory (`psd::diag::diag_step1_step2_trajectory_and_family_sanity`) + showed it converging *monotonically* toward Motzkin's boundary certificate + (minimum eigenvalue running from roughly `−1.6` to roughly `−0.0018` as the + floor annealed to `0`) without ever closing the last, asymptotically slow + stretch to exactly `0` — the textbook signature of alternating projection + stalling at a *tangential* (non-transversal) set intersection, which is + exactly what a *singular* witnessing Gram matrix (sitting on the PSD cone's + boundary rather than its interior) produces. The search now also tries + Douglas–Rachford splitting with over-relaxation + (`sdp::Family::douglas_rachford_from`) and a facial-reduction step + (`psd::facial_reduction_search`) — both are standard escapes for exactly + this stall — and with them, both `(x²+y²)·Motzkin(x,y)` (the affine, + 2-variable case, found via the full `sos_decompose` multiplier search, not + just a hand-fed pre-multiplied target) and `(x²+y²+z²)·Robinson(x,y,z)` + (via `psd_search` directly) are found and exactly re-verified — + `real::sos::tests::motzkin_certifies_via_a_reznick_multiplier` and + `psd::tests::psd_search_certifies_robinsons_form_with_a_reznick_multiplier` + check the identities by hand, independent of the search that proposed them. + + **What's still open:** the homogeneous 3-variable form of Motzkin, + `(x²+y²+z²)·(x⁴y²+x²y⁴−3x²y²z²+z⁶)`, still is not found + (`psd::tests::psd_search_does_not_yet_reach_homogeneous_motzkin_times_sum_of_squares`) + — a larger nullspace than the affine 2-variable case, and evidently still + hard enough for even Douglas–Rachford and facial reduction as currently + tuned. So a boundary-only certificate is not guaranteed to be found in + general; `E-SOS-002` still means "not found within this search", never "not + SOS" or "not non-negative". Everything reachable today is exact end to end: `verify()` re-expands every returned certificate with exact rational arithmetic, `to_lean()` - emits a sorry-free Lean sketch, and `PositivityCertificate.multiplier` is + emits a sorry-free Lean sketch, and `PositivityCertificate.multiplier()` is populated exactly when the certificate needed one (`None` for a direct SOS - decomposition). No new public API surface — `sos_decompose` and - `PositivityCertificate` are unchanged in shape; this is entirely a - strengthening of what the existing search covers before it refuses. + decomposition — a method rather than a field, since adding a field to this + already fully-public, exhaustively-constructible struct is a semver break + regardless of the field's own visibility; see the method's own doc comment). + No new public API surface — `sos_decompose` and `PositivityCertificate` are + unchanged in shape; this is entirely a strengthening of what the existing + search covers before it refuses. ### Performance diff --git a/alkahest-core/src/real/sos/mod.rs b/alkahest-core/src/real/sos/mod.rs index 37acec87..6aec0070 100644 --- a/alkahest-core/src/real/sos/mod.rs +++ b/alkahest-core/src/real/sos/mod.rs @@ -36,12 +36,20 @@ //! subcone (solvable exactly); the full PSD Gram cone, when DSOS fails (a //! strict superset, but only reachable via the sound-but-incomplete numeric //! search above); and a Reznick multiplier search `(Σxᵢ²)^N·p`, when even -//! that fails on `p` itself. None of these three is complete — the multiplier -//! search in particular does not yet reliably find certificates whose -//! witnessing Gram matrix is singular (sits exactly on the PSD cone's -//! boundary), which is the case for the textbook examples Motzkin and -//! Robinson (see `real::sos::tests::motzkin_reports_undecided_rather_than_a_false_certificate` -//! for the diagnosis). So [`SosError::NoCertificate`] means precisely *"no +//! that fails on `p` itself. None of these three is complete — a certificate +//! of a given shape may exist at a higher degree/budget than was searched, or +//! not exist in that shape at all, and `sos_decompose` cannot tell those +//! apart. The multiplier search specifically had to add Douglas–Rachford +//! splitting and facial reduction alongside its original annealed +//! alternating-projection search (see `real::sos::sdp::Family::douglas_rachford_from` +//! and `real::sos::psd`'s `facial_reduction_search`) because certificates +//! whose witnessing Gram matrix is *singular* — sitting exactly on the PSD +//! cone's boundary, as for the textbook examples Motzkin and Robinson — are a +//! well-known hard case for plain alternating projection (see +//! `real::sos::tests::motzkin_certifies_via_a_reznick_multiplier` for the +//! worked diagnosis); both of those examples are covered by the current +//! search, but a boundary-only certificate at some other degree is not +//! guaranteed to be. So [`SosError::NoCertificate`] means precisely *"no //! certificate of this shape was found at this degree/budget"*. It does //! **not** mean "not a sum of squares", and it does **not** mean "not //! non-negative". The three answers are kept distinct in the API on purpose — @@ -621,55 +629,63 @@ mod tests { } #[test] - fn motzkin_reports_undecided_rather_than_a_false_certificate() { + fn motzkin_certifies_via_a_reznick_multiplier() { let (pool, x, y) = setup(); // Motzkin: x^4·y^2 + x^2·y^4 − 3·x^2·y^2 + 1 is non-negative but is // the textbook example of a polynomial that is *not* itself a sum of // squares — Hilbert's 1888 theorem allows non-SOS PSD forms outside // ternary quartics, and Motzkin (1967) is the standard witness. // Multiplying by (x²+y²) is classically known to fix this (it is - // exactly the kind of case Reznick's theorem covers), but the - // witnessing Gram matrix for that fact is *singular* — it sits - // exactly on the boundary of the PSD cone, not in its interior — and - // [`crate::real::sos::psd`]'s numeric search (alternating projection - // with an annealed floor schedule and multiple random restarts) is - // demonstrably not a bug: a `psd::diag::diag_step3_planted_singular_example` - // planted boundary case with the same nullspace dimension *is* found - // and exactly re-verified, and the affine family constructed for - // Motzkin itself passes an independent sanity check - // (`psd::diag::diag_step1_step2_trajectory_and_family_sanity`). The - // search on Motzkin specifically converges monotonically (min - // eigenvalue runs from roughly −1.6 down to roughly −0.0018 as the - // floor anneals to 0) but does not close the last, asymptotically - // slow stretch to exactly 0 — the classic behaviour of alternating - // projection at a tangential (non-transversal) intersection. This is - // an honest search-budget limitation, not a soundness bug: recording - // `undecided` here, never a fabricated certificate, is the correct - // behaviour and is what this test checks. + // exactly the kind of case Reznick's theorem covers). The witnessing + // Gram matrix for that fact is *singular* — it sits exactly on the + // boundary of the PSD cone, not in its interior — which is why a + // plain annealed alternating-projection search (the first version of + // this feature) could get arbitrarily close (min eigenvalue from + // roughly −1.6 down to roughly −0.0018 as its floor annealed to 0) + // without ever closing the gap: the textbook symptom of a tangential + // (non-transversal) set intersection. `real::sos::psd`'s search now + // also tries Douglas–Rachford splitting and facial reduction, either + // of which is known to escape exactly this kind of stall, and it + // finds Motzkin's certificate here. let p = pool.add(vec![ pool.mul(vec![x, x, x, x, y, y]), pool.mul(vec![x, x, y, y, y, y]), pool.mul(vec![pool.integer(-3_i32), x, x, y, y]), pool.integer(1_i32), ]); - let err = sos_decompose(p, &[x, y], &pool, &SosOpts::default()) - .expect_err("Motzkin's multiplier certificate is not yet reached by this search"); - assert!(matches!(err, SosError::NoCertificate(_))); - assert_eq!(err.code(), "E-SOS-002"); + let cert = + sos_decompose(p, &[x, y], &pool, &SosOpts::default()).expect("Motzkin now certifies"); + assert_eq!(cert.kind, CertificateKind::Sos); + let sigma = cert + .multiplier() + .expect("Motzkin is not itself SOS, so this must be a multiplier certificate"); + // The exact identity actually checked: σ·p = Σ c_i q_i², in ℚ — the + // real soundness argument, independent of whatever the numeric + // search that proposed it actually converged to. + assert_eq!(sigma.mul(&cert.target), cert.expand()); + cert.verify().expect("re-verifies exactly end to end"); + + // Composes with to_lean: a self-contained, sorry-free Lean sketch. + let lean = cert + .to_lean() + .expect("multiplier certificates emit Lean too"); + assert!(!lean.contains("sorry")); + assert!(!lean.contains("admit")); + assert!(lean.contains("alkahest_multiplier_factor")); + assert!(lean.contains("ring")); } #[test] fn multiplier_search_reports_undecided_not_not_sos_when_out_of_budget() { let (pool, x, y) = setup(); - // Same Motzkin target as the previous test. Here the internal search - // is driven with a budget of *zero* multiplier powers directly — i.e. - // exactly the "search legitimately runs out of budget" case — and it - // must come back empty-handed rather than fabricate a certificate. - // (Motzkin also fails to certify at the production budget, per the - // previous test — this test's point is narrower: even independent of - // whether the production budget eventually finds Motzkin's - // certificate, a caller-supplied budget of zero must never - // manufacture one.) + // Same Motzkin target as the previous test, which now certifies at + // the production budget. Here the internal search is instead driven + // with a budget of *zero* multiplier powers directly — i.e. exactly + // the "search legitimately runs out of budget" case — and it must + // come back empty-handed rather than fabricate a certificate. This + // test's point is narrower than the previous one: independent of + // whether the production budget finds a given target's certificate, + // a caller-supplied budget of zero must never manufacture one. let p = pool.add(vec![ pool.mul(vec![x, x, x, x, y, y]), pool.mul(vec![x, x, y, y, y, y]), diff --git a/alkahest-core/src/real/sos/psd.rs b/alkahest-core/src/real/sos/psd.rs index 4fe05b9b..40d27153 100644 --- a/alkahest-core/src/real/sos/psd.rs +++ b/alkahest-core/src/real/sos/psd.rs @@ -33,7 +33,7 @@ use super::cert::SosPoly; use super::gram::monomial_basis; use super::linalg::{psd_decompose, solve_affine}; use super::ratpoly::{Exponents, RatPoly}; -use super::sdp::{min_eigenvalue, Family}; +use super::sdp::{min_eigenvalue, smallest_magnitude_eigenvectors, Family}; use rug::Rational; use std::collections::{BTreeMap, BTreeSet}; @@ -307,9 +307,49 @@ const ROUNDING_CANDIDATES: usize = 6; /// this keeps that bounded regardless of how large a monomial basis the /// caller asks for. A skip here returns `None` — "not found within /// budget", not "not SOS" — exactly like every other budget in this module. -const MAX_FREE_PARAMETERS: usize = 110; +const MAX_FREE_PARAMETERS: usize = 200; -/// Run the annealing schedule from a single starting point. +/// Over-relaxation parameters tried for the Douglas–Rachford polish, in +/// increasing order of overshoot. `1.0` is plain (non-relaxed) +/// Douglas–Rachford; the larger value is the standard mitigation for a +/// stall at a shallow tangential approach to the intersection (see +/// `Family::douglas_rachford_from`'s doc comment). Kept to two values +/// deliberately: [`DR_ITERS`] below is what actually closes the gap on hard +/// (boundary-only) instances — see its doc comment — and that only stays +/// affordable across [`DR_POLISH_CANDIDATES`] starting points *and* several +/// facial-reduction attempts ([`FACIAL_SEARCH_BUDGET`]) if the spread here +/// is kept small. +const DR_LAMBDAS: &[f64] = &[1.0]; + +/// Iterations run per Douglas–Rachford attempt. Reflections make *some* +/// progress every iteration on a boundary-only (tangential) intersection, +/// but the rate is still only sublinear there — a diagnostic run on the +/// homogeneous Motzkin family (165 free parameters) needed on the order of +/// `10^4` iterations to bring the minimum eigenvalue from `~-2·10⁻³` (where +/// the alternating-projection annealing schedule alone stalls) down to +/// `~-5·10⁻⁶`, which is close enough that rational rounding at +/// [`DENOM_CAPS`]'s larger denominators reliably lands exactly on the true +/// (small-denominator) certificate. This is run only on the best few +/// annealed candidates (see [`DR_POLISH_CANDIDATES`]), not every start, to +/// keep that cost affordable. +const DR_ITERS: usize = 15_000; + +/// How many of the (cheap, alternating-projection-only) annealed candidates +/// get the expensive Douglas–Rachford polish. Bounded well below the total +/// number of starts multistart annealing tries — see [`DR_ITERS`] for why +/// the polish itself is not cheap — since a candidate's annealed minimum +/// eigenvalue is already a good proxy for whether it is worth polishing: +/// the polish improves a near-feasible point, it does not rescue a +/// genuinely bad one. +const DR_POLISH_CANDIDATES: usize = 2; + +/// Run the alternating-projection annealing schedule from a single starting +/// point. This alone reaches most cases outright, and gets *close* even on +/// the hard boundary-only ones (on the homogeneous Motzkin case +/// specifically, `diag::diag_step1_step2_trajectory_and_family_sanity` +/// shows it running the minimum eigenvalue from about −1.6 to about +/// −0.0018 and no further) — closing that last "close but stalled" stretch +/// is what the Douglas–Rachford polish in [`multistart_anneal`] is for. fn anneal_from(family: &Family, start: Vec) -> Vec { let mut t = start; for &floor in FLOOR_SCHEDULE { @@ -321,13 +361,27 @@ fn anneal_from(family: &Family, start: Vec) -> Vec { } /// Try the annealing schedule from several starting points — the -/// deterministic `t = 0`, plus a handful of random restarts — and return -/// every result reached, best (highest minimum eigenvalue of `family.at(t)`) -/// first. See [`FLOOR_SCHEDULE`]'s and [`RANDOM_RESTARTS`]'s doc comments -/// for why both are needed: annealing handles boundary-only intersections -/// that a fixed floor stalls on, and multiple starts hedge against any -/// single trajectory converging to a merely-locally-nearest pair when the -/// family and the PSD cone do intersect elsewhere. +/// deterministic `t = 0`, plus a handful of random restarts — then hand the +/// best few results to a deep Douglas–Rachford polish (see +/// `Family::douglas_rachford_from` and [`DR_ITERS`]'s doc comment), and +/// return every result reached, best (highest minimum eigenvalue of +/// `family.at(t)`) first. +/// +/// See [`FLOOR_SCHEDULE`]'s and [`RANDOM_RESTARTS`]'s doc comments for why +/// both annealing and multiple starts are needed: annealing handles +/// boundary-only intersections that a fixed floor stalls on, and multiple +/// starts hedge against any single trajectory converging to a merely +/// locally-nearest pair when the family and the PSD cone do intersect +/// elsewhere. The Douglas–Rachford polish is applied only to the best few +/// annealed candidates, not every start — it is the standard reflection-based +/// upgrade for exactly the "close but stalled" signature annealing alone +/// leaves on a *tangential* (non-transversal) intersection (a singular +/// witnessing Gram matrix — the case for tight certificates like Motzkin's), +/// but running it to the depth that actually closes such a gap ([`DR_ITERS`]) +/// is too expensive to spend on every start indiscriminately; a candidate's +/// annealed eigenvalue is already a good proxy for which ones are worth it. +/// Polishing can never make a candidate worse — its own annealed point is +/// always kept as a floor. fn multistart_anneal(family: &Family, dim: usize) -> Vec> { let mut starts: Vec> = vec![vec![0.0; dim]]; let mut rng = SplitMix64::new(0xC0FFEE_D15EA5E5); @@ -345,6 +399,19 @@ fn multistart_anneal(family: &Family, dim: usize) -> Vec> { }) .collect(); results.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); + + for (eig, t) in results.iter_mut().take(DR_POLISH_CANDIDATES) { + for &lambda in DR_LAMBDAS { + if let Some(cand) = family.douglas_rachford_from(t.clone(), 0.0, lambda, DR_ITERS) { + let cand_eig = min_eigenvalue(&family.at(&cand)); + if cand_eig > *eig { + *eig = cand_eig; + *t = cand; + } + } + } + } + results.sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal)); results.into_iter().map(|(_, t)| t).collect() } @@ -360,40 +427,52 @@ fn multistart_anneal(family: &Family, dim: usize) -> Vec> { /// always sound: the returned [`SosPoly`] is checked to expand back to /// exactly `p` before it is returned, using the same exact rational /// arithmetic as everywhere else in this subsystem. -pub fn psd_search(target: &RatPoly, basis_deg: u32) -> Option { - let nvars = target.nvars(); - // A homogeneous target of degree exactly `2·basis_deg` needs only the - // monomials of degree *exactly* `basis_deg` in its Gram basis — mixing in - // lower-degree monomials can only ever contribute to coefficients the - // target does not have, since every product of two basis monomials of - // unequal degree still sums to `2·basis_deg` only when *both* already - // have degree `basis_deg`. This is standard (Blekherman–Parrilo–Thomas, - // Prop. 3.29): a homogeneous SOS decomposition can always be taken with - // homogeneous summands. Restricting here is not just an optimisation — - // the search is numeric, and a smaller basis is the difference between - // "converges" and "not within budget" on cases like Motzkin. - let basis: Vec = match target.is_homogeneous() { - Some(d) if d == 2 * basis_deg => monomial_basis(nvars, basis_deg) - .into_iter() - .filter(|e| e.iter().sum::() == basis_deg) - .collect(), - _ => monomial_basis(nvars, basis_deg), - }; - let n = basis.len(); - if n == 0 { - return if target.is_zero() { - Some(SosPoly::default()) - } else { - None - }; +/// Build `base + Σ t_k·dirs[k]` exactly, for a rational affine family whose +/// members are already unpacked `n×n` matrices (as opposed to +/// `linalg::AffineSolution::at`, which works in the packed upper-triangle +/// representation `gram_system` uses). +fn rat_family_at( + base: &[Vec], + dirs: &[Vec>], + t: &[Rational], +) -> Vec> { + let n = base.len(); + let mut q = base.to_vec(); + for (tk, dir) in t.iter().zip(dirs.iter()) { + if *tk == 0 { + continue; + } + for i in 0..n { + for j in 0..n { + if dir[i][j] != 0 { + q[i][j] += Rational::from(tk * &dir[i][j]); + } + } + } } + q +} - let (rows, rhs) = gram_system(target, &basis); - let sol = solve_affine(&rows, &rhs)?; - +/// Search the exact rational affine family `base_rat + Σ t_k·dirs_rat[k]` +/// (already unpacked `n×n` matrices) for a point that is both PSD and +/// reproduces `target` over `basis` — the shared core of [`psd_search`], +/// factored out so [`facial_reduction_search`] can run it again on a +/// *smaller* family without duplicating the numeric-search/rounding logic. +/// +/// Returns the certificate if found, together with the best numeric +/// candidate matrices tried (nearest-to-feasible first) regardless of +/// whether a certificate was found — [`facial_reduction_search`] reads +/// near-null eigenvectors off those candidates, so they are worth handing +/// back even from a `None` result. +fn search_rational_family( + nvars: usize, + basis: &[Exponents], + target: &RatPoly, + base_rat: &[Vec], + dirs_rat: &[Vec>], +) -> (Option, Vec>>) { let try_point = |t: &[Rational]| -> Option { - let packed = sol.at(t); - let q = unpack(n, &packed); + let q = rat_family_at(base_rat, dirs_rat, t); let decomp = psd_decompose(&q)?; let mut sos = SosPoly::default(); for (d, v) in decomp { @@ -420,23 +499,21 @@ pub fn psd_search(target: &RatPoly, basis_deg: u32) -> Option { }; // No freedom at all: the unique solution is the only candidate. - if sol.dimension() == 0 { - return try_point(&[]); + if dirs_rat.is_empty() { + return (try_point(&[]), Vec::new()); } - if sol.dimension() > MAX_FREE_PARAMETERS { - return None; + if dirs_rat.len() > MAX_FREE_PARAMETERS { + return (None, Vec::new()); } - let base: Vec> = unpack(n, &sol.particular) + let base: Vec> = base_rat .iter() .map(|row| row.iter().map(rat_to_f64).collect()) .collect(); - let dirs: Vec>> = sol - .nullspace + let dirs: Vec>> = dirs_rat .iter() .map(|dir| { - unpack(n, dir) - .iter() + dir.iter() .map(|row| row.iter().map(rat_to_f64).collect()) .collect() }) @@ -448,11 +525,11 @@ pub fn psd_search(target: &RatPoly, basis_deg: u32) -> Option { let family = Family::new(base, ortho_dirs); let dim = dirs.len(); - for s in multistart_anneal(&family, dim) - .into_iter() - .take(ROUNDING_CANDIDATES) - { - let Some(t) = back_substitute_upper(&r, &s) else { + let candidates = multistart_anneal(&family, dim); + let candidate_matrices: Vec>> = candidates.iter().map(|s| family.at(s)).collect(); + + for s in candidates.iter().take(ROUNDING_CANDIDATES) { + let Some(t) = back_substitute_upper(&r, s) else { continue; }; for &max_den in DENOM_CAPS { @@ -460,13 +537,240 @@ pub fn psd_search(target: &RatPoly, basis_deg: u32) -> Option { t.iter().map(|x| round_to_rational(*x, max_den)).collect(); let Some(t_rat) = t_rat else { continue }; if let Some(sos) = try_point(&t_rat) { - return Some(sos); + return (Some(sos), candidate_matrices); + } + } + } + (None, candidate_matrices) +} + +/// Number of near-null eigenvector directions imposed at once, tried in +/// increasing order — the guessed *corank* of the true certificate's +/// witnessing Gram matrix. +const FACIAL_CORANK_GUESSES: &[usize] = &[1, 2, 3]; + +/// Denominator caps tried when rounding a near-null eigenvector to an exact +/// rational direction. Coarser than [`DENOM_CAPS`]: facial reduction is a +/// *guess* at the true certificate's nullspace, and a wrong guess is cheap +/// to detect (`solve_affine`, or the final exact re-expansion check, simply +/// finds nothing) — so there is no benefit in trying many denominators here, +/// and [`DENOM_CAPS`]'s full spread is tried on the *reduced* family once a +/// guess produces one. +const FACIAL_DENOM_CAPS: &[i64] = &[64, 4096]; + +/// How many of the best numeric candidates (matrices) to read near-null +/// directions off of. +const FACIAL_CANDIDATES: usize = 2; + +/// Total number of reduced-family searches [`facial_reduction_search`] will +/// actually run (across every candidate/corank/denominator combination +/// tried) — each one reruns the full numeric search machinery on a smaller +/// family, so this bounds the added cost to a small multiple of one +/// [`search_rational_family`] call, regardless of how many guesses are +/// considered. +const FACIAL_SEARCH_BUDGET: usize = 3; + +/// Given a rational vector `v` — a *guessed* near-null direction of the +/// affine family's true PSD certificate — impose `Q(t)·v = 0` on the family +/// `base_rat + Σ t_k·dirs_rat[k]`, and return the resulting (generally +/// smaller) affine family in the same unpacked `n×n`-matrix representation. +/// `None` means the guess is inconsistent with the family — a routine "try +/// a different guess" outcome, not a bug: `dirs_rat`'s directions already +/// satisfy `z^T Q z = target` for any `t` by construction, and imposing +/// `Q(t)·v = 0` only ever cuts that same solution set down (to the sub-family +/// where a specific vector actually is a nullspace direction), never +/// enlarges or corrupts it — so any family this returns still reproduces +/// `target` exactly, for the same reason the original one does. +#[allow(clippy::type_complexity)] +fn restrict_by_near_null_vectors( + base_rat: &[Vec], + dirs_rat: &[Vec>], + vs: &[Vec], +) -> Option<(Vec>, Vec>>)> { + let n = base_rat.len(); + let dim = dirs_rat.len(); + let mut rows = Vec::with_capacity(vs.len() * n); + let mut rhs = Vec::with_capacity(vs.len() * n); + for v in vs { + for i in 0..n { + let mut row = vec![Rational::from(0); dim]; + for (k, dir) in dirs_rat.iter().enumerate() { + let mut acc = Rational::from(0); + for j in 0..n { + if dir[i][j] != 0 && v[j] != 0 { + acc += Rational::from(&dir[i][j] * &v[j]); + } + } + row[k] = acc; + } + let mut b = Rational::from(0); + for j in 0..n { + if base_rat[i][j] != 0 && v[j] != 0 { + b += Rational::from(&base_rat[i][j] * &v[j]); + } + } + rows.push(row); + rhs.push(-b); + } + } + let sol = solve_affine(&rows, &rhs)?; + + let mut new_base = base_rat.to_vec(); + for (k, tk) in sol.particular.iter().enumerate() { + if *tk == 0 { + continue; + } + for i in 0..n { + for j in 0..n { + if dirs_rat[k][i][j] != 0 { + new_base[i][j] += Rational::from(tk * &dirs_rat[k][i][j]); + } + } + } + } + let mut new_dirs = Vec::with_capacity(sol.nullspace.len()); + for coeffs in &sol.nullspace { + let mut nd = vec![vec![Rational::from(0); n]; n]; + for (k, ck) in coeffs.iter().enumerate() { + if *ck == 0 { + continue; + } + for i in 0..n { + for j in 0..n { + if dirs_rat[k][i][j] != 0 { + nd[i][j] += Rational::from(ck * &dirs_rat[k][i][j]); + } + } + } + } + new_dirs.push(nd); + } + Some((new_base, new_dirs)) +} + +/// Facial-reduction fallback for when [`search_rational_family`]'s direct +/// numeric search gets *close* to a PSD point but never lands one that +/// survives exact rational rounding — the textbook signature of a +/// boundary-only (singular) witnessing Gram matrix, which is exactly what +/// makes both plain alternating projection and Douglas–Rachford converge +/// only asymptotically rather than in finitely many steps (see the module +/// doc and `Family::douglas_rachford_from`'s doc comment). +/// +/// The fix, rather than searching harder in the same (degenerate) family, is +/// to search a *smaller* one: read the near-null eigenvector directions off +/// the best numeric candidates found so far (`smallest_magnitude_eigenvectors`), +/// round each to an exact rational, and impose it as an exact `Q(t)·v = 0` +/// constraint. If the guess is right, this restricts the family to (an +/// affine reparametrisation of) the face of the PSD cone the true +/// certificate actually lives on, where the intersection with the reduced +/// family is far less degenerate — and the same search-then-round-then-check +/// machinery in [`search_rational_family`] is run again on that smaller +/// family, this time with a real chance of landing exactly on it. If the +/// guess is wrong, `restrict_by_near_null_vectors` (or the reduced family's +/// own exact check) simply comes back empty and the next guess is tried — +/// every path through this function is either an exact rational identity or +/// nothing, never a fabricated certificate. +#[allow(clippy::too_many_arguments)] +fn facial_reduction_search( + nvars: usize, + basis: &[Exponents], + target: &RatPoly, + base_rat: &[Vec], + dirs_rat: &[Vec>], + candidates: &[Vec>], +) -> Option { + let max_corank = *FACIAL_CORANK_GUESSES.iter().max().unwrap_or(&0); + let mut budget = FACIAL_SEARCH_BUDGET; + for q_best in candidates.iter().take(FACIAL_CANDIDATES) { + let near_null = smallest_magnitude_eigenvectors(q_best, max_corank); + for &corank in FACIAL_CORANK_GUESSES { + if budget == 0 { + return None; + } + if corank > near_null.len() { + continue; + } + for &max_den in FACIAL_DENOM_CAPS { + let vs: Option>> = near_null[..corank] + .iter() + .map(|v| { + v.iter() + .map(|&x| round_to_rational(x, max_den)) + .collect::>>() + }) + .collect(); + let Some(vs) = vs else { continue }; + let Some((new_base, new_dirs)) = + restrict_by_near_null_vectors(base_rat, dirs_rat, &vs) + else { + continue; + }; + if new_dirs.len() >= dirs_rat.len() { + // The guess added no real constraint (rare, but possible + // if the rounded vector happens to be in the common + // nullspace of every direction) — searching the + // "reduced" family again would just repeat work. + continue; + } + if budget == 0 { + return None; + } + budget -= 1; + let (found, _deeper) = + search_rational_family(nvars, basis, target, &new_base, &new_dirs); + if found.is_some() { + return found; + } } } } None } +pub fn psd_search(target: &RatPoly, basis_deg: u32) -> Option { + let nvars = target.nvars(); + // A homogeneous target of degree exactly `2·basis_deg` needs only the + // monomials of degree *exactly* `basis_deg` in its Gram basis — mixing in + // lower-degree monomials can only ever contribute to coefficients the + // target does not have, since every product of two basis monomials of + // unequal degree still sums to `2·basis_deg` only when *both* already + // have degree `basis_deg`. This is standard (Blekherman–Parrilo–Thomas, + // Prop. 3.29): a homogeneous SOS decomposition can always be taken with + // homogeneous summands. Restricting here is not just an optimisation — + // the search is numeric, and a smaller basis is the difference between + // "converges" and "not within budget" on cases like Motzkin. + let basis: Vec = match target.is_homogeneous() { + Some(d) if d == 2 * basis_deg => monomial_basis(nvars, basis_deg) + .into_iter() + .filter(|e| e.iter().sum::() == basis_deg) + .collect(), + _ => monomial_basis(nvars, basis_deg), + }; + let n = basis.len(); + if n == 0 { + return if target.is_zero() { + Some(SosPoly::default()) + } else { + None + }; + } + + let (rows, rhs) = gram_system(target, &basis); + let sol = solve_affine(&rows, &rhs)?; + + let base_rat = unpack(n, &sol.particular); + let dirs_rat: Vec>> = sol.nullspace.iter().map(|d| unpack(n, d)).collect(); + + let (found, candidates) = search_rational_family(nvars, &basis, target, &base_rat, &dirs_rat); + if found.is_some() { + return found; + } + if candidates.is_empty() { + return None; + } + facial_reduction_search(nvars, &basis, target, &base_rat, &dirs_rat, &candidates) +} + #[cfg(test)] mod tests { use super::*; @@ -560,6 +864,41 @@ mod tests { than leaving it as a smoke test" ); } + + #[test] + fn psd_search_certifies_robinsons_form_with_a_reznick_multiplier() { + // Robinson's form: x^6+y^6+z^6 - (x^4y^2+x^2y^4+y^4z^2+y^2z^4+x^4z^2+x^2z^4) + 3x^2y^2z^2. + // A second textbook PSD-not-SOS example (distinct from Motzkin) whose + // multiplier certificate the Douglas-Rachford / facial-reduction + // search below now reaches. Direct search (no multiplier) still + // fails, as it must: Robinson's form is genuinely not SOS itself. + let mono = |ex: Vec, c: i64| RatPoly::monomial(3, ex, Rational::from(c)); + let mut r = mono(vec![6, 0, 0], 1); + r = r.add(&mono(vec![0, 6, 0], 1)); + r = r.add(&mono(vec![0, 0, 6], 1)); + r = r.add(&mono(vec![4, 2, 0], -1)); + r = r.add(&mono(vec![2, 4, 0], -1)); + r = r.add(&mono(vec![0, 4, 2], -1)); + r = r.add(&mono(vec![0, 2, 4], -1)); + r = r.add(&mono(vec![4, 0, 2], -1)); + r = r.add(&mono(vec![2, 0, 4], -1)); + r = r.add(&mono(vec![2, 2, 2], 3)); + assert_eq!(r.is_homogeneous(), Some(6)); + + assert!( + psd_search(&r, 3).is_none(), + "Robinson's form is not itself SOS, so a direct (unmultiplied) search must refuse" + ); + + let sigma = RatPoly::sum_of_squares(3); + let q = r.mul(&sigma); + assert_eq!(q.is_homogeneous(), Some(8)); + let sos = psd_search(&q, 4) + .expect("(x^2+y^2+z^2)*Robinson is a classical SOS example; the search should find it"); + // Exact re-expansion, over the rationals — the actual soundness + // argument, not the numeric search that proposed it. + assert_eq!(sos.to_poly(3), q); + } } #[cfg(test)] diff --git a/alkahest-core/src/real/sos/sdp.rs b/alkahest-core/src/real/sos/sdp.rs index 977ce299..ca77192d 100644 --- a/alkahest-core/src/real/sos/sdp.rs +++ b/alkahest-core/src/real/sos/sdp.rs @@ -110,6 +110,30 @@ pub fn min_eigenvalue(q: &[Vec]) -> f64 { vals.into_iter().fold(f64::INFINITY, f64::min) } +/// The `k` eigenvectors of `q` (symmetric) whose eigenvalues have the +/// smallest *magnitude*, nearest-to-zero first. +/// +/// Used by `super::psd::facial_reduction_search` to guess the nullspace of +/// a nearby singular PSD matrix from a numeric near-feasible point: when a +/// witnessing Gram matrix is singular (the boundary-only intersections that +/// make plain alternating projection or Douglas–Rachford converge only +/// asymptotically — see `Family::douglas_rachford_from`'s doc comment), the +/// eigenvectors of a near-optimal numeric point with the smallest-magnitude +/// eigenvalues are a good numeric proxy for the true certificate's exact +/// null directions, even while the numeric search has not yet closed the +/// last stretch to an exact zero eigenvalue. +pub fn smallest_magnitude_eigenvectors(q: &[Vec], k: usize) -> Vec> { + let (vals, vecs) = jacobi_eigen(q); + let mut idx: Vec = (0..vals.len()).collect(); + idx.sort_by(|&a, &b| { + vals[a] + .abs() + .partial_cmp(&vals[b].abs()) + .unwrap_or(std::cmp::Ordering::Equal) + }); + idx.into_iter().take(k).map(|i| vecs[i].clone()).collect() +} + fn dot(a: &[Vec], b: &[Vec]) -> f64 { a.iter() .zip(b.iter()) @@ -117,6 +141,19 @@ fn dot(a: &[Vec], b: &[Vec]) -> f64 { .sum() } +/// Elementwise `alpha·a + beta·b` for two matrices of the same shape. +fn mat_combo(alpha: f64, a: &[Vec], beta: f64, b: &[Vec]) -> Vec> { + a.iter() + .zip(b.iter()) + .map(|(ra, rb)| { + ra.iter() + .zip(rb.iter()) + .map(|(x, y)| alpha * x + beta * y) + .collect() + }) + .collect() +} + /// Solve a small symmetric positive definite system by Cholesky. fn solve_spd(g: &[Vec], rhs: &[f64]) -> Option> { let n = g.len(); @@ -255,6 +292,89 @@ impl Family { } Some(t) } + + /// Douglas–Rachford splitting between `{Q : Q ⪰ floor·I}` and this + /// affine family, with over-relaxation parameter `lambda`. + /// + /// [`Self::search_from`] (plain alternating projection) only ever + /// *projects* onto each set in turn, and is well known to converge + /// arbitrarily slowly — stalling asymptotically close to, but short of, + /// the true intersection — when the two sets meet **tangentially** + /// (non-transversally): exactly the situation when the witnessing Gram + /// matrix is singular, i.e. sits on the relative boundary of the PSD + /// cone rather than its interior. That is the textbook case for tight + /// SOS certificates such as Motzkin's, and is what + /// `super::psd::diag::diag_step1_step2_trajectory_and_family_sanity` + /// shows happening in practice (min eigenvalue creeping from about −1.6 + /// to about −0.0018 and no further as the floor is annealed to 0). + /// + /// Douglas–Rachford instead *reflects* through each set + /// (`R = 2·P − id`) rather than merely projecting, and combines the two + /// reflections: + /// + /// ```text + /// Q_{k+1} = Q_k + λ·(P_psd(R_family(Q_k)) − P_family(Q_k)) + /// ``` + /// + /// which is the standard Lions–Mercier (1979) splitting; `λ = 1` is + /// plain (non-relaxed) Douglas–Rachford, and `λ ∈ (1, 2)` is + /// over-relaxation, which gives the iterate room to overshoot past a + /// shallow tangential approach instead of creeping toward it — the + /// standard reason Douglas–Rachford (over-relaxed or not) is the usual + /// upgrade path from alternating projection for exactly this failure + /// mode in the convex-feasibility and semidefinite-programming + /// literature (see e.g. Bauschke & Combettes on reflection methods, or + /// Douglas–Rachford applied directly to semidefinite feasibility + /// problems). + /// + /// Unlike [`Self::search_from`], the iterate here is an *ambient* + /// symmetric matrix, not necessarily a point of the family — both + /// `P_family` (via `parameters_of`) and `P_psd` (via `project_psd`) + /// accept any symmetric matrix as input, which is what + /// lets the reflections be taken at all. The parameter vector returned + /// is `P_family` of the final iterate. As with [`Self::search_from`], + /// convergence is not guaranteed and not required: the result is a + /// **suggestion**, checked exactly by the caller, never trusted here. + pub fn douglas_rachford_from( + &self, + start: Vec, + floor: f64, + lambda: f64, + iters: usize, + ) -> Option> { + let mut q = self.at(&start); + for _ in 0..iters { + let pa_t = self.parameters_of(&q)?; + let pa = self.at(&pa_t); + // Reflect through the affine family: R_family(q) = 2·pa − q. + let r_family = mat_combo(2.0, &pa, -1.0, &q); + let pb = project_psd(&r_family, floor); + // q_{k+1} = q + λ·(pb − pa). + let step = mat_combo(1.0, &pb, -1.0, &pa); + let q_next = mat_combo(1.0, &q, lambda, &step); + + // Track movement of the *iterate itself* (`q`), not the shadow + // sequence `P_family(q)`: the shadow can go numerically static + // for a stretch (or even cycle) while `q` is still moving, which + // would make an early exit keyed on the shadow declare + // convergence falsely. + let moved: f64 = q_next + .iter() + .zip(q.iter()) + .flat_map(|(ra, rb)| ra.iter().zip(rb.iter())) + .map(|(a, b)| (a - b).abs()) + .fold(0.0, f64::max); + q = q_next; + if moved < 1e-13 { + break; + } + } + let final_t = self.parameters_of(&q)?; + if final_t.iter().any(|v| !v.is_finite()) { + return None; + } + Some(final_t) + } } #[cfg(test)] diff --git a/alkahest-skill/alkahest.md b/alkahest-skill/alkahest.md index 5dba36f9..328ae2c2 100644 --- a/alkahest-skill/alkahest.md +++ b/alkahest-skill/alkahest.md @@ -1437,5 +1437,5 @@ reg.coverage_report_markdown() # same, rendered as a Markdown table 29. **`cert.specialize_at_root_of_unity(d, n)` is the decision that carries a `q_zeilberger` verdict to `q = ζ_d`, and it is three-valued** (since 3.9). A proved `Q(q)` recurrence does not by itself license setting `q` to a primitive `d`-th root of unity — a coefficient or a sum value can have a pole there, and specialising anyway is the `q`-analogue of the A279013 failure (item 22): a certificate that re-checks perfectly while the specialised claim is false. The hypotheses (no pole in any `a_i(qⁿ)` or `S(n+i)` at `ζ_d`) are decided **exactly**, by polynomial divisibility by `Φ_d(q)` over `Q` in the cyclotomic field `Q(ζ_d) = Q[q]/(Φ_d(q))` — never numerically — and `cyclotomic_polynomial(pool, d)` exposes `Φ_d(q)` itself so a caller can redo the check by hand. `status` is `"specializes"` (proved, and re-checked as an exact identity in `Q(ζ_d)` before being returned), `"obstructed"` (a pole was **exhibited** — `sum_value`/`coefficient` raise, but `sum_valuation(i)` is still available since the negative valuation *is* the obstruction — and this is not a claim the specialised identity is false, only that this route is blocked), or `"unknown"` (the generic boundary verdict was already `"unknown"`, so there is nothing to specialise). Three things a `"specializes"` verdict does **not** by itself mean, each with its own accessor: `is_vacuous` (every coefficient died — always true at `d = 1`, the `q → 1` limit — so the recurrence is `0 = 0`, still true, but empty), `leading_coefficient_survives` (`False` means the specialised recurrence no longer determines the last value from the earlier ones), and `support_shrinks` (`q`-Lucas killing terms — `[2;1]_q = 1 + q` is non-zero in `Q(q)` and zero at `ζ_2` — reported via `effective_support`, which can shrink but never grow). `sum_valuation(i)` is the `q`-supercongruence content itself: the exact integer `v` with `Φ_d(q)^v ∥ S(n+i)`, so `v ≥ r` is precisely `Φ_d(q)^r | S(n)`. -30. **`sos_decompose` now tries the full PSD Gram cone and a Reznick multiplier search before refusing, but still cannot certify the hardest classical boundary cases** (since 3.9). Past diagonal dominance (`E-SOS-002` from DSOS alone) it searches the general PSD Gram cone, and past that — when `p` itself is not SOS — tries `(x_1²+…+x_n²)^N·p` for `N = 1..4` and searches *that* cone; a witness for `p < 0` still refuses separately with `E-SOS-003`, unaffected. Every certificate this returns is exact end to end: the numeric search only ever proposes a Gram matrix, which is rounded to nearby rationals and re-expanded to check it equals the target exactly before anything is returned — a `Some`/returned certificate is always sound regardless of what the float search converged to. Budget exhaustion is still `E-SOS-002`, undecided, never "not SOS" — say so, don't paraphrase it as a disproof. **What it does not yet do:** Motzkin's polynomial and Robinson's form — the textbook PSD-not-SOS examples whose multiplier certificates are *singular* Gram matrices sitting exactly on the boundary of the PSD cone — are not found by the current search; it's a diagnosed convergence limitation of the annealed alternating-projection method on tangential intersections (verified sound on a planted boundary example of the same size), not a soundness bug, and the tests record `undecided` rather than a false certificate. Raise `basis_degree`, or fall back to `alkahest.decide`, exactly as for any other `E-SOS-002`. +30. **`sos_decompose` tries the full PSD Gram cone and a Reznick multiplier search before refusing, and now certifies Motzkin and Robinson's form too** (since 3.9). Past diagonal dominance (`E-SOS-002` from DSOS alone) it searches the general PSD Gram cone, and past that — when `p` itself is not SOS — tries `(x_1²+…+x_n²)^N·p` for `N = 1..4` and searches *that* cone; a witness for `p < 0` still refuses separately with `E-SOS-003`, unaffected. Every certificate this returns is exact end to end: the numeric search only ever proposes a Gram matrix, which is rounded to nearby rationals and re-expanded to check it equals the target exactly before anything is returned — a `Some`/returned certificate is always sound regardless of what the float search converged to. Budget exhaustion is still `E-SOS-002`, undecided, never "not SOS" — say so, don't paraphrase it as a disproof. **The textbook PSD-not-SOS examples whose multiplier certificates are *singular* Gram matrices sitting exactly on the boundary of the PSD cone** — Motzkin's polynomial and Robinson's form — used to be out of reach for the original annealed alternating-projection search (a diagnosed convergence limitation at tangential PSD-cone intersections, not a soundness bug); the search now also tries Douglas–Rachford splitting with over-relaxation and a facial-reduction step, and with them both examples are found and exactly re-verified. **What's still open:** the homogeneous 3-variable form of Motzkin (a larger nullspace than the affine 2-variable case) is still not reached, so a boundary-only certificate is not guaranteed to be found in general — `E-SOS-002` still means "not found within this search", never "not SOS". Raise `basis_degree`, or fall back to `alkahest.decide`, exactly as for any other `E-SOS-002`. 31. **Double sums need `experimental.telescope2d`, not `zeilberger`** (since 3.9). `zeilberger`/`q_zeilberger` reach a sum over *one* index; `telescope2d(term, n, j, k)` is the Apagodu–Zeilberger generalization to a proper hypergeometric `F(n,j,k)` with **two** bound indices `j`, `k`, returning `a_0(n), …, a_J(n)` and *two* certificates `cert1`, `cert2` with `Σ_i a_i(n)·F(n+i,j,k) = Δ_j(cert1·F) + Δ_k(cert2·F)`, re-checked exactly in `Q(n,j,k)`. Three real, stated scope limits, not unfinished polish: (1) the certificate ansatz uses a *fixed* denominator built from `F`'s own shift-ratio denominators rather than a minimal 2-D Gosper normal form, so a search that finds nothing raises `E-HOLO-041` and does not prove no certificate exists; (2) `cert.boundary_status(j_lo, j_hi, k_lo, k_hi)` only accepts **constant** (not `n`-dependent) rectangles — for a natural range like `j = 0..n`, pick a fixed bound safely larger than any `n` you check and let `F`'s own combinatorial vanishing do the rest, exactly as the module's own worked example does; (3) the boundary of a rectangle is **four one-dimensional strip sums along its edges, not four corner-point evaluations** — a naive corner-evaluation formula is simply wrong — and this version only proves the sufficient (not necessary) condition that each strip vanishes identically, so `boundary_status` can return `"unknown"` for a boundary that is genuinely `0` but not by that pointwise route; it never guesses `"vanishes"`. There is no inhomogeneous `"nonzero"` verdict yet — an unresolved strip is always `"unknown"`. `E-HOLO-040` is the class refusal (not proper hypergeometric in `n, j, k`), `E-HOLO-042` a malformed call (`n`, `j`, `k` not distinct). diff --git a/docs/mdbook/src/positivity.md b/docs/mdbook/src/positivity.md index e8e87bc1..d9558ba0 100644 --- a/docs/mdbook/src/positivity.md +++ b/docs/mdbook/src/positivity.md @@ -57,12 +57,20 @@ wrong by collapsing: > and it carries a witness point. `E-SOS-002` is *not* a claim that the polynomial is not a sum of squares, and -certainly not that it is negative. The canonical example is the **Motzkin -polynomial** `x⁴y² + x²y⁴ − 3x²y² + 1`, which is non-negative everywhere but -provably not a sum of squares. Asked to decompose it, this module refuses with -`E-SOS-002` — it does not report it as negative, and it does not invent a -decomposition. Choi–Lam and Robinson refuse the same way, and for the same -reason: **the refusal is a property of the search, not of the polynomial.** +certainly not that it is negative. The canonical illustration is the +**Motzkin polynomial** `x⁴y² + x²y⁴ − 3x²y² + 1`, which is non-negative +everywhere but provably not a sum of squares *itself* — asked to decompose it +*directly* (no multiplier), this module refuses with `E-SOS-002`, correctly: +it does not report it as negative, and it does not invent a decomposition. +`sos_decompose`'s full pipeline does not stop there, though (see "What the +search actually covers" below) — it also tries multiplying by a power of +`x²+y²` before giving up, and that succeeds for Motzkin, so the *end-to-end* +call returns a certificate, not a refusal. The homogeneous 3-variable form of +Motzkin still refuses even through the full pipeline (multiplier search +included) — that is the actual reachable illustration of a genuine +`E-SOS-002` from this module today: **the refusal is a property of the +search, not of the polynomial**, and which polynomials it applies to shifts +as the search grows more complete. The three-way branch a loop should write: @@ -115,19 +123,28 @@ The search tries three things, in order, before refusing: `E-SOS-002` at the end of all three is phrased as a statement about the search, not the polynomial — the search's incompleteness, at any step. -**A known, diagnosed gap in step 3:** Motzkin's polynomial and Robinson's -form — the textbook PSD-not-SOS examples — are not yet certified even with a -multiplier, because their witnessing Gram matrices are *singular*, sitting -exactly on the boundary of the PSD cone rather than its interior. The -annealed search converges toward that boundary (monotonically, confirmed by a -diagnostic trajectory) but does not reliably close the last, asymptotically -slow stretch — the textbook behaviour of alternating projection at a -tangential (non-transversal) intersection. This was checked to be a -convergence limitation and not a bug in the search machinery: an independent -sanity check confirms the affine Gram-matrix family is constructed correctly, -and a synthetic planted example with a singular Gram matrix of the same size -*is* found and exactly re-verified. The tests for Motzkin record `undecided` -rather than a false certificate. +**Step 3's search has to work harder than plain alternating projection**, +because the multiplier certificates it exists for are frequently *tight* — +Motzkin's polynomial and Robinson's form (the textbook PSD-not-SOS examples) +both have witnessing Gram matrices that are *singular*, sitting exactly on +the boundary of the PSD cone rather than its interior. A first version of +this search (annealed alternating projection with several random restarts) +converged toward that boundary monotonically (confirmed by a diagnostic +trajectory) but never reliably closed the last, asymptotically slow stretch — +the textbook behaviour of alternating projection at a tangential +(non-transversal) set intersection. The search now also tries +Douglas–Rachford splitting with over-relaxation and a facial-reduction step — +both standard escapes for exactly this stall — and with them, both +`(x²+y²)·Motzkin(x,y)` and `(x²+y²+z²)·Robinson(x,y,z)` are found and +exactly re-verified. **What's still open:** the *homogeneous 3-variable* +form of Motzkin, `(x²+y²+z²)·(x⁴y²+x²y⁴−3x²y²z²+z⁶)`, is not — a larger +nullspace than the affine 2-variable case, and still hard enough for the +current search as tuned. This was checked to be a genuine, scoped search +limitation and not a bug in the machinery: an independent sanity check +confirms the affine Gram-matrix family is constructed correctly, and a +synthetic planted example with a singular Gram matrix of the same size *is* +found and exactly re-verified. The remaining test for the 3-variable +Motzkin form records `undecided` rather than a false certificate. ## Constrained certificates @@ -228,16 +245,18 @@ statement is good enough. Shipped: exact rational SOS over the DSOS generator cone, a general PSD Gram search (floating-point proposal, exact rational verification) for cases DSOS -alone refuses, a Reznick multiplier search (`(Σxᵢ²)^N·p` for `N ≤ 4`) on top -of that, Handelman certificates on basic semialgebraic sets, exact -verification, and Lean export. - -Not yet shipped: reliable certification of the hardest classical -boundary-case examples (Motzkin, Robinson — see above; the multiplier search -finds them for neither, diagnosed as an alternating-projection convergence -limitation at a tangential PSD-cone intersection, not a soundness gap), a -proper interior-point or facial-reduction-based solver that would close that -gap, and Putinar-style certificates with genuine SOS — rather than -non-negative constant — multipliers on the *constraints*. -`CertificateKind::Putinar` exists in the certificate type so those can be -added without a shape change. +alone refuses — with Douglas–Rachford splitting and a facial-reduction step +alongside the original annealed alternating projection, specifically so +boundary-only (singular Gram matrix) certificates are reachable — a Reznick +multiplier search (`(Σxᵢ²)^N·p` for `N ≤ 4`) on top of that (finds both +Motzkin's polynomial and Robinson's form), Handelman certificates on basic +semialgebraic sets, exact verification, and Lean export. + +Not yet shipped: reliable certification of *every* boundary-case example — +the homogeneous 3-variable form of Motzkin specifically is still out of +reach (see above), so this is a real but narrower gap than "Motzkin doesn't +certify" was in the prior release — a proper interior-point solver that +would close it more systematically, and Putinar-style certificates with +genuine SOS — rather than non-negative constant — multipliers on the +*constraints*. `CertificateKind::Putinar` exists in the certificate type so +those can be added without a shape change. diff --git a/tests/test_sos.py b/tests/test_sos.py index 9de7208d..d2ce6269 100644 --- a/tests/test_sos.py +++ b/tests/test_sos.py @@ -102,12 +102,12 @@ def test_negative_polynomial_yields_a_witness_not_a_shrug(): assert "< 0" in str(excinfo.value) -def test_motzkin_refuses_without_claiming_negativity(): - """The Motzkin polynomial is non-negative but not a sum of squares. - - The one thing the implementation must not do is call it negative. It has to - come back with E-SOS-002 — "no certificate of this shape" — which is an - honest statement about the search, not about the polynomial. +def test_motzkin_certifies_via_a_multiplier(): + """The Motzkin polynomial is non-negative but not a sum of squares + *itself* — the textbook example (Hilbert 1888) of that phenomenon. + Multiplying by (x^2+y^2) is classically known to fix it (Reznick's + theorem), and `sos_decompose` finds that certificate by searching for a + multiplier automatically: the call succeeds, it does not refuse. """ pool = ak.ExprPool() x, y = pool.symbol("x"), pool.symbol("y") @@ -119,12 +119,21 @@ def test_motzkin_refuses_without_claiming_negativity(): + pool.integer(1) ) - with pytest.raises(ak.SosError) as excinfo: - ak.sos_decompose(p, [x, y]) + cert = ak.sos_decompose(p, [x, y]) - assert excinfo.value.code == "E-SOS-002" - # The remediation must tell an agent that this is not a proof of non-SOS. - assert "decide" in excinfo.value.remediation + assert cert.kind == "sos" + # `verify()` re-expands the identity exactly in Q and confirms it holds — + # the real soundness argument, not a numeric search's own confidence. + assert cert.verify() is True + # The identity is `(target) * (multiplier) = rhs` for a multiplier + # certificate, which `identity` renders as two parenthesised factors on + # the left of `=` — distinct from a direct certificate's plain `target = + # rhs`. This is the only Python-level signal that a multiplier was used + # (there is no separate `.multiplier` accessor at this surface). + assert cert.identity.count("=") == 1 + lhs, _rhs = cert.identity.split("=", 1) + assert lhs.strip().startswith("(") + assert lhs.count(")") >= 2 def test_non_polynomial_is_refused():