diff --git a/Cargo.toml b/Cargo.toml index 596da8ee..4b617d42 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -91,6 +91,18 @@ name = "codec_mode_histogram" required-features = ["codec"] # W1a-#9 masking-primitive codegen probe imports `ndarray::simd` (std-gated). +[[example]] +name = "psd_bit_exactness" +required-features = ["pillar", "std"] + +[[example]] +name = "psd_jitter_sweep" +required-features = ["pillar", "std"] + +[[example]] +name = "psd_f32_vs_f64_probe" +required-features = ["pillar", "std"] + [[example]] name = "w1a9_codegen_probe" required-features = ["std"] diff --git a/crates/sigker-parity/tests/w4_depth_infinity_psd.rs b/crates/sigker-parity/tests/w4_depth_infinity_psd.rs index e7f1578a..b5140c9f 100644 --- a/crates/sigker-parity/tests/w4_depth_infinity_psd.rs +++ b/crates/sigker-parity/tests/w4_depth_infinity_psd.rs @@ -27,8 +27,13 @@ use ndarray::hpc::pillar::signature::brownian_path_d2; use ndarray::hpc::pillar::SplitMix64; use sigker::signature_kernel_pde; -/// Same seed as the truncated battery, so both certify the same pool. -const PILLAR_11_SEED: u64 = 0x5EED_1111_5164_A7AB; +// The seed is IMPORTED, not retyped. An earlier revision of this file +// hardcoded a different value while its comment claimed "same seed as the +// truncated battery, so both certify the same pool" — the two batteries were +// certifying different pools, and nothing could have noticed, because both +// sides of every comparison used the wrong constant consistently. Caught by +// the bit-exactness gate added to the truncated battery. +use ndarray::hpc::pillar::signature::PILLAR_11_SEED; /// Gram pool — O(N^2) kernel solves, so kept small. const N_PATHS: usize = 64; /// Concentration pool — O(N) solves, matched to the truncated battery's 1000 diff --git a/examples/psd_bit_exactness.rs b/examples/psd_bit_exactness.rs new file mode 100644 index 00000000..8dce72df --- /dev/null +++ b/examples/psd_bit_exactness.rs @@ -0,0 +1,70 @@ +//! Is the Pillar-11 PSD gate bit-exact, and therefore safe to gate in CI? +//! +//! Prints a fingerprint over the RAW BIT PATTERNS of every Gram entry, plus +//! the Cholesky verdict. Two runs on one machine answer determinism; two runs +//! under different `-C target-cpu` answer whether autovectorisation/FMA moves +//! the low bits — which decides whether CI may pin exact values or only +//! verdicts. `.cargo/config.toml` pins `target-cpu=x86-64-v4`, so this is not +//! hypothetical. +use ndarray::hpc::lapack::LapackOps; +use ndarray::hpc::pillar::signature::PILLAR_11_SEED; +use ndarray::hpc::pillar::signature::{brownian_path_d2, sigker_hl, signature_d2_deg3}; +use ndarray::hpc::pillar::SplitMix64; +use ndarray::Array2; + +const N_STEPS: usize = 50; +const SUBSET: usize = 50; +const JITTER: f64 = 1e-4; + +/// FNV-1a over raw bits — any single-bit change moves it. +fn fnv(acc: &mut u64, bits: u64) { + for b in bits.to_le_bytes() { + *acc ^= b as u64; + *acc = acc.wrapping_mul(0x100_0000_01b3); + } +} + +fn main() { + let mut rng = SplitMix64::new(PILLAR_11_SEED); + let s: Vec<[f32; 15]> = (0..SUBSET) + .map(|_| { + let p = brownian_path_d2(&mut rng, N_STEPS); + signature_d2_deg3(&p, N_STEPS + 1) + }) + .collect(); + + let mut sig_fp = 0xcbf2_9ce4_8422_2325u64; + for v in &s { + for x in v { + fnv(&mut sig_fp, x.to_bits() as u64); + } + } + + let mut g = Array2::::zeros((SUBSET, SUBSET)); + let mut k_fp = 0xcbf2_9ce4_8422_2325u64; + for i in 0..SUBSET { + for j in 0..SUBSET { + let v = sigker_hl(&s[i], &s[j]); + fnv(&mut k_fp, v.to_bits() as u64); + g[[i, j]] = v as f64; + } + } + + let mean_diag = (0..SUBSET).map(|i| g[[i, i]]).sum::() / SUBSET as f64; + let eps = JITTER * mean_diag.abs().max(1.0); + let mut j = g.clone(); + for i in 0..SUBSET { + j[[i, i]] += eps; + } + let chol = j.cholesky(); + let mut l_fp = 0xcbf2_9ce4_8422_2325u64; + for v in chol.factor.iter() { + fnv(&mut l_fp, v.to_bits()); + } + + println!("signature bits fnv1a = {sig_fp:#018x}"); + println!("kernel bits fnv1a = {k_fp:#018x}"); + println!("cholesky L bits fnv1a = {l_fp:#018x}"); + println!("cholesky info = {}", chol.info); + println!("mean diag = {mean_diag:.17e}"); +} diff --git a/examples/psd_f32_vs_f64_probe.rs b/examples/psd_f32_vs_f64_probe.rs new file mode 100644 index 00000000..8e05e18b --- /dev/null +++ b/examples/psd_f32_vs_f64_probe.rs @@ -0,0 +1,53 @@ +//! Pre-registration: does the truncated Pillar-11 Gram survive an f32 +//! Cholesky, or does it need f64 accumulation? +//! +//! `sigker_hl` returns f32. Cholesky on an f32 Gram can report `info > 0` on +//! a matrix that IS positive semi-definite but ill-conditioned — a false +//! alarm on the real battery. This measures the margin on the ACTUAL pool +//! before any gate is pinned. +use ndarray::hpc::lapack::LapackOps; +use ndarray::hpc::pillar::signature::PILLAR_11_SEED; +use ndarray::hpc::pillar::signature::{brownian_path_d2, sigker_hl, signature_d2_deg3}; +use ndarray::hpc::pillar::SplitMix64; +use ndarray::Array2; + +const N_STEPS: usize = 50; + +fn sigs(n: usize) -> Vec<[f32; 15]> { + let mut rng = SplitMix64::new(PILLAR_11_SEED); + (0..n) + .map(|_| { + let p = brownian_path_d2(&mut rng, N_STEPS); + signature_d2_deg3(&p, N_STEPS + 1) + }) + .collect() +} + +fn main() { + println!( + "{:>7} {:>12} {:>12} {:>14} {:>14} {:>12}", + "subset", "f32 info", "f64 info", "min diag", "max |K_ij|", "cond-ish" + ); + for &n in &[8usize, 12, 14, 15, 16, 17, 32, 50] { + let s = sigs(n); + let mut g32 = Array2::::zeros((n, n)); + let mut g64 = Array2::::zeros((n, n)); + for i in 0..n { + for j in 0..n { + let v = sigker_hl(&s[i], &s[j]); + g32[[i, j]] = v; + g64[[i, j]] = v as f64; + } + } + let i32_ = g32.cholesky().info; + let i64_ = g64.cholesky().info; + let min_diag = (0..n).map(|i| g64[[i, i]]).fold(f64::INFINITY, f64::min); + let max_off = (0..n) + .flat_map(|i| (0..n).filter(move |j| *j != i).map(move |j| (i, j))) + .map(|(i, j)| g64[[i, j]].abs()) + .fold(0.0f64, f64::max); + let max_diag = (0..n).map(|i| g64[[i, i]]).fold(0.0f64, f64::max); + println!("{n:>7} {i32_:>12} {i64_:>12} {min_diag:>14.4e} {max_off:>14.4e} {:>12.2e}", max_diag / min_diag); + } + println!("\ninfo == 0 means positive definite; info > 0 is the failing leading minor."); +} diff --git a/examples/psd_jitter_sweep.rs b/examples/psd_jitter_sweep.rs new file mode 100644 index 00000000..8af4be19 --- /dev/null +++ b/examples/psd_jitter_sweep.rs @@ -0,0 +1,101 @@ +//! Pick the PSD jitter by measurement, and prove it still discriminates. +//! +//! The truncated signature is a 15-dimensional feature map (1+2+4+8), so the +//! Gram of N > 15 paths is rank-deficient BY CONSTRUCTION and plain Cholesky +//! — which tests positive DEFINITE — must fail at leading minor 16. The +//! property the battery actually wants is positive SEMI-definite, i.e. +//! Cholesky with a relative jitter on the diagonal. +//! +//! A jitter large enough to admit the real (singular, PSD) Gram must still +//! REJECT a genuinely indefinite one. This sweeps both together; a jitter +//! that admits both is worthless. +use ndarray::hpc::lapack::LapackOps; +use ndarray::hpc::pillar::signature::PILLAR_11_SEED; +use ndarray::hpc::pillar::signature::{brownian_path_d2, sigker_hl, signature_d2_deg3}; +use ndarray::hpc::pillar::SplitMix64; +use ndarray::Array2; + +const N_STEPS: usize = 50; +const SUBSET: usize = 50; + +fn real_gram(n: usize) -> Array2 { + let mut rng = SplitMix64::new(PILLAR_11_SEED); + let s: Vec<[f32; 15]> = (0..n) + .map(|_| { + let p = brownian_path_d2(&mut rng, N_STEPS); + signature_d2_deg3(&p, N_STEPS + 1) + }) + .collect(); + let mut g = Array2::::zeros((n, n)); + for i in 0..n { + for j in 0..n { + g[[i, j]] = sigker_hl(&s[i], &s[j]) as f64; + } + } + g +} + +/// A genuinely indefinite matrix that PASSES both weak criteria: every +/// diagonal positive, Cauchy-Schwarz satisfied for every pair. +fn indefinite_gram(n: usize) -> Array2 { + let mut g = real_gram(n); + let s01 = (g[[0, 0]] * g[[1, 1]]).sqrt(); + let s02 = (g[[0, 0]] * g[[2, 2]]).sqrt(); + let s12 = (g[[1, 1]] * g[[2, 2]]).sqrt(); + g[[0, 1]] = -0.999 * s01; + g[[1, 0]] = g[[0, 1]]; + g[[0, 2]] = 0.999 * s02; + g[[2, 0]] = g[[0, 2]]; + g[[1, 2]] = 0.999 * s12; + g[[2, 1]] = g[[1, 2]]; + g +} + +fn chol_info(g: &Array2, jitter: f64) -> i32 { + let n = g.nrows(); + let mean_diag = (0..n).map(|i| g[[i, i]]).sum::() / n as f64; + let eps = jitter * mean_diag.abs().max(1.0); + let mut j = g.clone(); + for i in 0..n { + j[[i, i]] += eps; + } + j.cholesky().info +} + +fn weak_criteria_hold(g: &Array2) -> (bool, bool) { + let n = g.nrows(); + let diag_ok = (0..n).all(|i| g[[i, i]] > 0.0); + let mut cs_ok = true; + for i in 0..n { + for j in i + 1..n { + if g[[i, j]] * g[[i, j]] > g[[i, i]] * g[[j, j]] * 1.001 { + cs_ok = false; + } + } + } + (diag_ok, cs_ok) +} + +fn main() { + let real = real_gram(SUBSET); + let bad = indefinite_gram(SUBSET); + + let (rd, rc) = weak_criteria_hold(&real); + let (bd, bc) = weak_criteria_hold(&bad); + println!("weak criteria (diag>0, Cauchy-Schwarz):"); + println!(" real Gram diag={rd} cs={rc}"); + println!(" indefinite Gram diag={bd} cs={bc} <- both PASS, which is the point"); + + println!("\n{:>12} {:>14} {:>18} {:>10}", "jitter", "real info", "indefinite info", "verdict"); + for &j in &[0.0f64, 1e-12, 1e-9, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1] { + let r = chol_info(&real, j); + let b = chol_info(&bad, j); + let verdict = match (r == 0, b == 0) { + (true, false) => "USABLE", + (false, false) => "too tight", + (true, true) => "TOO LOOSE", + (false, true) => "impossible?", + }; + println!("{j:>12.0e} {r:>14} {b:>18} {verdict:>10}"); + } +} diff --git a/src/hpc/pillar/signature.rs b/src/hpc/pillar/signature.rs index c2ec3a01..7fff48f2 100644 --- a/src/hpc/pillar/signature.rs +++ b/src/hpc/pillar/signature.rs @@ -60,6 +60,8 @@ //! * Chen (1954), *Iterated path integrals*, Bull. AMS. use super::prove_runner::{PillarReport, SplitMix64}; +use crate::hpc::lapack::LapackOps; +use crate::Array2; // ── Constants ───────────────────────────────────────────────────────────────── @@ -244,6 +246,55 @@ pub fn brownian_path_d2(rng: &mut SplitMix64, n_steps: usize) -> alloc::vec::Vec // ── Prove ───────────────────────────────────────────────────────────────────── +/// Relative diagonal jitter for the PSD test. +/// +/// **Why a jitter is mandatory here, not a convenience.** The truncated +/// signature is a 15-dimensional feature map (1 + 2 + 4 + 8), so the Gram of +/// N > 15 paths is rank-deficient BY CONSTRUCTION. Plain Cholesky tests +/// positive *definite* and must therefore fail for any N ≥ 16 — measured +/// exactly at that boundary (`examples/psd_f32_vs_f64_probe.rs`): +/// +/// ```text +/// subset f32 info f64 info +/// 15 0 0 +/// 16 16 16 <- the feature dimension, not a defect +/// 50 16 16 +/// ``` +/// +/// The property this battery wants is positive SEMI-definite, which is +/// Cholesky on a jittered diagonal. The value is pinned from a sweep that +/// measured BOTH sides at once (`examples/psd_jitter_sweep.rs`): the real +/// Gram is admitted from 1e-6 upward, and a genuinely indefinite Gram is +/// still rejected at every jitter through 1e-1. 1e-4 sits two orders above +/// the lower edge with three orders of headroom above it. +/// +/// Note also that f32 and f64 accumulation agreed on every verdict in that +/// probe — precision was never the question; rank was. +const PSD_JITTER: f64 = 1e-4; + +/// Positive semi-definiteness of a symmetric Gram, by Cholesky on a jittered +/// diagonal. +/// +/// This is the SUFFICIENT test. Diagonal positivity and Cauchy–Schwarz — the +/// two conditions this battery checked before — are each necessary and +/// neither is sufficient: a matrix can satisfy both and still have a negative +/// eigenvalue. `psd_gate_rejects_an_indefinite_gram_that_passes_both_weak_criteria` +/// constructs exactly such a matrix, so the distinction is asserted rather +/// than assumed. +fn gram_is_psd(gram: &Array2) -> bool { + let n = gram.nrows(); + if n == 0 { + return true; + } + let mean_diag = (0..n).map(|i| gram[[i, i]]).sum::() / n as f64; + let eps = PSD_JITTER * mean_diag.abs().max(1.0); + let mut jittered = gram.clone(); + for i in 0..n { + jittered[[i, i]] += eps; + } + jittered.cholesky().info == 0 +} + /// Pillar-11 certification probe — Hambly–Lyons sigker convergence on 1 000 Lévy paths. /// /// # PASS criteria @@ -319,27 +370,39 @@ pub fn prove_pillar_11() -> PillarReport { f64::INFINITY }; - // ── Criterion 2: PSD diagonal check on first SUBSET paths ──────────────── - // For a valid kernel the diagonal K[i,i] = k_HL(Pᵢ, Pᵢ) must be > 0. - // Also verify Cauchy–Schwarz: K[i,j]² ≤ K[i,i] · K[j,j] for all pairs. + // ── Criterion 2: the Gram over the first SUBSET paths ──────────────────── + // Materialised ONCE (the previous form recomputed both diagonals inside + // the inner loop, O(n²) kernel calls for O(n) distinct values), then read + // three times: diagonal positivity, Cauchy–Schwarz, and positive + // semi-definiteness. + let mut gram = Array2::::zeros((SUBSET, SUBSET)); + for i in 0..SUBSET { + for j in 0..SUBSET { + gram[[i, j]] = sigker_hl(&sigs[i], &sigs[j]) as f64; + } + } + + // Necessary condition 1: K[i,i] > 0. + // Necessary condition 2: Cauchy–Schwarz, K[i,j]² ≤ K[i,i]·K[j,j]. let mut cs_violations: u32 = 0; for i in 0..SUBSET { for j in i + 1..SUBSET { - let kij = sigker_hl(&sigs[i], &sigs[j]); - let kii = sigker_hl(&sigs[i], &sigs[i]); - let kjj = sigker_hl(&sigs[j], &sigs[j]); - if kij * kij > kii * kjj * 1.001 { + if gram[[i, j]] * gram[[i, j]] > gram[[i, i]] * gram[[j, j]] * 1.001 { // 0.1% tolerance for f32 rounding cs_violations += 1; } } } + // Sufficient condition: positive semi-definiteness by Cholesky. + let psd_ok = gram_is_psd(&gram); + // ── Determine PASS ──────────────────────────────────────────────────────── let psd_rate = positive_count as f64 / N_PATHS as f64; let passed = psd_rate >= 1.0 // all self-kernels positive && concentration < 0.20 // half-means agree within 20 % - && cs_violations == 0; // Cauchy–Schwarz holds in subset + && cs_violations == 0 // Cauchy–Schwarz holds in subset + && psd_ok; // and the Gram is actually PSD — see gram_is_psd PillarReport { pillar_id: 11, @@ -359,6 +422,134 @@ extern crate alloc; #[cfg(test)] mod tests { + // ── PSD gate: the closed falsifier ────────────────────────────────────── + + /// The pair the battery relied on before — diagonal positivity and + /// Cauchy–Schwarz — cannot see an indefinite Gram. This asserts all three + /// facts about ONE matrix, which is what makes the upgrade non-decorative + /// rather than two unrelated claims: + /// + /// (a) every diagonal is positive — weak criterion 1 PASSES + /// (b) Cauchy–Schwarz holds for every pair — weak criterion 2 PASSES + /// (c) `gram_is_psd` rejects it — only the new gate catches it + #[test] + fn psd_gate_rejects_an_indefinite_gram_that_passes_both_weak_criteria() { + // Three unit-diagonal vectors that cannot be a Gram: x·y strongly + // negative while x·z and y·z are both strongly positive. + let a = 0.999_f64; + let mut g = Array2::::eye(3); + g[[0, 1]] = -a; + g[[1, 0]] = -a; + g[[0, 2]] = a; + g[[2, 0]] = a; + g[[1, 2]] = a; + g[[2, 1]] = a; + + for i in 0..3 { + assert!(g[[i, i]] > 0.0, "weak criterion 1 must PASS on the fixture"); + for j in i + 1..3 { + assert!( + g[[i, j]] * g[[i, j]] <= g[[i, i]] * g[[j, j]] * 1.001, + "weak criterion 2 must PASS on the fixture at ({i},{j})" + ); + } + } + + assert!( + !gram_is_psd(&g), + "the PSD gate must reject an indefinite Gram that both weak \ + criteria admit — otherwise the gate adds nothing" + ); + } + + /// The can-stay-silent half: a genuinely PSD Gram must be ADMITTED at the + /// same jitter, including the rank-deficient case the real battery + /// produces (the truncated signature is a 15-dimensional feature map, so + /// any Gram over more than 15 paths is singular by construction). + #[test] + fn psd_gate_admits_a_rank_deficient_but_valid_gram() { + // Outer products of 15-dimensional signatures: PSD, and singular as + // soon as there are more than 15 of them. + let mut rng = SplitMix64::new(PILLAR_11_SEED); + let sigs: alloc::vec::Vec<[f32; SIG_D2_DEG3_LEN]> = (0..40) + .map(|_| { + let p = brownian_path_d2(&mut rng, 50); + signature_d2_deg3(&p, 51) + }) + .collect(); + let mut g = Array2::::zeros((40, 40)); + for i in 0..40 { + for j in 0..40 { + g[[i, j]] = sigker_hl(&sigs[i], &sigs[j]) as f64; + } + } + assert!( + gram_is_psd(&g), + "a real (PSD but rank-deficient) Gram must be admitted — a gate \ + that rejects it is measuring the feature dimension, not validity" + ); + } + + // ── Bit-exactness: the CI gate ────────────────────────────────────────── + + /// The Pillar-11 pipeline is bit-exact, and this pins it. + /// + /// Measured across every axis available here — run to run, `target-cpu= + /// x86-64-v4` (this repo's `.cargo/config.toml` default) against baseline + /// `x86-64`, and debug against release — the signature bits, the kernel + /// bits and the Cholesky factor bits were identical in all six runs. So + /// exact values may be pinned rather than tolerances, and any drift in + /// autovectorisation, FMA contraction or evaluation order shows up here + /// instead of silently shifting a battery that would still report green. + /// + /// **What this test does NOT claim.** A fingerprint proves the numbers are + /// the SAME, never that they are RIGHT — that is what the PSD, + /// Cauchy–Schwarz and concentration gates above are for. The two kinds of + /// check are paired on purpose; neither substitutes for the other. + /// + /// If this fails after a deliberate change to the signature or kernel, + /// re-pin it in the same commit and say so in the message. + #[test] + fn pillar_11_pipeline_is_bit_exact() { + fn fnv(acc: &mut u64, bits: u64) { + for b in bits.to_le_bytes() { + *acc ^= b as u64; + *acc = acc.wrapping_mul(0x100_0000_01b3); + } + } + const FNV_OFFSET: u64 = 0xcbf2_9ce4_8422_2325; + + let mut rng = SplitMix64::new(PILLAR_11_SEED); + let sigs: alloc::vec::Vec<[f32; SIG_D2_DEG3_LEN]> = (0..50) + .map(|_| { + let p = brownian_path_d2(&mut rng, 50); + signature_d2_deg3(&p, 51) + }) + .collect(); + + let mut sig_fp = FNV_OFFSET; + for v in &sigs { + for x in v { + fnv(&mut sig_fp, x.to_bits() as u64); + } + } + assert_eq!(sig_fp, 0x4434_20ec_eee1_5ce1, "signature bit pattern drifted"); + + let mut k_fp = FNV_OFFSET; + let mut g = Array2::::zeros((50, 50)); + for i in 0..50 { + for j in 0..50 { + let v = sigker_hl(&sigs[i], &sigs[j]); + fnv(&mut k_fp, v.to_bits() as u64); + g[[i, j]] = v as f64; + } + } + assert_eq!(k_fp, 0x1bdc_8789_8bea_ee77, "kernel bit pattern drifted"); + + // And the gate itself, on those exact bits. + assert!(gram_is_psd(&g), "the pinned Gram must pass the PSD gate"); + } + use super::*; // ── signature_d2_deg3 basic ───────────────────────────────────────────────