|
| 1 | +//! W4 — PSD at depth-infinity (census M-3). |
| 2 | +//! |
| 3 | +//! ndarray's `prove_pillar_11` certifies the TRUNCATED (d=2, deg-3) kernel's |
| 4 | +//! stability. The property kernel machines actually rely on — PSD-ness of the |
| 5 | +//! Gram — is uncertified for the depth-infinity kernel everywhere. This leg |
| 6 | +//! re-runs that machinery over `sigker::signature_kernel_pde` values on the |
| 7 | +//! SAME Brownian pool, per ruling Q2 (cross-repo call, not a second f32 |
| 8 | +//! Goursat port; a port waits for W5's trigger). |
| 9 | +//! |
| 10 | +//! # Strengthened relative to the truncated battery |
| 11 | +//! |
| 12 | +//! `prove_pillar_11`'s "PSD" criteria are diagonal positivity plus |
| 13 | +//! Cauchy-Schwarz. Both are NECESSARY conditions, neither is sufficient: a |
| 14 | +//! matrix can satisfy them and still have a negative eigenvalue. This leg |
| 15 | +//! adds the sufficient test — a Cholesky factorization, which exists iff the |
| 16 | +//! matrix is positive definite — and pairs it with the falsifier that the |
| 17 | +//! weaker criteria cannot supply (see `a_non_psd_gram_is_rejected`). |
| 18 | +//! |
| 19 | +//! # Tolerance |
| 20 | +//! |
| 21 | +//! The Gram is built in f64 from a first-order solver, so the diagonal is not |
| 22 | +//! exact; Cholesky is run with a relative jitter of `JITTER` times the mean |
| 23 | +//! diagonal, the standard numerical allowance. The falsifier below confirms |
| 24 | +//! that this jitter does not swallow a genuinely indefinite matrix. |
| 25 | +
|
| 26 | +use ndarray::hpc::pillar::signature::brownian_path_d2; |
| 27 | +use ndarray::hpc::pillar::SplitMix64; |
| 28 | +use sigker::signature_kernel_pde; |
| 29 | + |
| 30 | +/// Same seed as the truncated battery, so both certify the same pool. |
| 31 | +const PILLAR_11_SEED: u64 = 0x5EED_1111_5164_A7AB; |
| 32 | +/// Gram pool — O(N^2) kernel solves, so kept small. |
| 33 | +const N_PATHS: usize = 64; |
| 34 | +/// Concentration pool — O(N) solves, matched to the truncated battery's 1000 |
| 35 | +/// so the two numbers are directly comparable (see that leg's doc). |
| 36 | +const N_CONC: usize = 1000; |
| 37 | +const N_STEPS: usize = 50; |
| 38 | +const JITTER: f64 = 1e-9; |
| 39 | + |
| 40 | +/// The f32 hardware path buffer, in sigker's point-list shape. |
| 41 | +fn as_points(flat: &[f32], n_points: usize) -> Vec<Vec<f64>> { |
| 42 | + (0..n_points) |
| 43 | + .map(|k| vec![flat[2 * k] as f64, flat[2 * k + 1] as f64]) |
| 44 | + .collect() |
| 45 | +} |
| 46 | + |
| 47 | +fn brownian_pool(seed: u64, n_paths: usize) -> Vec<Vec<Vec<f64>>> { |
| 48 | + let mut rng = SplitMix64::new(seed); |
| 49 | + (0..n_paths) |
| 50 | + .map(|_| { |
| 51 | + let p = brownian_path_d2(&mut rng, N_STEPS); |
| 52 | + as_points(&p, N_STEPS + 1) |
| 53 | + }) |
| 54 | + .collect() |
| 55 | +} |
| 56 | + |
| 57 | +fn gram(pool: &[Vec<Vec<f64>>]) -> Vec<Vec<f64>> { |
| 58 | + let n = pool.len(); |
| 59 | + let mut k = vec![vec![0.0f64; n]; n]; |
| 60 | + for i in 0..n { |
| 61 | + for j in i..n { |
| 62 | + let v = signature_kernel_pde(&pool[i], &pool[j]); |
| 63 | + k[i][j] = v; |
| 64 | + k[j][i] = v; // the kernel is symmetric by construction |
| 65 | + } |
| 66 | + } |
| 67 | + k |
| 68 | +} |
| 69 | + |
| 70 | +/// Cholesky with relative jitter. `Some(_)` iff the matrix is positive |
| 71 | +/// definite to that tolerance — the SUFFICIENT PSD test the truncated |
| 72 | +/// battery lacks. |
| 73 | +fn cholesky(k: &[Vec<f64>], jitter: f64) -> Option<usize> { |
| 74 | + let n = k.len(); |
| 75 | + let mean_diag = (0..n).map(|i| k[i][i]).sum::<f64>() / n as f64; |
| 76 | + let eps = jitter * mean_diag.abs().max(1.0); |
| 77 | + let mut l = vec![vec![0.0f64; n]; n]; |
| 78 | + for i in 0..n { |
| 79 | + for j in 0..=i { |
| 80 | + let mut s = k[i][j]; |
| 81 | + if i == j { |
| 82 | + s += eps; |
| 83 | + } |
| 84 | + // Dot of the two already-computed row prefixes. Both borrows are |
| 85 | + // immutable and end before the write below. |
| 86 | + let dot: f64 = l[i][..j].iter().zip(&l[j][..j]).map(|(a, b)| a * b).sum(); |
| 87 | + s -= dot; |
| 88 | + if i == j { |
| 89 | + if s <= 0.0 { |
| 90 | + return Some(i); // failed at this leading minor |
| 91 | + } |
| 92 | + l[i][i] = s.sqrt(); |
| 93 | + } else { |
| 94 | + l[i][j] = s / l[j][j]; |
| 95 | + } |
| 96 | + } |
| 97 | + } |
| 98 | + None |
| 99 | +} |
| 100 | + |
| 101 | +#[test] |
| 102 | +fn the_depth_infinity_gram_is_positive_definite() { |
| 103 | + let pool = brownian_pool(PILLAR_11_SEED, N_PATHS); |
| 104 | + let k = gram(&pool); |
| 105 | + |
| 106 | + // Necessary conditions first — the same two the truncated battery uses, |
| 107 | + // so a failure here is directly comparable to that report. |
| 108 | + let mut cs_violations = 0usize; |
| 109 | + for i in 0..N_PATHS { |
| 110 | + assert!(k[i][i] > 0.0, "depth-inf self-kernel K[{i},{i}] = {} is not positive", k[i][i]); |
| 111 | + for j in i + 1..N_PATHS { |
| 112 | + if k[i][j] * k[i][j] > k[i][i] * k[j][j] * 1.001 { |
| 113 | + cs_violations += 1; |
| 114 | + } |
| 115 | + } |
| 116 | + } |
| 117 | + assert_eq!(cs_violations, 0, "Cauchy-Schwarz violated at depth-infinity"); |
| 118 | + |
| 119 | + // The sufficient one. |
| 120 | + assert!( |
| 121 | + cholesky(&k, JITTER).is_none(), |
| 122 | + "depth-infinity Gram is NOT positive definite (Cholesky failed at \ |
| 123 | + leading minor {:?}) — M-3 would be falsified", |
| 124 | + cholesky(&k, JITTER) |
| 125 | + ); |
| 126 | + println!("W4: depth-inf Gram over {N_PATHS} Brownian paths — diag > 0, Cauchy-Schwarz clean, Cholesky OK"); |
| 127 | +} |
| 128 | + |
| 129 | +/// Concentration is measured at `N_CONC`, NOT at `N_PATHS`. |
| 130 | +/// |
| 131 | +/// The 0.20 bound comes from the truncated battery, which runs 1000 paths. |
| 132 | +/// Applying it to the 64-path Gram pool would silently be a DIFFERENT gate — |
| 133 | +/// measured, at N = 64 the truncated kernel itself concentrates to 0.3461 and |
| 134 | +/// would fail its own bound. Half-mean agreement is a sample-size statistic |
| 135 | +/// before it is a kernel property, so the comparison is only meaningful at |
| 136 | +/// matched N (`examples/w4_concentration_sweep.rs`): |
| 137 | +/// |
| 138 | +/// ```text |
| 139 | +/// depth-INFINITY depth-3 TRUNCATED |
| 140 | +/// N concentr N concentr |
| 141 | +/// 64 0.2892 64 0.3461 |
| 142 | +/// 128 0.2156 1000 0.0481 |
| 143 | +/// 256 0.1884 |
| 144 | +/// 512 0.0053 |
| 145 | +/// 1000 0.0038 |
| 146 | +/// ``` |
| 147 | +/// |
| 148 | +/// At matched N = 1000 the depth-infinity kernel concentrates BETTER than the |
| 149 | +/// truncated one it extends (0.0038 vs 0.0481). The Gram/Cholesky leg above |
| 150 | +/// stays at 64 paths because it costs O(N^2) solves; this leg is O(N). |
| 151 | +#[test] |
| 152 | +fn the_depth_infinity_self_kernel_concentrates() { |
| 153 | + let pool = brownian_pool(PILLAR_11_SEED, N_CONC); |
| 154 | + let self_k: Vec<f64> = pool.iter().map(|p| signature_kernel_pde(p, p)).collect(); |
| 155 | + let concentration = half_mean_gap(&self_k); |
| 156 | + println!( |
| 157 | + "W4: depth-inf self-kernel concentration {concentration:.4} over {N_CONC} paths \ |
| 158 | + (bound 0.20, the truncated battery's; truncated scores 0.0481 at the same N)" |
| 159 | + ); |
| 160 | + assert!( |
| 161 | + concentration < 0.20, |
| 162 | + "depth-inf self-kernel half-means disagree by {concentration:.3} at N = {N_CONC}" |
| 163 | + ); |
| 164 | +} |
| 165 | + |
| 166 | +/// Can-stay-silent's partner: a pool that genuinely does NOT concentrate must |
| 167 | +/// be caught at the SAME N, so the green above is a property of the kernel |
| 168 | +/// and not of the bound being loose. The fixture rescales the second half of |
| 169 | +/// the pool, which is exactly the drift a half-mean statistic exists to see. |
| 170 | +#[test] |
| 171 | +fn a_non_concentrating_pool_is_caught() { |
| 172 | + let pool = brownian_pool(PILLAR_11_SEED, N_CONC); |
| 173 | + let self_k: Vec<f64> = pool |
| 174 | + .iter() |
| 175 | + .enumerate() |
| 176 | + .map(|(i, p)| { |
| 177 | + if i < N_CONC / 2 { |
| 178 | + signature_kernel_pde(p, p) |
| 179 | + } else { |
| 180 | + let scaled: Vec<Vec<f64>> = p |
| 181 | + .iter() |
| 182 | + .map(|pt| pt.iter().map(|v| v * 1.35).collect()) |
| 183 | + .collect(); |
| 184 | + signature_kernel_pde(&scaled, &scaled) |
| 185 | + } |
| 186 | + }) |
| 187 | + .collect(); |
| 188 | + let concentration = half_mean_gap(&self_k); |
| 189 | + assert!( |
| 190 | + concentration >= 0.20, |
| 191 | + "ANTI-VACUITY FAILED: a pool whose second half is rescaled 1.35x still \ |
| 192 | + concentrated to {concentration:.4} — the bound cannot see drift" |
| 193 | + ); |
| 194 | + println!("W4 anti-vacuity: rescaled-half pool reaches concentration {concentration:.4} (bound 0.20)"); |
| 195 | +} |
| 196 | + |
| 197 | +/// Relative gap between the two half-means — the truncated battery's statistic. |
| 198 | +fn half_mean_gap(v: &[f64]) -> f64 { |
| 199 | + let h = v.len() / 2; |
| 200 | + let m1 = v[..h].iter().sum::<f64>() / h as f64; |
| 201 | + let m2 = v[h..].iter().sum::<f64>() / (v.len() - h) as f64; |
| 202 | + let combined = v.iter().sum::<f64>() / v.len() as f64; |
| 203 | + if combined > 0.0 { |
| 204 | + (m1 - m2).abs() / combined |
| 205 | + } else { |
| 206 | + f64::INFINITY |
| 207 | + } |
| 208 | +} |
| 209 | + |
| 210 | +// ════════════════════════════════════════════════════════════════════════════ |
| 211 | +// Falsifier pair for the PSD gate. |
| 212 | +// ════════════════════════════════════════════════════════════════════════════ |
| 213 | + |
| 214 | +/// Can-fire: an indefinite matrix of the same shape and scale must be |
| 215 | +/// REJECTED. Without this, `cholesky` returning `None` would be evidence of |
| 216 | +/// nothing — and it is precisely the test the truncated battery cannot run, |
| 217 | +/// because diagonal-positivity + Cauchy-Schwarz both HOLD on this matrix. |
| 218 | +#[test] |
| 219 | +fn a_non_psd_gram_is_rejected() { |
| 220 | + let pool = brownian_pool(PILLAR_11_SEED, 8); |
| 221 | + let mut k = gram(&pool); |
| 222 | + // Flip the sign of one off-diagonal pair. Diagonals are untouched (still |
| 223 | + // positive) and the magnitudes are unchanged (so Cauchy-Schwarz still |
| 224 | + // holds) — only definiteness is destroyed. |
| 225 | + let scale = (k[0][0] * k[1][1]).sqrt(); |
| 226 | + k[0][1] = -0.999 * scale; |
| 227 | + k[1][0] = -0.999 * scale; |
| 228 | + k[0][2] = 0.999 * (k[0][0] * k[2][2]).sqrt(); |
| 229 | + k[2][0] = k[0][2]; |
| 230 | + k[1][2] = 0.999 * (k[1][1] * k[2][2]).sqrt(); |
| 231 | + k[2][1] = k[1][2]; |
| 232 | + |
| 233 | + // The weak criteria still pass — that is the point. |
| 234 | + for i in 0..3 { |
| 235 | + assert!(k[i][i] > 0.0, "diagonal must stay positive in the fixture"); |
| 236 | + for j in i + 1..3 { |
| 237 | + assert!(k[i][j] * k[i][j] <= k[i][i] * k[j][j] * 1.001, "Cauchy-Schwarz must still hold in the fixture"); |
| 238 | + } |
| 239 | + } |
| 240 | + // The strong one must catch it. |
| 241 | + assert!(cholesky(&k, JITTER).is_some(), "ANTI-VACUITY FAILED: an indefinite Gram passed the Cholesky gate"); |
| 242 | + println!( |
| 243 | + "W4 anti-vacuity: indefinite Gram passes diag>0 AND Cauchy-Schwarz, rejected by Cholesky at minor {:?}", |
| 244 | + cholesky(&k, JITTER) |
| 245 | + ); |
| 246 | +} |
0 commit comments