WP28: pamle library (models, splits, Newton variants, CSL/AEE/Cox-Snell), symbolic checks, 29 tests, literature audits A-B
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
20 changed files +1,292 −0
added
literature/PRIOR_ART_MATRIX.csv
+13 −0
@@ -0,0 +1,13 @@ | ||
| 1 | +citation_key,authors,year,title,venue,doi_or_url,verified_full_text,assumptions,number_of_rounds,upstream_bytes,local_computation,statistical_rate,proof_type,applicable_models,overlap_with_WP28,novelty_threat,notes,page_numbers_checked | |
| 2 | +zhang2013communication,"Zhang, Yuchen; Duchi, John C.; Wainwright, Martin J.",2013,"Communication-Efficient Algorithms for Statistical Optimization","Journal of Machine Learning Research 14:3321-3363",https://www.jmlr.org/papers/v14/zhang13b.html (arXiv:1209.4129),"yes (official JMLR PDF)","A1 compact convex Theta with theta* interior; A2 local strong convexity of population risk at theta*; A3 8th moments of gradient and Hessian deviation + Hessian L(x)-Lipschitz with 8th moments; A4 (SAVGM only) third derivatives Lipschitz with E[M^8]<inf","1","AVGM: one d-vector per machine; SAVGM: corrected d-vector (two local fits)","one ERM on n points (+ one on rn points for SAVGM)","Thm 1 (AVGM): MSE = O(1/N + K^2/N^2); optimal rate iff K <~ sqrt(N) (eq.10). Thm 4 (SAVGM, (theta1 - r theta2)/(1-r)): MSE <= ((2+3r)/(1-r)^2)(1/N)E||H^-1 grad||^2 + c n^-3/(r(1-r)^2) + O((1-r)^-2 m^-1 n^-2); optimal rate iff K <~ N^(2/3) (eq.11)","non-asymptotic MSE bounds via Taylor expansion of FOC + 8th-moment concentration","general smooth M-estimation / ERM, fixed d","(i) partial: SAVGM is a single-ratio subsampling bias correction with inflated leading variance constant, not a jackknife; (ii) none; (iii) none","medium for (i); low for (ii)/(iii)","No Richardson/jackknife; multiple subsamples left open (p.3330); Newton only as local solver","3321-3332, 3337, 3340-3345" | |
| 3 | +huang2019distributed,"Huang, Cheng; Huo, Xiaoming",2019,"A distributed one-step estimator","Mathematical Programming 174(1-2):41-76",https://doi.org/10.1007/s10107-019-01369-0 (arXiv:1511.01443v2),"arXiv v2 yes; Springer text NOT VERIFIED","compact convex Theta; population Hessian negative definite at theta0; 8th moments; Hessian Lipschitz; concave criterion; iid equal split; fixed d","2","round 1: d-vector; round 2: d gradient + d^2 Hessian at the broadcast average","one local M-estimation; one gradient+Hessian at theta^(0); centre one d x d solve","Thm 3.4: one-step from the plain average is sqrt(N)-normal and equivalent to the centralized M-estimator when k = O(sqrt N). Thm 3.5: MSE <= 2Tr(Sigma)/N + O(N^-2) + O(k^4 N^-4). Lemma B.5: local bias C/n","asymptotic normality + non-asymptotic MSE bound","M-estimators and MLE, fixed d","(ii) direct: one global Newton step already proved equivalent, but only under k = O(sqrt N) and from an uncorrected average; (i) none; (iii) none (q=1 only)","high for (ii) unless positioned as bias-corrected input + enlarged frontier","Hessian = average of local Hessians at the averaged estimate; sharper k=o(N^(3/4)) implied by k^4N^-4 not stated by the authors","arXiv v2 pp. 1-3, 5-8, 11, 15-16, 24-26" | |
| 4 | +lin2011aggregated,"Lin, Nan; Xi, Ruibin",2011,"Aggregated estimating equation estimation","Statistics and Its Interface 4(1):73-83",https://doi.org/10.4310/SII.2011.v4.n1.a8,"yes (publisher reprint)","(C1)-(C6): smooth psi, p.d. Jacobians, strong consistency, bounded n^-1 A_k, bounded second derivatives, tail bound with alpha in (1/4,1/2); fixed p","1","p-vector + p x p local Jacobian at the LOCAL estimate per block","one EE solve per block; centre one p x p solve","AEE = (sum A_k)^-1 sum A_k beta_k. Thm 4.1: ||AEE - EE_full|| <= C(||worst local error||^2 + ...). Thm 4.2: sqrt(N)||AEE - EE_full|| -> 0 if K = O(n^gamma), gamma < 1/3 (K << N^(1/4))","asymptotic (Taylor expansion of aggregated score + tail bounds)","estimating equations, GLM/QLE, fixed p","none for (i)-(iii); Hessian-weighted one-shot baseline with restrictive K growth","low","No variance estimator; OLAP motivation","73-76, 78-79" | |
| 5 | +rosenblatt2016optimality,"Rosenblatt, Jonathan D.; Nadler, Boaz",2016,"On the optimality of averaging in distributed statistical learning","Information and Inference: A Journal of the IMA 5(4):379-404",https://doi.org/10.1093/imaiai/iaw013 (arXiv:1407.2724v2),"arXiv v2 yes; OUP pagination NOT VERIFIED","Set 1: van der Vaart 5.23-type; Set 2: Rilstone-Srivastava-Ullah derivatives to order 3/4; Set 3: Gaussian design p/n -> kappa; m FIXED in fixed-p results","1","one p-vector per machine","one local M-estimation","Thm 1 first-order equivalence (m fixed); Thm 2 B2(theta_bar) = delta/n = m B2(theta_hat_N); Thm 3 second-order MSE; HD Lemma 1: MSE ratio 1 + kappa(r2/r1)(1-1/m)","asymptotic stochastic expansions","M-estimators; Gaussian linear models (HD)","(i) gives the target of the correction (bias = m x centralized bias) and warns about variance; no correction proposed; (ii)/(iii) none","low (supporting)","Remark 2 lists bias-reduction routes","arXiv v2 pp. 1-12, 15-18" | |
| 6 | +jordan2019csl,"Jordan, Michael I.; Lee, Jason D.; Yang, Yun",2019,Communication-Efficient Distributed Statistical Inference,"JASA 114(526):668-681",10.1080/01621459.2018.1429274,"full text arXiv 1605.07689v3 (JASA typeset 403; published theorem numbers NOT VERIFIED)","PA compact convex Theta; PB mu- I <= Hessian(theta*) <= mu+ I; PC identifiability; PD 16th-moment bounds and Lipschitz Hessian; initial estimator O_p(n^-1/2)","1 round per iteration; ceil(log k/log n) iterations; 1 round suffices if ||theta_bar-theta*||=o_P(sqrt(n/N)) (Cor 3.4)","d numbers per machine per round (gradient); no Hessian","machine 1 minimizes the surrogate or takes one quadratic step with its own Hessian","||theta_tilde-theta_hat|| <= C(||theta_bar-theta_hat||^2 + (N^-1/2 + n^-1/2)||theta_bar-theta_hat||) (Thm 3.1, 3.3); contraction C/sqrt(n) per round; translation: K=o(N^1/2) local-MLE start, K=o(N^2/3) average start","finite-sample high-probability bounds + asymptotic normality corollary","regular M-estimators / MLE; l1; quasi-posterior","single-machine Hessian; no bias correction; precedent for the inexact-Hessian factor n^-1/2 (claim iv)",moderate,"JASA numbering unknown","arXiv pp. 7-17" | |
| 7 | +fan2023cease,"Fan, Jianqing; Guo, Yongyi; Wang, Kaizheng",2023,Communication-Efficient Accurate Statistical Estimation,"JASA 118(542):1000-1010",10.1080/01621459.2021.1969238,"full text arXiv 1906.04870v2 + PMC10281708 (numbering identical); JASA typeset 403","rho-strong convexity on a ball; Hessian homogeneity delta; Lipschitz Hessian; GLM canonical link, sub-Gaussian covariates, n >= c p","1 round per iteration (Alg 3) or 2 (Alg 4); 1 step from the plain average suffices iff n^3 >> N(kappa^2 p log N)(p+kappa^2 log p) (p16)","p numbers per machine per round; no Hessian","gradient-enhanced proximal local minimization on machine 1 (Alg 3) or on all machines + averaging (Alg 4)","||theta_t-theta_hat||=O_P(eta^{t/2}||theta_0-theta_hat||), eta=kappa^2 p log N/n (Thm 4.1); linear model: averaged step p^2/n^2 => K=o(N^3/4); single machine p^3/2/n^3/2 => K=o(N^2/3)","deterministic contraction + Hessian concentration","regularized ERM, GLM, linear regression","explicit one-step exponents 2/3 and 3/4 (local / averaged local Hessians) from the plain average; no bias correction; no jackknife comparison","moderate-high for claim (ii) architecture; none for jackknife/7/8","Alg 4 = DANE extension; proximal point view","arXiv pp. 4-18; PMC lines 230-262" | |
| 8 | +chen2022fone,"Chen, Xi; Liu, Weidong; Zhang, Yichen",2022,First-Order Newton-Type Estimator for Distributed Estimation and Inference,"JASA 117(540):1858-1874",10.1080/01621459.2021.1891925,"full text arXiv 1811.11368v2; JASA typeset 403","p=O(n^kappa1); batch >> p log n; local strong convexity at theta*; N=O(n^A); initial O_P(n^-delta1); loss may be non-smooth","K rounds, K constant depending on delta's (Thm 3.3)","O(p) per machine per round; no Hessian ever","machine 1 runs mini-batch first-order iterations approximating the Newton step","DC-SGD needs L=O(sqrt(N/p)) (Thm 3.1, necessary on examples); Dist-FONE: per-round factor n^-delta0 (Thm 3.3)","asymptotic O_P rates by induction over rounds","convex possibly non-smooth losses; diverging p","no bias correction; no K-exponent for fixed rounds; citable sqrt(N) barrier statement",low,"computation on largest machine","arXiv pp. 3, 7-17" | |
| 9 | +shamir2014dane,"Shamir, Ohad; Srebro, Nathan; Zhang, Tong",2014,Communication-Efficient Distributed Optimization using an Approximate Newton-type Method,"Proc. 31st ICML, PMLR 32(2):1000-1008",https://proceedings.mlr.press/v32/shamir14.html,"full text PMLR PDF + arXiv 1312.7853v4","quadratic losses lambda <= Hess <= L; matrix concentration of local Hessians; strongly convex smooth local objectives (Thm 4)","2 averaging communications per iteration; t=O((L/lambda)^2 log(dm)/n log(1/eps)) iterations (Thm 3)","d numbers per machine per averaging; no Hessian","local proximal subproblem with gradient correction","||w^t-w_hat|| <= ||I - H_tilde^-1 H||^t ||w^0-w_hat||, ||I-H_tilde^-1 H|| <= 4 beta^2/lambda^2, beta ~ sqrt(log(dm)/n) (Thm 2, Lemmas 1-2); Thm 1 one-shot averaging fails when lambda <= O(1/sqrt n)","optimization (eps-suboptimality); mirror descent view","smooth strongly convex ERM; quadratics","averaged local Hessians cost ~1/n per round (claim iv precedent); no statistical K-exponent; no bias correction",low,"Thm 5 single-machine variant = CSL precursor; App. A remark on a bias-corrected average","PMLR pp. 1003-1006; arXiv pp. 5-10, 16" | |
| 10 | +wang2017edsl,"Wang, Jialei; Kolar, Mladen; Srebro, Nathan; Zhang, Tong",2017,Efficient Distributed Learning with Sparsity,"Proc. 34th ICML, PMLR 70:3636-3645",https://proceedings.mlr.press/v70/wang17f.html,"full text PMLR PDF + arXiv 1605.07991v1 (numbering differs)","L-smooth loss; restricted strong convexity on machine 1; sub-Gaussian designs; n >~ s^2 log p","1 round per iteration; 1 round matches the lasso iff n >~ m s^2 log p; O(log m) rounds otherwise","p numbers per machine per round","master solves the shifted l1 problem (= subsampled Newton when lambda=0, quadratic)","contraction a_n ~ s sqrt(log p/n) per round (Thm 6); geometric after t rounds (Cor 7)","high-probability recursive bounds","sparse linear/logistic, graphical models","single-machine Hessian; l1; no bias correction; K <~ sqrt(N) for one round",low,"concurrent with CSL","PMLR pp. 3637-3641; arXiv pp. 2-7" | |
| 11 | +battey2018distributed,"Battey, Heather; Fan, Jianqing; Liu, Han; Lu, Junwei; Zhu, Ziwei",2018,Distributed testing and estimation under sparse high dimensional models,"Annals of Statistics 46(3):1352-1382",10.1214/17-AOS1587,"full text arXiv 1509.05457v1 + PMC6051757 (low-dim theorems moved to supplement); journal page mapping NOT VERIFIED","lasso/SCAD/MCP error rates; bounded eigenvalues, sub-Gaussian rows; GLM analogues; finite fourth moments","1 round (one-shot)","d numbers per subsample","local lasso + debiasing (nodewise lassos); low-dim local OLS/GLM","Wald normality k=o((s log d)^-1 sqrt n) (Thm 3.3); lasso rate retained iff k=O(sqrt(n/(s^2 log d))) (Thm 4.3); low-dim GLM k=o(sqrt n/(d v log n)) (arXiv Thm 4.13)","high-probability decomposition + CLT","sparse high-dim linear & GLM; low-dim linear & GLM","no bias correction beyond lasso debiasing; establishes the sqrt(N) baseline (r=1,q=0); Remark 4.5 = same leading term, remainders differ",low,"n = total sample size in their notation","arXiv pp. 1-2, 8, 13-17; PMC lines 92-99, 368-442, 531-537" | |
| 12 | +volgushev2019quantile,"Volgushev, Stanislav; Chao, Shih-Kang; Cheng, Guang",2019,Distributed inference for quantile regression processes,"Annals of Statistics 47(3):1634-1662",10.1214/18-AOS1730,"full text arXiv 1701.06088v3 + NSF-PAR accepted manuscript (numbering identical)","bounded design with eigenvalue bounds; conditional density bounded with bounded derivative; Holder smoothness of tau -> beta(tau)","1 round (one-shot)","m numbers per machine per quantile level","local quantile regression fits","sup_tau ||beta_bar - beta_or|| = O_P(S log N/N + S^1/4 (log N)^7/4/N^3/4) + o_P(N^-1/2) (Thm 3.1); lower bound >= C S/N for some model (Thm 3.2); sufficient S=o(N^1/2/log N), NECESSARY S=o(N^1/2) (Cor 3.3)","new Bahadur representation; explicit worst-case construction","linear quantile regression; series models; conditional distribution functions","template for sharp (necessary + sufficient) partition-count conditions; no correction; no Newton",low,"necessity is worst-case over a model class","arXiv pp. 1, 6, 9-11" | |
| 13 | +shang2017limits,"Shang, Zuofeng; Cheng, Guang",2017,Computational Limits of A Distributed Algorithm for Smoothing Spline,"JMLR 18(108):1-37",https://jmlr.org/papers/v18/16-289.html,"full text JMLR PDF","periodic m-th order smoothing spline; evenly spaced design; iid errors; Sobolev ball","1 round (one-shot averaging)","one function estimate per machine","local smoothing spline with common lambda","minimax optimal iff s=O(N^2m/(2m+1)) (Thm 3.1); worst-case squared bias >= C(a_m n^-2m - 8N^-1) (Thm 3.2) => sharp; testing s=O(N^(4m-1)/(4m+1)) sharp","non-asymptotic bounds in eigen-basis; phase transition","nonparametric regression; testing","sharp maximal number of machines set by local bias vs 1/N variance; same mechanism, nonparametric",low,"improves ZDW bound","JMLR pp. 1, 5-7, 9-10" | |
added
literature/audit_A_zdw_hh_linxi_rn.md
+34 −0
@@ -0,0 +1,34 @@ | ||
| 1 | +# Literature audit A — ZDW 2013, Huang–Huo 2019, Lin–Xi 2011, Rosenblatt–Nadler 2016 (theorem level) | |
| 2 | + | |
| 3 | +Method (verification task, 2026-10-28 18:47–18:53): full-text PDFs read (JMLR official PDF; arXiv 1511.01443v2; arXiv 1407.2724v2; Chalmers-hosted reprint of Lin & Xi with running header "Statistics and Its Interface Volume 4 (2011) 73–83"). Bibliographic records cross-checked with CrossRef. Springer's page for Huang–Huo is bot-blocked: theorem numbers for that paper refer to arXiv v2 (Nov 2015); the published text (36 pp. vs 31 pp.) is NOT VERIFIED. | |
| 4 | + | |
| 5 | +## 1. Zhang, Duchi & Wainwright (2013), JMLR 14:3321–3363 (arXiv 1209.4129), full text read | |
| 6 | +Notation: N = mn samples, m machines (= our K), n = N/m (= our m). θ̄₁ = AVGM (eq. 2, p. 3324). | |
| 7 | +- Assumptions (Sec. 2.1, 3.1, pp. 3324–3326): A1 compact convex Θ, θ* interior; A2 local strong convexity of population risk at θ*; A3 8th moments of gradient and Hessian deviation, Hessian L(x)-Lipschitz with 8th moments; A4 (SAVGM only) third-derivative Lipschitz with E[M⁸] < ∞. Example 1 (p. 3326): without Hessian Lipschitz continuity the local bias is Ω(n^{-1/2}). | |
| 8 | +- Theorem 1 (AVGM, p. 3327): E‖θ̄−θ*‖² ≤ (2/(nm))E‖H^{-1}∇f‖² + (c/(λ²n²))[H² log d + (L²G²/λ²)E‖H^{-1}∇f‖²] + O(m^{-1}n^{-2}) + O(n^{-3}); i.e. O(1/N + K²/N²). Corollary 3 compares with the Hájek–Le Cam bound tr I^{-1}/N. The n^{-2} term is attributed to the O(1/n) local bias (p. 3329). | |
| 9 | +- SAVGM (Sec. 2.2.2, Theorem 4, p. 3329): subsample of size ⌈rn⌉ per machine, θ_SAVGM = (θ̄₁ − rθ̄₂)/(1−r); MSE ≤ ((2+3r)/(1−r)²)(1/(nm))E‖H^{-1}∇f‖² + c[…]n^{-3}/(r(1−r)²) + O((1−r)^{-2}m^{-1}n^{-2}). Leading constant inflated by (2+3r)/(1−r)²; r → 0 recommended (r ∝ m^{3/2}/N). | |
| 10 | +- K-vs-N conditions (p. 3330, eqs. 10–11): AVGM optimal rate iff m ≲ √N; SAVGM iff m ≲ N^{2/3}. These are MSE-rate statements, not first-order equivalence to the full MLE. | |
| 11 | +- No Richardson / higher-order extrapolation, no jackknife (words absent); multiple subsamples left as an open question about variance (p. 3330). No one-step Newton comparison (Newton appears only as the LIBSVM local solver, p. 3337). Section 3.5 is an SGD variant. | |
| 12 | +- Rounds: 1; uplink: one d-vector (two local fits for SAVGM). | |
| 13 | + | |
| 14 | +## 2. Huang & Huo (2019), Math. Programming 174(1–2):41–76, DOI 10.1007/s10107-019-01369-0 (arXiv 1511.01443v2 read) | |
| 15 | +- One-step estimator (Sec. 2.4, eq. 2.9): θ^(0) = simple average of local M-estimators; θ^(1) = θ^(0) − [M̈(θ^(0))]^{-1} Ṁ(θ^(0)) with M the GLOBAL empirical criterion: machines compute local gradient and local Hessian AT THE BROADCAST AVERAGE and send both (d + d² numbers); two rounds. | |
| 16 | +- Assumptions (Sec. 3.1) copied from ZDW (compactness, negative-definite population Hessian, 8th moments, Hessian Lipschitz, concavity); no third-derivative condition. | |
| 17 | +- Theorem 3.4: when k = O(√N), √N(θ^(1) − θ0) ⇒ N(0, Σ): "asymptotically equivalent to the centralized M-estimator". Remark: proof only needs a √N-consistent start. | |
| 18 | +- Theorem 3.5: E‖θ^(1) − θ0‖² ≤ 2 Tr Σ/N + O(N^{-2}) + O(k⁴N^{-4}). Corollary 3.6: MLE efficiency. Lemma B.5: local bias ≤ C̃/n + O(n^{-2}); Lemma C.1 needs k = O(n) = O(√N) for √N-consistency of the average. | |
| 19 | +- No bias correction, no jackknife (SAVGM only as a simulation competitor, Sec. 4). q = 1 only, no repeated rounds. Authors claim priority for the distributed one-step idea. | |
| 20 | +- Note (auditor's derivation, not stated by the authors): the k⁴N^{-4} remainder is o(N^{-1}) whenever k = o(N^{3/4}). | |
| 21 | + | |
| 22 | +## 3. Lin & Xi (2011), Statistics and Its Interface 4(1):73–83, DOI 10.4310/SII.2011.v4.n1.a8, full text read | |
| 23 | +- AEE: A_k = −Σ_i ∂ψ(z_ki, β̂_nk)/∂β (local Jacobian at the LOCAL estimate), β̃ = (Σ A_k)^{-1} Σ A_k β̂_nk (eq. 4). One round, p + p² numbers per block. Motivated by OLAP/data streams. | |
| 24 | +- Theorem 4.1 (K bounded): ‖AEE − EE_full‖ ≤ C₁(‖worst local error‖² + ‖full error‖²). Theorem 4.2: √N‖AEE − EE_full‖ → 0 if K = O(n^γ), 0 < γ < min{1−2α, 4α−1} ≤ 1/3 (⇒ K ≪ N^{1/4}; conversion is the auditor's). Theorem 5.2: GLM/QLE normality under the same condition. | |
| 25 | +- No bias expansion, no correction, no Newton step, no variance estimator. | |
| 26 | + | |
| 27 | +## 4. Rosenblatt & Nadler (2016), Information and Inference 5(4):379–404, DOI 10.1093/imaiai/iaw013 (arXiv 1407.2724v2 read) | |
| 28 | +- Fixed p, m FIXED: Theorem 1 first-order equivalence of the average to the centralized estimator; Theorem 2 (second-order bias, Rilstone–Srivastava–Ullah expansions): B₂(θ̄) = δ/n = m·B₂(θ̂_N); Theorem 3 second-order MSE with excess ((m−1)/m)δδᵀ/n² + …; Remark 2 lists bias-reduction routes (Kim 2006, Liu & Ihler 2014, SAVGM, DDPCA, DANE) and warns bias is traded for variance. High-dimensional regime p/n → κ: Lemma 1 first-order loss linear in m. | |
| 29 | +- No jackknife, no Newton, no K-vs-N growth condition. | |
| 30 | + | |
| 31 | +## Summary of stated K-vs-N conditions | |
| 32 | +ZDW13: √N (AVGM), N^{2/3} (SAVGM), MSE-rate sense. HH19: k = O(√N) for the one-step estimator (Thm 3.4). LX11: K ≪ N^{1/4}. RN16: m fixed. None states a jackknife of local MLEs, a Newton step from a bias-corrected average, an N^{7/8} threshold, or an (r, q) frontier. | |
| 33 | + | |
| 34 | +Local artefacts kept by the verification task: /tmp/wp28lit/{zdw2013,hh2019,rn2016}.{pdf,txt}, /tmp/wp28lit/linxi/. | |
added
literature/audit_B_csl_cease_fone_dane_edsl_battey_volgushev_shang.md
+40 −0
@@ -0,0 +1,40 @@ | ||
| 1 | +# Literature audit B — CSL, CEASE, FONE, DANE, EDSL, Battey et al., Volgushev et al., Shang–Cheng (theorem level) | |
| 2 | + | |
| 3 | +Method (verification task, 2026-10-28 18:47–18:58): full texts read (arXiv PDFs; PMLR PDFs for DANE/EDSL; JMLR PDF for Shang–Cheng; NSF-PAR accepted manuscript for Volgushev et al.; PMC author manuscripts for Battey et al. and Fan–Guo–Wang). Crossref resolved all DOIs. Typeset JASA versions returned 403: published JASA theorem numbers NOT VERIFIED (arXiv / accepted-manuscript numbering used; for CEASE and Battey the accepted manuscripts were checked and match arXiv). Local text copies: /tmp/litaudit/. | |
| 4 | +Notation: K = number of blocks, m = N/K (WP notation) where stated as "translation". | |
| 5 | + | |
| 6 | +## 1. Jordan, Lee & Yang (2019), JASA 114(526):668–681, DOI 10.1080/01621459.2018.1429274 (arXiv 1605.07689v3 read) | |
| 7 | +- CSL surrogate (eq. 8): L̃(θ) = L₁(θ) − ⟨θ, ∇L₁(θ̄) − ∇L_N(θ̄)⟩, machine 1's loss only; all higher-order global derivatives replaced by machine 1's. Approximation error (eq. 6) O(n^{-1/2}‖θ−θ̄‖² + ‖θ−θ̄‖³). | |
| 8 | +- Communication: gradients only, d numbers per machine per round; no Hessians. | |
| 9 | +- Theorem 3.1: ‖θ̃ − θ̂‖ ≤ C₂{‖θ̄−θ̂‖² + (‖θ̂−θ*‖ + |||∇²L₁(θ*) − ∇²L_N(θ*)|||)‖θ̄−θ̂‖} w.p. ≥ 1 − C₁kn^{-8}; i.e. O_p(n^{-1/2})·‖θ̄−θ̂‖ + ‖θ̄−θ̂‖². Theorem 3.3: quadratic one-step θ̃_H = θ̄ − ∇²L₁(θ̄)^{-1}∇L_N(θ̄) obeys the same bound; only an O(n^{-1/2}) start is required. | |
| 10 | +- Multi-round (Algorithm 1): contraction C₃/√n per round; ⌈log k/log n⌉ iterations. | |
| 11 | +- Corollary 3.4: √N-normality if ‖θ̄−θ*‖ = o_P(√(n/N)). Translation: local-MLE start ⇒ K = o(N^{1/2}); plain-average start ⇒ K = o(N^{2/3}) (not stated by the authors). No bias correction. | |
| 12 | + | |
| 13 | +## 2. Fan, Guo & Wang (2023) CEASE, JASA 118(542):1000–1010, DOI 10.1080/01621459.2021.1969238 (arXiv 1906.04870v2 + PMC10281708; numbering identical) | |
| 14 | +- Algorithm 1 = CSL; Algorithm 2 = every machine minimizes its gradient-enhanced loss and the centre averages (DANE-type); Algorithm 3 (CEASE) = Alg. 1 + proximal term; Algorithm 4 = Alg. 2 + proximal term (two rounds per iteration). p numbers per machine per round; no Hessians. | |
| 15 | +- Theorem 3.1 / Corollary 3.2: contraction ≍ δ/ρ (Alg. 3), (δ/ρ)² (Alg. 4) with Hessian homogeneity δ. | |
| 16 | +- Theorem 4.1 (GLM, canonical link, sub-Gaussian covariates, n ≥ c₁p): ‖θ_t − θ̂‖ = O_P(η^{t/2}‖θ₀−θ̂‖), η = κ²(log N)p/n; Alg. 4: η^{t−t₀}. | |
| 17 | +- p. 16 remark (verbatim): starting from the one-shot average, "the one-step estimator θ₁ obtained by Algorithm 3 behaves the same as the global minimizer θ̂ if n³ ≫ N(κ²p log N)(p + κ² log p)"; t steps: n^{t+2} ≫ N(κ²p log N)^t(p + κ² log p). Translation (p, κ fixed, logs dropped): one single-machine-Hessian round from the average ⇒ K = o(N^{2/3}); t rounds ⇒ K = o(N^{(t+1)/(t+2)}). Linear-regression example (pp. 17–18): one averaged local-Hessian step ⇒ ‖θ₁−θ̂‖ = O_P(p²/n²) ⇒ K = o(N^{3/4}); un-averaged ⇒ K = o(N^{2/3}). | |
| 18 | +- Compares with CSL, DANE, classical one-step (Bickel 1975); no SAVGM/jackknife comparison. | |
| 19 | + | |
| 20 | +## 3. Chen, Liu & Zhang (2022) FONE, JASA 117(540):1858–1874, DOI 10.1080/01621459.2021.1891925 (arXiv 1811.11368v2) | |
| 21 | +- Gradients only; machine 1 approximates the Newton step by mini-batch first-order iterations; O(p) per machine per round; no Hessian. | |
| 22 | +- Theorem 3.1 (DC-SGD): plain averaging needs L = O(√(N/p)) machines; "averaging only reduces the variance but not the bias term" (p. 3). Theorem 3.3: per-round improvement n^{-δ₀}; rounds, not K, absorb the partition penalty (N = O(n^A) assumed). | |
| 23 | + | |
| 24 | +## 4. Shamir, Srebro & Zhang (2014) DANE, ICML, PMLR 32(2):1000–1008 (arXiv 1312.7853v4) | |
| 25 | +- Two averaging communications per iteration, d numbers each; no Hessians. Theorem 2 + Lemmas 1–2: contraction ‖I − H̃^{-1}H‖ ≤ 4β²/λ² with β ≍ √(L² log(dm)/n): averaged local Hessians cost ≍ 1/n per round. Theorem 1: one-shot averaging fails when λ ≤ O(1/√n); Appendix A remark that the example extends to a bias-corrected average. Theorem 5: single-machine variant (CSL precursor). | |
| 26 | + | |
| 27 | +## 5. Wang, Kolar, Srebro & Zhang (2017) EDSL, ICML, PMLR 70:3636–3645 (arXiv 1605.07991v1; numbering differs) | |
| 28 | +- CSL surrogate with ℓ₁ penalty, master = machine 1; contraction a_n ≍ s√(log p/n) per round (Thm 6); one round matches the centralized lasso iff n ≳ m s² log p (K ≲ √N up to s, log p); O(log m) rounds otherwise. Concurrent with CSL. | |
| 29 | + | |
| 30 | +## 6. Battey, Fan, Liu, Lu & Zhu (2018), AoS 46(3):1352–1382, DOI 10.1214/17-AOS1587 (arXiv 1509.05457v1 + PMC6051757) | |
| 31 | +- One-shot; high-dim: k = o((s log d)^{-1}√n) for Wald normality (Thm 3.3), lasso rate retained iff k = O(√(n/(s² log d))) (Thm 4.3). Low-dim (supplement): linear k = O(nd/(d ∨ log n)²); GLM k = o(√n/(d ∨ log n)) (√N barrier). Remark 4.5: average and full estimator "share the same leading term in their Taylor expansions; the difference is only the difference of two remainder terms". No bias correction beyond lasso debiasing. | |
| 32 | + | |
| 33 | +## 7. Volgushev, Chao & Cheng (2019), AoS 47(3):1634–1662, DOI 10.1214/18-AOS1730 (arXiv 1701.06088v3 + NSF-PAR) | |
| 34 | +- Theorem 3.1: sup_τ‖β̄ − β̂_or‖ = O_P(S log N/N + S^{1/4}(log N)^{7/4}/N^{3/4}) + o_P(N^{-1/2}). Theorem 3.2 (lower bound): some (P, Z) in the class with ‖β̄ − β̂_or‖ ≥ CS/N with positive probability. Corollary 3.3: sufficient S = o(N^{1/2}/log N), NECESSARY S = o(N^{1/2}). Necessity is worst-case over the model class. Template for sharp partition-count conditions. | |
| 35 | + | |
| 36 | +## 8. Shang & Cheng (2017), JMLR 18(108):1–37 | |
| 37 | +- Theorem 3.1: minimax optimality iff s = O(N^{2m/(2m+1)}); Theorem 3.2: worst-case squared bias ≥ C(a_m n^{-2m} − 8N^{-1}) ⇒ sharp. Mechanism: local bias n^{-2m} vs variance 1/N ("phase transition in the number of machines"). Nonparametric. | |
| 38 | + | |
| 39 | +## Summary of overlap with the WP28 claims (auditor's assessment) | |
| 40 | +(i) Jackknife of local MLEs: in none of the eight; nearest precedent ZDW Theorem 4 (SAVGM) and DANE's remark. (ii) one aggregated-Hessian Newton step from a bias-corrected average with K = o(N^{7/8}): no paper states it; closest Huang–Huo (same architecture, plain average, stated k = O(√N), implicit 3/4), CEASE p. 16 (2/3 single-machine; 3/4 averaged local Hessians in the linear model), JLY Cor. 3.4. (iii) general (r, q) frontier: not stated anywhere; anchor cells (1,0) = 1/2 (Battey, ZDW, FONE Thm 3.1, Volgushev with sharpness), (1,1) = 3/4 (CEASE linear; Huang–Huo bound), single-machine multi-round sequence (t+1)/(t+2) (CEASE). (iv) inexact Hessians: precedents JLY Thm 3.1/3.3 (machine-1 Hessian ⇒ factor n^{-1/2}), DANE Thm 2 (averaged local Hessians ⇒ factor 1/n), CEASE δ vs δ²; none expresses this as a change in the admissible K exponent. Byte accounting: none beyond O(p)/O(d²) remarks. | |
added
pytest.ini
+4 −0
@@ -0,0 +1,4 @@ | ||
| 1 | +[pytest] | |
| 2 | +pythonpath = src | |
| 3 | +testpaths = tests | |
| 4 | +addopts = -q | |
added
requirements.txt
+8 −0
@@ -0,0 +1,8 @@ | ||
| 1 | +numpy==2.4.4 | |
| 2 | +scipy==1.17.1 | |
| 3 | +pandas==3.0.2 | |
| 4 | +matplotlib==3.10.9 | |
| 5 | +sympy==1.14.0 | |
| 6 | +pytest==9.1.1 | |
| 7 | +pyyaml==6.0.3 | |
| 8 | +psutil==7.2.2 | |
modified
research_log.md
+2 −0
@@ -9,3 +9,5 @@ | ||
| 9 | 9 | - 18:47 Four literature-verification tasks launched (ZDW/Huang–Huo/Lin–Xi/Rosenblatt–Nadler; CSL/CEASE/FONE/DANE/EDSL/Battey/ |
| 10 | 10 | Volgushev/Shang–Cheng; jackknife & bias-correction classics; 2013–2026 priority-threat sweep). |
| 11 | 11 | - Decision: identifier WP28 kept (see SOURCE_PROVENANCE.md collision note). |
| 12 | +- 18:55 Audit A received and archived (literature/audit_A_*.md, PRIOR_ART_MATRIX.csv rows 1-4). Key facts: HH19 proves one-step equivalence only under k=O(√N) from an uncorrected average; ZDW SAVGM = single-ratio subsampling correction with inflated variance constant, no jackknife/higher order; LX11 K≪N^{1/4}; RN16 m fixed. | |
| 13 | +- 19:00 Audit B archived (8 papers). Fixed sign of the analytic Cox–Snell correction (b_1 = +½ I^{-1} T; exponential check b_1 = θ0). | |
added
src/pamle/__init__.py
+14 −0
@@ -0,0 +1,14 @@ | ||
| 1 | +"""pamle — reference implementations for UQO WP28: | |
| 2 | +distributed maximum likelihood with local split-sample bias correction and global Newton corrections. | |
| 3 | + | |
| 4 | +Modules | |
| 5 | +------- | |
| 6 | +jackknife Vandermonde extrapolation weights, amplification constants, conditioning. | |
| 7 | +partitioning Block and nested sub-block index construction (equal and unequal sizes). | |
| 8 | +likelihoods Model classes (Exponential, Logistic, Poisson) with batched Newton fitting. | |
| 9 | +algorithms Pilot estimators, exact/inexact Newton corrections, CSL, AEE, analytic bias correction. | |
| 10 | +communication Byte and round accounting for the protocols. | |
| 11 | +diagnostics Monte Carlo summaries, confidence intervals, coverage. | |
| 12 | +protocol Coordinator/worker implementation over TCP for real multi-process / multi-node benchmarks. | |
| 13 | +""" | |
| 14 | +__version__ = "0.1.0" | |
added
src/pamle/algorithms.py
+224 −0
@@ -0,0 +1,224 @@ | ||
| 1 | +"""Distributed estimators (in-memory, statistically faithful to the coordinator/worker protocol). | |
| 2 | + | |
| 3 | +All routines take a batched dataset `db` = {'X': (K, m, p), 'y': (K, m)} (or {'x': (K, m)}) so that block k is | |
| 4 | +`db[...][k]`. Worker-side quantities are computed block by block; coordinator-side quantities are block | |
| 5 | +averages. Nothing here measures communication; see `communication.py` (accounting) and `protocol.py` (real runs). | |
| 6 | +""" | |
| 7 | +from __future__ import annotations | |
| 8 | +import numpy as np | |
| 9 | +from .jackknife import vandermonde_weights, DEFAULT_SPLITS | |
| 10 | +from .likelihoods import STATUS_OK | |
| 11 | + | |
| 12 | + | |
| 13 | +# ----------------------------------------------------------------------------- local fits | |
| 14 | +def local_split_fits(model, db, splits=(1,), **fit_kw): | |
| 15 | + """Fit the MLE on every sub-block of every block for each split multiplier s in `splits`. | |
| 16 | + | |
| 17 | + Returns dict s -> (theta (K, s, p), status (K, s), hessian (K, s, p, p)). | |
| 18 | + """ | |
| 19 | + K = next(iter(db.values())).shape[0] | |
| 20 | + m = model.n_obs(db) | |
| 21 | + out = {} | |
| 22 | + for s in splits: | |
| 23 | + if m % s: | |
| 24 | + raise ValueError(f"m={m} not divisible by s={s}") | |
| 25 | + sub = {k: v.reshape((K * s, m // s) + v.shape[2:]) for k, v in db.items()} | |
| 26 | + th, info = model.fit_batched(sub, **fit_kw) | |
| 27 | + out[s] = (th.reshape(K, s, -1), info["status"].reshape(K, s), info["hessian"].reshape(K, s, th.shape[1], th.shape[1])) | |
| 28 | + return out | |
| 29 | + | |
| 30 | + | |
| 31 | +def worker_corrected_estimates(fits, splits): | |
| 32 | + """J_{r,k} = sum_j w_j * mean_l theta_hat(sub-block l of level j) -> (K, p); also validity mask (K,).""" | |
| 33 | + w = vandermonde_weights(splits) | |
| 34 | + J = 0.0; ok = None | |
| 35 | + for wj, s in zip(w, splits): | |
| 36 | + th, st, _ = fits[s] | |
| 37 | + J = J + wj * th.mean(axis=1) | |
| 38 | + good = np.all(st == STATUS_OK, axis=1) | |
| 39 | + ok = good if ok is None else (ok & good) | |
| 40 | + return J, ok | |
| 41 | + | |
| 42 | + | |
| 43 | +def pilot_estimator(J, block_sizes=None, ok=None): | |
| 44 | + """Coordinator average of the worker estimates (weights m_k/N if sizes given). Invalid workers can be dropped.""" | |
| 45 | + J = np.asarray(J) | |
| 46 | + if block_sizes is None: | |
| 47 | + wts = np.ones(J.shape[0]) | |
| 48 | + else: | |
| 49 | + wts = np.asarray(block_sizes, float) | |
| 50 | + if ok is not None: | |
| 51 | + wts = wts * ok | |
| 52 | + wts = wts / wts.sum() | |
| 53 | + return np.einsum('k,kp->p', wts, J) | |
| 54 | + | |
| 55 | + | |
| 56 | +# ----------------------------------------------------------------------------- aggregated derivatives | |
| 57 | +def aggregated_score_hessian(model, db, theta, block_sizes=None): | |
| 58 | + """Workers evaluate (score, Hessian) of their block at the common theta; coordinator averages (weighted by m_k).""" | |
| 59 | + K = next(iter(db.values())).shape[0] | |
| 60 | + th = np.broadcast_to(np.asarray(theta, float), (K, len(theta))) | |
| 61 | + g, H = model.score_hessian(th, db) | |
| 62 | + if block_sizes is None: | |
| 63 | + return g.mean(axis=0), H.mean(axis=0), g, H | |
| 64 | + w = np.asarray(block_sizes, float); w = w / w.sum() | |
| 65 | + return np.einsum('k,kp->p', w, g), np.einsum('k,kpq->pq', w, H), g, H | |
| 66 | + | |
| 67 | + | |
| 68 | +def newton_step(theta, g, H, damping=None): | |
| 69 | + """theta - H^{-1} g via a linear solve (never an explicit inverse). Returns (new theta, diagnostics).""" | |
| 70 | + try: | |
| 71 | + step = np.linalg.solve(-H, g) | |
| 72 | + ok = np.all(np.isfinite(step)) | |
| 73 | + except np.linalg.LinAlgError: | |
| 74 | + step = np.zeros_like(g); ok = False | |
| 75 | + diag = {"solve_ok": ok, "cond": float(np.linalg.cond(-H)) if ok else np.inf, "step_norm": float(np.linalg.norm(step))} | |
| 76 | + if damping is not None and ok: | |
| 77 | + step = step * damping | |
| 78 | + return theta + step, diag | |
| 79 | + | |
| 80 | + | |
| 81 | +def global_newton(model, db, theta0, q=1, hessian_mode="exact", block_sizes=None, stale_hessian=None, | |
| 82 | + quantize_bits=None, single_machine=0, rng=None): | |
| 83 | + """q globally aggregated Newton corrections from theta0. | |
| 84 | + | |
| 85 | + hessian_mode: | |
| 86 | + 'exact' workers send score and Hessian at the broadcast iterate (p + p(p+1)/2 numbers each round); | |
| 87 | + 'stale' Hessian fixed = `stale_hessian` (e.g. average of local Hessians at local MLEs, sent once with the pilot); | |
| 88 | + workers send only scores each round (p numbers); | |
| 89 | + 'single' Hessian of worker `single_machine` evaluated at the iterate (no Hessian uplink; CSL-type curvature); | |
| 90 | + 'quantized' exact aggregated Hessian rounded to `quantize_bits` significant bits per entry (uniform relative | |
| 91 | + quantization), scores at full precision. | |
| 92 | + Returns list of iterates [theta0, theta1, ..., thetaq] and per-step diagnostics. | |
| 93 | + """ | |
| 94 | + th = np.asarray(theta0, float).copy(); path = [th.copy()]; diags = [] | |
| 95 | + for _ in range(q): | |
| 96 | + g, H, g_k, H_k = aggregated_score_hessian(model, db, th, block_sizes) | |
| 97 | + if hessian_mode == "stale": | |
| 98 | + H = np.asarray(stale_hessian) | |
| 99 | + elif hessian_mode == "single": | |
| 100 | + H = H_k[single_machine] | |
| 101 | + elif hessian_mode == "quantized": | |
| 102 | + H = quantize_relative(H, quantize_bits) | |
| 103 | + elif hessian_mode != "exact": | |
| 104 | + raise ValueError(hessian_mode) | |
| 105 | + th, d = newton_step(th, g, H) | |
| 106 | + d["grad_norm"] = float(np.linalg.norm(g)) | |
| 107 | + path.append(th.copy()); diags.append(d) | |
| 108 | + return path, diags | |
| 109 | + | |
| 110 | + | |
| 111 | +def quantize_relative(A, bits): | |
| 112 | + """Round each entry to `bits` significant binary digits (mantissa truncation), keeping the sign and exponent.""" | |
| 113 | + if bits is None: | |
| 114 | + return A | |
| 115 | + A = np.asarray(A, float) | |
| 116 | + mant, exp = np.frexp(A) | |
| 117 | + scale = 2.0 ** bits | |
| 118 | + return np.ldexp(np.round(mant * scale) / scale, exp) | |
| 119 | + | |
| 120 | + | |
| 121 | +# ----------------------------------------------------------------------------- baselines | |
| 122 | +def aee_lin_xi(theta_local, hess_local): | |
| 123 | + """Lin & Xi (2011) aggregated estimating equation: (sum_k A_k)^{-1} sum_k A_k theta_k with A_k = -H_k(theta_k).""" | |
| 124 | + A = -np.asarray(hess_local) # (K,p,p) | |
| 125 | + rhs = np.einsum('kpq,kq->p', A, theta_local) | |
| 126 | + return np.linalg.solve(A.sum(axis=0), rhs) | |
| 127 | + | |
| 128 | + | |
| 129 | +def csl(model, db, theta_init, rounds=1, machine=0, block_sizes=None, fit_kw=None): | |
| 130 | + """Jordan, Lee & Yang (2019) communication-efficient surrogate likelihood, maximized on `machine`. | |
| 131 | + | |
| 132 | + Round t: workers send gradients at theta_t (p numbers each); machine 0 maximizes | |
| 133 | + L_tilde(theta) = L_0(theta) - <grad L_0(theta_t) - grad L_N(theta_t), theta> | |
| 134 | + by Newton's method on its own block (local Hessian only). Returns iterates. | |
| 135 | + """ | |
| 136 | + fit_kw = fit_kw or {} | |
| 137 | + th = np.asarray(theta_init, float).copy(); path = [th.copy()] | |
| 138 | + d0 = {k: v[machine][None] for k, v in db.items()} | |
| 139 | + for _ in range(rounds): | |
| 140 | + gN, _, g_k, _ = aggregated_score_hessian(model, db, th, block_sizes) | |
| 141 | + shift = g_k[machine] - gN # p-vector | |
| 142 | + # maximize L_0(theta) - shift' theta with Newton on machine 0 | |
| 143 | + x = th.copy() | |
| 144 | + for it in range(100): | |
| 145 | + g, H = model.score_hessian(x[None], d0) | |
| 146 | + g = g[0] - shift; H = H[0] | |
| 147 | + step = np.linalg.solve(-H, g) | |
| 148 | + # backtracking on the surrogate objective | |
| 149 | + f = lambda z: model.loglik(z[None], d0)[0] - shift @ z | |
| 150 | + f0 = f(x); t = 1.0 | |
| 151 | + while f(x + t * step) < f0 - 1e-12 and t > 1e-8: | |
| 152 | + t *= 0.5 | |
| 153 | + x = x + t * step | |
| 154 | + if np.max(np.abs(g)) < 1e-10: | |
| 155 | + break | |
| 156 | + th = x; path.append(th.copy()) | |
| 157 | + return path | |
| 158 | + | |
| 159 | + | |
| 160 | +def analytic_bias_correction(model, db, theta_local): | |
| 161 | + """Cox–Snell / Cordeiro–McCullagh O(1/m) bias correction of each local MLE (GLM canonical link, or Exponential): | |
| 162 | + b_hat = (1/2) I^{-1} [n^{-1} sum_i nabla^2 psi_i (I^{-1})] per block; theta_corrected = theta - b_hat / m. | |
| 163 | + (Derivation: theta_hat - theta0 = -H^{-1} Psi - H^{-1}(V delta) - (1/2) H^{-1} nabla^2 Psi[delta,delta] + ...; with H = -I, | |
| 164 | + E[delta delta^T] = I^{-1}/n and the V-term vanishing conditionally on the design for canonical GLMs.) | |
| 165 | + Equivalent to the general second-order bias of M-estimators with the first (score-covariance) term dropped, | |
| 166 | + which is exact for canonical GLMs conditional on the design. | |
| 167 | + """ | |
| 168 | + K, m = next(iter(db.values())).shape[:2] | |
| 169 | + g, H = model.score_hessian(theta_local, db) | |
| 170 | + I = -H | |
| 171 | + Iinv = np.linalg.inv(I) # p small; explicit inverse acceptable for a diagnostic quantity | |
| 172 | + T = model.third_contraction(theta_local, db, Iinv) # (K,p) = n^{-1} sum nabla^2 psi_i [I^{-1}] | |
| 173 | + b = 0.5 * np.einsum("kpq,kq->kp", Iinv, T) # per-block O(1) bias coefficient b_1 = (1/2) I^{-1} T | |
| 174 | + return theta_local - b / m, b | |
| 175 | + | |
| 176 | + | |
| 177 | +# ----------------------------------------------------------------------------- one call = all estimators on one sample | |
| 178 | +def estimate_all(model, db, splits_max=3, q_max=2, block_sizes=None, fit_kw=None, extras=True): | |
| 179 | + """Compute, on the same partitioned sample, the full family of estimators used in the Monte Carlo. | |
| 180 | + | |
| 181 | + Returns dict name -> theta (p,), plus a dict of diagnostics (validity masks, Newton diagnostics, local fit stats). | |
| 182 | + Names: raw (r=1,q=0) ... 'r{r}_q{q}' for r=1..splits_max, q=0..q_max with splits (1,2,4,..); | |
| 183 | + 'central' centralized MLE; 'aee' Lin–Xi; 'csl1','csl2' CSL rounds from the raw average; | |
| 184 | + 'cs_q{q}' analytic Cox–Snell corrected average + q Newton; 'stale_r{r}_q{q}' stale-Hessian variants. | |
| 185 | + """ | |
| 186 | + fit_kw = fit_kw or {} | |
| 187 | + splits = DEFAULT_SPLITS[splits_max] | |
| 188 | + fits = local_split_fits(model, db, splits, **fit_kw) | |
| 189 | + K, m = next(iter(db.values())).shape[:2] | |
| 190 | + out, diag = {}, {} | |
| 191 | + diag["local_status_counts"] = {s: np.bincount(fits[s][1].ravel(), minlength=4).tolist() for s in splits} | |
| 192 | + # centralized MLE on the pooled data | |
| 193 | + full = {k: v.reshape((1, K * m) + v.shape[2:]) for k, v in db.items()} | |
| 194 | + th_c, info_c = model.fit_batched(full, tol=1e-12, maxit=100) | |
| 195 | + out["central"] = th_c[0]; diag["central_status"] = int(info_c["status"][0]) | |
| 196 | + theta_full_local, status_full, H_local = fits[1] | |
| 197 | + H_bar = H_local[:, 0].mean(axis=0) # average local Hessian at local MLEs (sent with the pilot) | |
| 198 | + for r in range(1, splits_max + 1): | |
| 199 | + sp = splits[:r] | |
| 200 | + J, ok = worker_corrected_estimates({s: fits[s] for s in sp}, sp) | |
| 201 | + diag[f"valid_r{r}"] = ok | |
| 202 | + th0 = pilot_estimator(J, block_sizes, ok=None) # keep all workers (failures tracked separately) | |
| 203 | + path, nd = global_newton(model, db, th0, q=q_max, block_sizes=block_sizes) | |
| 204 | + for q in range(q_max + 1): | |
| 205 | + out[f"r{r}_q{q}"] = path[q] | |
| 206 | + diag[f"newton_r{r}"] = nd | |
| 207 | + if extras: | |
| 208 | + path_s, _ = global_newton(model, db, th0, q=q_max, hessian_mode="stale", stale_hessian=H_bar, block_sizes=block_sizes) | |
| 209 | + for q in range(1, q_max + 1): | |
| 210 | + out[f"stale_r{r}_q{q}"] = path_s[q] | |
| 211 | + if extras: | |
| 212 | + out["aee"] = aee_lin_xi(theta_full_local[:, 0], H_local[:, 0]) | |
| 213 | + raw = out["r1_q0"] | |
| 214 | + cpath = csl(model, db, raw, rounds=2, block_sizes=block_sizes) | |
| 215 | + out["csl1"], out["csl2"] = cpath[1], cpath[2] | |
| 216 | + try: | |
| 217 | + th_cs, b = analytic_bias_correction(model, db, theta_full_local[:, 0]) | |
| 218 | + th0 = pilot_estimator(th_cs, block_sizes) | |
| 219 | + pth, _ = global_newton(model, db, th0, q=q_max, block_sizes=block_sizes) | |
| 220 | + for q in range(q_max + 1): | |
| 221 | + out[f"cs_q{q}"] = pth[q] | |
| 222 | + except NotImplementedError: | |
| 223 | + pass | |
| 224 | + return out, diag | |
added
src/pamle/communication.py
+52 −0
@@ -0,0 +1,52 @@ | ||
| 1 | +"""Accounting of rounds, messages and transmitted numbers/bytes for the protocols (star topology).""" | |
| 2 | +from __future__ import annotations | |
| 3 | + | |
| 4 | + | |
| 5 | +def protocol_accounting(K: int, p: int, r: int = 1, q: int = 1, hessian_mode: str = "exact", bytes_per_scalar: int = 8, | |
| 6 | + send_local_hessian: bool = False, header_bytes: int = 0): | |
| 7 | + """Numbers and bytes moved for: local pilot (phase 1) + q derivative phases. | |
| 8 | + | |
| 9 | + Phase 1 uplink per worker: p (corrected local estimate) [+ p(p+1)/2 if the local Hessian is shipped: stale/AEE]. | |
| 10 | + Each derivative phase: downlink broadcast of p to every worker, uplink of p (score) [+ p(p+1)/2 Hessian if exact]. | |
| 11 | + hessian_mode: 'exact' | 'stale' | 'single' (no Hessian uplink; machine 0 solves) | 'none' (q = 0). | |
| 12 | + Returns dict with rounds, messages, uplink/downlink numbers and bytes, local fits per worker. | |
| 13 | + """ | |
| 14 | + h = p * (p + 1) // 2 | |
| 15 | + up1 = p + (h if (send_local_hessian or hessian_mode == "stale") else 0) | |
| 16 | + if q == 0: | |
| 17 | + up_d, down_d = 0, 0 | |
| 18 | + elif hessian_mode == "exact" or hessian_mode == "quantized": | |
| 19 | + up_d, down_d = q * (p + h), q * p | |
| 20 | + elif hessian_mode in ("stale", "single"): | |
| 21 | + up_d, down_d = q * p, q * p | |
| 22 | + else: | |
| 23 | + raise ValueError(hessian_mode) | |
| 24 | + uplink_numbers = K * (up1 + up_d) | |
| 25 | + downlink_numbers = K * down_d | |
| 26 | + rounds = 1 + q # synchronized exchange phases (uplink of estimates, then q derivative phases) | |
| 27 | + messages = K * (1 + 2 * q) # K uploads + per phase: K broadcasts + K uploads | |
| 28 | + local_fits = 2 ** r - 1 # splits (1,2,4,...): 1 + 2 + 4 + ... | |
| 29 | + local_data_passes = r # each level touches every local observation once | |
| 30 | + return {"rounds": rounds, "messages": messages, "uplink_numbers": uplink_numbers, "downlink_numbers": downlink_numbers, | |
| 31 | + "uplink_bytes": uplink_numbers * bytes_per_scalar + messages * header_bytes, | |
| 32 | + "downlink_bytes": downlink_numbers * bytes_per_scalar, | |
| 33 | + "local_fits_per_worker": local_fits, "local_data_passes": local_data_passes} | |
| 34 | + | |
| 35 | + | |
| 36 | +def sufficient_exponent(r: int, q: int) -> float: | |
| 37 | + """alpha(r,q) = 1 - 1/(2^(q+1) r): K = o(N^alpha) is the sufficient regime of Theorem 2 (exact Newton).""" | |
| 38 | + return 1.0 - 1.0 / (2 ** (q + 1) * r) | |
| 39 | + | |
| 40 | + | |
| 41 | +def sufficient_exponent_stale(r: int, q: int) -> float: | |
| 42 | + """Sufficient exponent with a fixed (stale) averaged local Hessian: effective depth r + q (Theorem 4, case (b)).""" | |
| 43 | + if q == 0: | |
| 44 | + return sufficient_exponent(r, 0) | |
| 45 | + return 1.0 - 1.0 / (2.0 * (r + q)) | |
| 46 | + | |
| 47 | + | |
| 48 | +def sufficient_exponent_single(r: int, q: int) -> float: | |
| 49 | + """Sufficient exponent with a single machine's Hessian (CSL-type curvature): effective depth r + q/2.""" | |
| 50 | + if q == 0: | |
| 51 | + return sufficient_exponent(r, 0) | |
| 52 | + return 1.0 - 1.0 / (2.0 * r + q) | |
added
src/pamle/diagnostics.py
+32 −0
@@ -0,0 +1,32 @@ | ||
| 1 | +"""Monte Carlo summaries, confidence intervals and coverage.""" | |
| 2 | +from __future__ import annotations | |
| 3 | +import numpy as np | |
| 4 | +from scipy import stats | |
| 5 | + | |
| 6 | + | |
| 7 | +def summarize(errors: np.ndarray): | |
| 8 | + """errors: (B,) or (B,p) array -> dict with bias, sd, rmse, MC standard errors.""" | |
| 9 | + e = np.asarray(errors, float) | |
| 10 | + B = e.shape[0] | |
| 11 | + mean = e.mean(axis=0); sd = e.std(axis=0, ddof=1) | |
| 12 | + rmse = np.sqrt(np.mean(e ** 2, axis=0)) | |
| 13 | + return {"n": B, "bias": mean, "sd": sd, "rmse": rmse, "bias_se": sd / np.sqrt(B), | |
| 14 | + "rmse_se": np.sqrt(np.var(e ** 2, axis=0, ddof=1) / B) / (2 * np.maximum(rmse, 1e-300)), | |
| 15 | + "p95_abs": np.quantile(np.abs(e), 0.95, axis=0)} | |
| 16 | + | |
| 17 | + | |
| 18 | +def wald_ci(theta, cov, N, level=0.95): | |
| 19 | + z = stats.norm.ppf(0.5 + level / 2) | |
| 20 | + se = np.sqrt(np.diag(cov) / N) | |
| 21 | + return theta - z * se, theta + z * se | |
| 22 | + | |
| 23 | + | |
| 24 | +def coverage(theta_hat, cov_hat, N, theta0, levels=(0.90, 0.95, 0.99)): | |
| 25 | + """theta_hat (B,p), cov_hat (B,p,p) [asymptotic covariance of sqrt(N)(theta-theta0)] -> coverage per level and param.""" | |
| 26 | + out = {} | |
| 27 | + for lv in levels: | |
| 28 | + z = stats.norm.ppf(0.5 + lv / 2) | |
| 29 | + se = np.sqrt(np.einsum('bpp->bp', cov_hat) / N) | |
| 30 | + hit = np.abs(theta_hat - theta0) <= z * se | |
| 31 | + out[lv] = {"coverage": hit.mean(axis=0), "mean_length": (2 * z * se).mean(axis=0)} | |
| 32 | + return out | |
added
src/pamle/jackknife.py
+70 −0
@@ -0,0 +1,70 @@ | ||
| 1 | +"""Split-granularity extrapolation weights (generalized jackknife / Richardson extrapolation). | |
| 2 | + | |
| 3 | +For split multipliers s_1 < ... < s_r (s_j sub-blocks of size m/s_j per worker; s_1 = 1 is the full block), | |
| 4 | +the weights solve the Vandermonde system | |
| 5 | + sum_j w_j = 1, sum_j w_j s_j^a = 0, a = 1, ..., r-1, | |
| 6 | +so that a bias expansion b_1/n + ... + b_{r-1}/n^{r-1} evaluated at n = m/s_j cancels through order m^{-(r-1)}. | |
| 7 | +""" | |
| 8 | +from __future__ import annotations | |
| 9 | +from fractions import Fraction | |
| 10 | +import numpy as np | |
| 11 | + | |
| 12 | +DEFAULT_SPLITS = {1: (1,), 2: (1, 2), 3: (1, 2, 4), 4: (1, 2, 4, 8), 5: (1, 2, 4, 8, 16)} | |
| 13 | + | |
| 14 | + | |
| 15 | +def vandermonde_weights_exact(splits): | |
| 16 | + """Exact rational weights (list of Fraction) for the split multipliers `splits`.""" | |
| 17 | + s = [Fraction(int(x)) for x in splits] | |
| 18 | + r = len(s) | |
| 19 | + # Solve by Gaussian elimination on exact rationals. | |
| 20 | + A = [[sj ** a for sj in s] for a in range(r)] | |
| 21 | + b = [Fraction(1)] + [Fraction(0)] * (r - 1) | |
| 22 | + M = [row[:] + [bi] for row, bi in zip(A, b)] | |
| 23 | + for col in range(r): | |
| 24 | + piv = next(i for i in range(col, r) if M[i][col] != 0) | |
| 25 | + M[col], M[piv] = M[piv], M[col] | |
| 26 | + pv = M[col][col] | |
| 27 | + M[col] = [x / pv for x in M[col]] | |
| 28 | + for i in range(r): | |
| 29 | + if i != col and M[i][col] != 0: | |
| 30 | + f = M[i][col] | |
| 31 | + M[i] = [x - f * y for x, y in zip(M[i], M[col])] | |
| 32 | + return [M[i][r] for i in range(r)] | |
| 33 | + | |
| 34 | + | |
| 35 | +def vandermonde_weights(splits): | |
| 36 | + """Float64 weights.""" | |
| 37 | + return np.array([float(w) for w in vandermonde_weights_exact(splits)]) | |
| 38 | + | |
| 39 | + | |
| 40 | +def leading_bias_coefficient(splits): | |
| 41 | + """beta_r = sum_j w_j s_j^r: multiplies b_r m^{-r} in the residual bias of J_r (exact rational).""" | |
| 42 | + w = vandermonde_weights_exact(splits); r = len(splits) | |
| 43 | + return sum(wj * Fraction(int(sj)) ** r for wj, sj in zip(w, splits)) | |
| 44 | + | |
| 45 | + | |
| 46 | +def amplification_constant(splits): | |
| 47 | + """A_r = sum_j |w_j| s_j: multiplies sqrt(C)/m in the L2 bound of the corrected remainder (Lemma 3).""" | |
| 48 | + w = vandermonde_weights_exact(splits) | |
| 49 | + return sum(abs(wj) * Fraction(int(sj)) for wj, sj in zip(w, splits)) | |
| 50 | + | |
| 51 | + | |
| 52 | +def weight_l1(splits): | |
| 53 | + return sum(abs(w) for w in vandermonde_weights_exact(splits)) | |
| 54 | + | |
| 55 | + | |
| 56 | +def vandermonde_condition_number(splits): | |
| 57 | + s = np.asarray(splits, dtype=float); r = len(s) | |
| 58 | + A = np.vstack([s ** a for a in range(r)]) | |
| 59 | + return float(np.linalg.cond(A)) | |
| 60 | + | |
| 61 | + | |
| 62 | +def combine_levels(level_means, weights): | |
| 63 | + """J = sum_j w_j T_{s_j}. level_means: array (r, ..., p) of per-level averages; weights: (r,).""" | |
| 64 | + w = np.asarray(weights, dtype=float) | |
| 65 | + return np.tensordot(w, np.asarray(level_means), axes=(0, 0)) | |
| 66 | + | |
| 67 | + | |
| 68 | +def level_averages(estimates_by_level): | |
| 69 | + """estimates_by_level[j]: array (K, s_j, p) of sub-block estimates at level j -> (r, K, p) level means.""" | |
| 70 | + return np.stack([np.mean(E, axis=1) for E in estimates_by_level], axis=0) | |
added
src/pamle/likelihoods.py
+246 −0
@@ -0,0 +1,246 @@ | ||
| 1 | +"""Likelihood models with batched (per-block) Newton fitting. | |
| 2 | + | |
| 3 | +Data layout | |
| 4 | +----------- | |
| 5 | +A dataset is a dict of arrays whose leading axis indexes observations: GLMs use {'X': (n, p), 'y': (n,)}; | |
| 6 | +the exponential model uses {'x': (n,)}. A *batched* dataset has one extra leading axis for blocks: | |
| 7 | +{'X': (B, n, p), 'y': (B, n)}. All model methods accept batched data and batched parameters (B, p). | |
| 8 | + | |
| 9 | +Status codes returned by `fit_batched`: 0 converged, 1 max iterations, 2 singular/indefinite Hessian, | |
| 10 | +3 non-finite or diverged (e.g. logistic separation). | |
| 11 | +""" | |
| 12 | +from __future__ import annotations | |
| 13 | +import numpy as np | |
| 14 | + | |
| 15 | +STATUS_OK, STATUS_MAXIT, STATUS_SINGULAR, STATUS_DIVERGED = 0, 1, 2, 3 | |
| 16 | + | |
| 17 | + | |
| 18 | +def gather_blocks(data: dict, blocks: np.ndarray) -> dict: | |
| 19 | + """Index a dataset by an integer array `blocks` of shape (K, m) -> batched dataset (K, m, ...).""" | |
| 20 | + return {k: v[blocks] for k, v in data.items()} | |
| 21 | + | |
| 22 | + | |
| 23 | +def subset(data: dict, idx) -> dict: | |
| 24 | + return {k: v[idx] for k, v in data.items()} | |
| 25 | + | |
| 26 | + | |
| 27 | +def _batched_solve(A, b): | |
| 28 | + """Solve A x = b for batched A (B,p,p), b (B,p); returns x and a boolean mask of failures.""" | |
| 29 | + B, p = b.shape | |
| 30 | + x = np.zeros_like(b); fail = np.zeros(B, dtype=bool) | |
| 31 | + try: | |
| 32 | + x = np.linalg.solve(A, b[..., None])[..., 0] | |
| 33 | + bad = ~np.all(np.isfinite(x), axis=1) | |
| 34 | + except np.linalg.LinAlgError: | |
| 35 | + bad = np.ones(B, dtype=bool) | |
| 36 | + if np.any(bad): | |
| 37 | + for i in np.flatnonzero(bad): | |
| 38 | + try: | |
| 39 | + x[i] = np.linalg.solve(A[i], b[i]); fail[i] = not np.all(np.isfinite(x[i])) | |
| 40 | + except np.linalg.LinAlgError: | |
| 41 | + fail[i] = True; x[i] = 0.0 | |
| 42 | + return x, fail | |
| 43 | + | |
| 44 | + | |
| 45 | +class Model: | |
| 46 | + """Abstract model. Subclasses implement the vectorised primitives below.""" | |
| 47 | + name = "abstract" | |
| 48 | + p = None | |
| 49 | + | |
| 50 | + # ---- primitives (batched: theta (B,p), data batched) ---- | |
| 51 | + def loglik(self, theta, data): # -> (B,) | |
| 52 | + raise NotImplementedError | |
| 53 | + def score_hessian(self, theta, data): # -> (B,p), (B,p,p) averages per block | |
| 54 | + raise NotImplementedError | |
| 55 | + def third_contraction(self, theta, data, Sigma): | |
| 56 | + """n^{-1} sum_i nabla^2 psi_i(theta)[Sigma] -> (B,p); psi_i = per-observation score.""" | |
| 57 | + raise NotImplementedError | |
| 58 | + def outer_scores(self, theta, data): | |
| 59 | + """n^{-1} sum_i psi_i psi_i' -> (B,p,p) (for sandwich covariances).""" | |
| 60 | + raise NotImplementedError | |
| 61 | + def simulate(self, rng, n, theta0, **kw) -> dict: | |
| 62 | + raise NotImplementedError | |
| 63 | + def init(self, data): # starting value (B,p) | |
| 64 | + raise NotImplementedError | |
| 65 | + def n_obs(self, data): | |
| 66 | + raise NotImplementedError | |
| 67 | + def degenerate(self, theta, data): | |
| 68 | + """Boolean mask (B,): a stationary point that should NOT be accepted as an MLE (e.g. logistic separation: | |
| 69 | + fitted probabilities numerically 0/1 so the score vanishes while the Hessian collapses).""" | |
| 70 | + return np.zeros(theta.shape[0], dtype=bool) | |
| 71 | + | |
| 72 | + # ---- generic Newton with step halving ---- | |
| 73 | + def fit_batched(self, data, init=None, tol=1e-9, maxit=60, max_norm=1e4): | |
| 74 | + theta = np.array(self.init(data) if init is None else init, dtype=float) | |
| 75 | + if theta.ndim == 1: | |
| 76 | + theta = theta[None, :] | |
| 77 | + B = theta.shape[0] | |
| 78 | + status = np.full(B, STATUS_MAXIT, dtype=int) | |
| 79 | + iters = np.zeros(B, dtype=int) | |
| 80 | + active = np.ones(B, dtype=bool) | |
| 81 | + ll = self.loglik(theta, data) | |
| 82 | + for it in range(maxit): | |
| 83 | + g, H = self.score_hessian(theta, data) | |
| 84 | + gmax = np.max(np.abs(g), axis=1) | |
| 85 | + conv = active & (gmax < tol) | |
| 86 | + if np.any(conv): | |
| 87 | + deg = conv & self.degenerate(theta, data) | |
| 88 | + status[deg] = STATUS_DIVERGED; conv &= ~deg | |
| 89 | + status[conv] = STATUS_OK; active &= ~(conv | (status == STATUS_DIVERGED)) | |
| 90 | + if not np.any(active): | |
| 91 | + break | |
| 92 | + step, fail = _batched_solve(-H, g) | |
| 93 | + if np.any(fail & active): | |
| 94 | + status[fail & active] = STATUS_SINGULAR; active &= ~fail | |
| 95 | + # ascent check: if the direction is not an ascent direction, fall back to gradient | |
| 96 | + asc = np.einsum('bp,bp->b', step, g) | |
| 97 | + bad = active & (asc <= 0) | |
| 98 | + step[bad] = g[bad] | |
| 99 | + t = np.ones(B) | |
| 100 | + new = theta + t[:, None] * step | |
| 101 | + ll_new = self.loglik(new, data) | |
| 102 | + for _ in range(30): | |
| 103 | + worse = active & (~np.isfinite(ll_new) | (ll_new < ll - 1e-12)) | |
| 104 | + if not np.any(worse): | |
| 105 | + break | |
| 106 | + t[worse] *= 0.5 | |
| 107 | + new = theta + t[:, None] * step | |
| 108 | + ll_new = self.loglik(new, data) | |
| 109 | + theta = np.where(active[:, None], new, theta) | |
| 110 | + ll = np.where(active, ll_new, ll) | |
| 111 | + iters[active] += 1 | |
| 112 | + div = active & (~np.all(np.isfinite(theta), axis=1) | (np.linalg.norm(theta, axis=1) > max_norm)) | |
| 113 | + status[div] = STATUS_DIVERGED; active &= ~div | |
| 114 | + g, H = self.score_hessian(theta, data) | |
| 115 | + info = {"status": status, "iters": iters, "grad_max": np.max(np.abs(g), axis=1), | |
| 116 | + "hessian": H, "cond": np.linalg.cond(-H)} | |
| 117 | + return theta, info | |
| 118 | + | |
| 119 | + def fit(self, data, init=None, **kw): | |
| 120 | + db = {k: v[None] for k, v in data.items()} | |
| 121 | + th, info = self.fit_batched(db, init=None if init is None else np.asarray(init)[None], **kw) | |
| 122 | + return th[0], {k: (v[0] if isinstance(v, np.ndarray) else v) for k, v in info.items()} | |
| 123 | + | |
| 124 | + | |
| 125 | +class Exponential(Model): | |
| 126 | + """X ~ Exp(rate theta); average log-likelihood log(theta) - theta * xbar; MLE 1/xbar.""" | |
| 127 | + name = "exponential"; p = 1 | |
| 128 | + | |
| 129 | + def n_obs(self, data): | |
| 130 | + return data['x'].shape[-1] | |
| 131 | + def simulate(self, rng, n, theta0, **kw): | |
| 132 | + return {'x': rng.exponential(scale=1.0 / float(np.ravel(theta0)[0]), size=n)} | |
| 133 | + def init(self, data): | |
| 134 | + return 1.0 / np.mean(data['x'], axis=-1)[:, None] | |
| 135 | + def loglik(self, theta, data): | |
| 136 | + th = theta[:, 0] | |
| 137 | + return np.log(th) - th * np.mean(data['x'], axis=-1) | |
| 138 | + def score_hessian(self, theta, data): | |
| 139 | + th = theta[:, 0] | |
| 140 | + g = (1.0 / th - np.mean(data['x'], axis=-1))[:, None] | |
| 141 | + H = (-1.0 / th ** 2)[:, None, None] | |
| 142 | + return g, H | |
| 143 | + def third_contraction(self, theta, data, Sigma): | |
| 144 | + th = theta[:, 0] | |
| 145 | + return (2.0 / th ** 3 * Sigma[:, 0, 0])[:, None] | |
| 146 | + def outer_scores(self, theta, data): | |
| 147 | + th = theta[:, 0][:, None] | |
| 148 | + psi = 1.0 / th - data['x'] | |
| 149 | + return np.mean(psi ** 2, axis=-1)[:, None, None] | |
| 150 | + def mle(self, data): | |
| 151 | + """Closed form, batched.""" | |
| 152 | + return (1.0 / np.mean(data['x'], axis=-1))[:, None] | |
| 153 | + | |
| 154 | + | |
| 155 | +class _GLM(Model): | |
| 156 | + """Canonical-link GLM with per-observation score x_i (y_i - mu_i), observed Hessian -sum x_i x_i' mu'(eta_i)/n.""" | |
| 157 | + def n_obs(self, data): | |
| 158 | + return data['y'].shape[-1] | |
| 159 | + def _eta(self, theta, data): | |
| 160 | + return np.einsum('bnp,bp->bn', data['X'], theta) | |
| 161 | + def score_hessian(self, theta, data): | |
| 162 | + X = data['X']; y = data['y']; n = y.shape[-1] | |
| 163 | + eta = self._eta(theta, data) | |
| 164 | + mu, w = self._mu_w(eta) | |
| 165 | + g = np.einsum('bnp,bn->bp', X, y - mu) / n | |
| 166 | + H = -np.einsum('bnp,bn,bnq->bpq', X, w, X) / n | |
| 167 | + return g, H | |
| 168 | + def third_contraction(self, theta, data, Sigma): | |
| 169 | + X = data['X']; n = X.shape[1] | |
| 170 | + eta = self._eta(theta, data) | |
| 171 | + mu2 = self._mu2(eta) # mu''(eta) | |
| 172 | + quad = np.einsum('bnp,bpq,bnq->bn', X, Sigma, X) # x_i' Sigma x_i | |
| 173 | + return -np.einsum('bnp,bn,bn->bp', X, mu2, quad) / n | |
| 174 | + def outer_scores(self, theta, data): | |
| 175 | + X = data['X']; y = data['y']; n = y.shape[-1] | |
| 176 | + mu, _ = self._mu_w(self._eta(theta, data)) | |
| 177 | + r = y - mu | |
| 178 | + return np.einsum('bnp,bn,bnq->bpq', X, r * r, X) / n | |
| 179 | + def fisher(self, theta, data): | |
| 180 | + """Expected information = -observed Hessian for canonical links (B,p,p).""" | |
| 181 | + _, H = self.score_hessian(theta, data); return -H | |
| 182 | + | |
| 183 | + | |
| 184 | +class Logistic(_GLM): | |
| 185 | + name = "logistic" | |
| 186 | + def __init__(self, p): | |
| 187 | + self.p = p | |
| 188 | + def init(self, data): | |
| 189 | + B = data['y'].shape[0] | |
| 190 | + return np.zeros((B, self.p)) | |
| 191 | + def _mu_w(self, eta): | |
| 192 | + mu = 1.0 / (1.0 + np.exp(-eta)); return mu, mu * (1.0 - mu) | |
| 193 | + def _mu2(self, eta): | |
| 194 | + mu = 1.0 / (1.0 + np.exp(-eta)); return mu * (1.0 - mu) * (1.0 - 2.0 * mu) | |
| 195 | + def degenerate(self, theta, data, eta_max=30.0): | |
| 196 | + """Separation safeguard: some fitted |eta| > eta_max (probability within 1e-13 of 0 or 1) or a numerically | |
| 197 | + singular information matrix. Legitimate designs in this paper never reach |eta| = 30.""" | |
| 198 | + eta = self._eta(theta, data) | |
| 199 | + _, H = self.score_hessian(theta, data) | |
| 200 | + return (np.max(np.abs(eta), axis=-1) > eta_max) | (np.linalg.cond(-H) > 1e12) | |
| 201 | + def loglik(self, theta, data): | |
| 202 | + eta = self._eta(theta, data); y = data['y'] | |
| 203 | + return np.mean(y * eta - np.logaddexp(0.0, eta), axis=-1) | |
| 204 | + def simulate(self, rng, n, theta0, cov=None, intercept=True, **kw): | |
| 205 | + theta0 = np.asarray(theta0, float); p = theta0.size | |
| 206 | + q = p - 1 if intercept else p | |
| 207 | + Z = rng.standard_normal((n, q)) if cov is None else rng.multivariate_normal(np.zeros(q), cov, size=n) | |
| 208 | + X = np.hstack([np.ones((n, 1)), Z]) if intercept else Z | |
| 209 | + eta = X @ theta0 | |
| 210 | + y = (rng.random(n) < 1.0 / (1.0 + np.exp(-eta))).astype(float) | |
| 211 | + return {'X': X, 'y': y} | |
| 212 | + | |
| 213 | + | |
| 214 | +class Poisson(_GLM): | |
| 215 | + name = "poisson" | |
| 216 | + def __init__(self, p): | |
| 217 | + self.p = p | |
| 218 | + def init(self, data): | |
| 219 | + B = data['y'].shape[0] | |
| 220 | + th = np.zeros((B, self.p)); th[:, 0] = np.log(np.maximum(np.mean(data['y'], axis=-1), 1e-3)) | |
| 221 | + return th | |
| 222 | + def _mu_w(self, eta): | |
| 223 | + mu = np.exp(np.minimum(eta, 50.0)); return mu, mu | |
| 224 | + def _mu2(self, eta): | |
| 225 | + return np.exp(np.minimum(eta, 50.0)) | |
| 226 | + def loglik(self, theta, data): | |
| 227 | + eta = self._eta(theta, data); y = data['y'] | |
| 228 | + return np.mean(y * eta - np.exp(np.minimum(eta, 50.0)), axis=-1) | |
| 229 | + def simulate(self, rng, n, theta0, cov=None, intercept=True, **kw): | |
| 230 | + theta0 = np.asarray(theta0, float); p = theta0.size | |
| 231 | + q = p - 1 if intercept else p | |
| 232 | + Z = rng.standard_normal((n, q)) if cov is None else rng.multivariate_normal(np.zeros(q), cov, size=n) | |
| 233 | + X = np.hstack([np.ones((n, 1)), Z]) if intercept else Z | |
| 234 | + y = rng.poisson(np.exp(X @ theta0)).astype(float) | |
| 235 | + return {'X': X, 'y': y} | |
| 236 | + | |
| 237 | + | |
| 238 | +def make_model(name: str, p: int = 1) -> Model: | |
| 239 | + name = name.lower() | |
| 240 | + if name == "exponential": | |
| 241 | + return Exponential() | |
| 242 | + if name == "logistic": | |
| 243 | + return Logistic(p) | |
| 244 | + if name == "poisson": | |
| 245 | + return Poisson(p) | |
| 246 | + raise ValueError(name) | |
added
src/pamle/partitioning.py
+48 −0
@@ -0,0 +1,48 @@ | ||
| 1 | +"""Index construction for blocks and nested sub-blocks.""" | |
| 2 | +from __future__ import annotations | |
| 3 | +import numpy as np | |
| 4 | + | |
| 5 | + | |
| 6 | +def equal_blocks(N: int, K: int): | |
| 7 | + """Contiguous equal blocks; requires K | N. Returns array (K, m) of indices.""" | |
| 8 | + if N % K: | |
| 9 | + raise ValueError(f"K={K} does not divide N={N}") | |
| 10 | + m = N // K | |
| 11 | + return np.arange(N).reshape(K, m) | |
| 12 | + | |
| 13 | + | |
| 14 | +def random_blocks(N: int, K: int, rng: np.random.Generator): | |
| 15 | + """Random equal blocks (uniform permutation), (K, m).""" | |
| 16 | + return rng.permutation(N).reshape(K, N // K) | |
| 17 | + | |
| 18 | + | |
| 19 | +def unequal_blocks(sizes): | |
| 20 | + """Contiguous blocks with the given sizes -> list of index arrays.""" | |
| 21 | + sizes = [int(s) for s in sizes] | |
| 22 | + out, start = [], 0 | |
| 23 | + for s in sizes: | |
| 24 | + out.append(np.arange(start, start + s)); start += s | |
| 25 | + return out | |
| 26 | + | |
| 27 | + | |
| 28 | +def nested_sub_blocks(block: np.ndarray, s: int): | |
| 29 | + """Split one block (length m) into s equal contiguous sub-blocks: (s, m/s).""" | |
| 30 | + m = block.shape[-1] | |
| 31 | + if m % s: | |
| 32 | + raise ValueError(f"split multiplier {s} does not divide block size {m}") | |
| 33 | + return block.reshape(block.shape[:-1] + (s, m // s)) | |
| 34 | + | |
| 35 | + | |
| 36 | +def check_splits(m: int, splits): | |
| 37 | + for s in splits: | |
| 38 | + if m % s: | |
| 39 | + raise ValueError(f"block size m={m} not divisible by split multiplier s={s}") | |
| 40 | + return True | |
| 41 | + | |
| 42 | + | |
| 43 | +def skewed_sizes(N: int, K: int, ratio: float): | |
| 44 | + """Unequal block sizes: sizes grow geometrically so that max/min ~= ratio; sum = N (integer-rounded).""" | |
| 45 | + g = np.geomspace(1.0, ratio, K) | |
| 46 | + sizes = np.floor(g / g.sum() * N).astype(int) | |
| 47 | + sizes[-1] += N - sizes.sum() | |
| 48 | + return sizes | |
added
tests/test_communication_accounting.py
+22 −0
@@ -0,0 +1,22 @@ | ||
| 1 | +from pamle.communication import protocol_accounting, sufficient_exponent, sufficient_exponent_stale, sufficient_exponent_single | |
| 2 | + | |
| 3 | + | |
| 4 | +def test_exponents_table(): | |
| 5 | + assert sufficient_exponent(1, 0) == 0.5 and sufficient_exponent(1, 1) == 0.75 and sufficient_exponent(1, 2) == 0.875 | |
| 6 | + assert sufficient_exponent(2, 0) == 0.75 and sufficient_exponent(2, 1) == 0.875 and sufficient_exponent(3, 1) == 11 / 12 | |
| 7 | + assert abs(sufficient_exponent_stale(2, 1) - 5 / 6) < 1e-15 and abs(sufficient_exponent_single(1, 1) - 2 / 3) < 1e-15 | |
| 8 | + assert sufficient_exponent_stale(1, 1) == sufficient_exponent(1, 1) # r=1: both 3/4 | |
| 9 | + | |
| 10 | + | |
| 11 | +def test_accounting_exact_vs_stale(): | |
| 12 | + K, p = 100, 10 | |
| 13 | + e = protocol_accounting(K, p, r=2, q=1, hessian_mode="exact") | |
| 14 | + s = protocol_accounting(K, p, r=2, q=1, hessian_mode="stale") | |
| 15 | + assert e["uplink_numbers"] == K * (p + p + p * (p + 1) // 2) | |
| 16 | + assert s["uplink_numbers"] == K * (p + p * (p + 1) // 2 + p) | |
| 17 | + e2 = protocol_accounting(K, p, r=1, q=2, hessian_mode="exact") | |
| 18 | + s2 = protocol_accounting(K, p, r=1, q=2, hessian_mode="stale") | |
| 19 | + assert s2["uplink_numbers"] < e2["uplink_numbers"] # stale Hessian saves bytes only when q >= 2 | |
| 20 | + assert e["rounds"] == 2 and e2["rounds"] == 3 | |
| 21 | + assert protocol_accounting(K, p, r=3, q=0)["local_fits_per_worker"] == 7 | |
| 22 | + assert protocol_accounting(K, p, r=1, q=0)["uplink_numbers"] == K * p | |
added
tests/test_estimators.py
+78 −0
@@ -0,0 +1,78 @@ | ||
| 1 | +import numpy as np | |
| 2 | +from pamle.likelihoods import Exponential, Logistic, gather_blocks | |
| 3 | +from pamle.algorithms import (local_split_fits, worker_corrected_estimates, pilot_estimator, global_newton, | |
| 4 | + aee_lin_xi, csl, analytic_bias_correction, estimate_all, quantize_relative) | |
| 5 | + | |
| 6 | + | |
| 7 | +def _exp_design(seed=0, N=2 ** 16, K=2 ** 8): | |
| 8 | + rng = np.random.default_rng(seed); ex = Exponential() | |
| 9 | + d = ex.simulate(rng, N, 1.0) | |
| 10 | + return ex, d, gather_blocks(d, np.arange(N).reshape(K, N // K)), N, K | |
| 11 | + | |
| 12 | + | |
| 13 | +def test_newton_identity_on_data(): | |
| 14 | + ex, d, db, N, K = _exp_design() | |
| 15 | + th_c = 1 / d['x'].mean() | |
| 16 | + for x in (0.7, 1.0, 1.4): | |
| 17 | + path, diag = global_newton(ex, db, np.array([x]), q=1) | |
| 18 | + assert abs(path[1][0] - (2 * x - x * x * d['x'].mean())) < 1e-12 | |
| 19 | + assert abs((path[1][0] - th_c) + (x - th_c) ** 2 / th_c) < 1e-12 | |
| 20 | + | |
| 21 | + | |
| 22 | +def test_pilot_bias_orders_exponential(): | |
| 23 | + """Average over many samples: raw bias ~ theta0/m, jackknife bias ~ -2 theta0/m^2 (exact formulas).""" | |
| 24 | + rng = np.random.default_rng(1); ex = Exponential(); m = 32; K = 2000; reps = 200 | |
| 25 | + raw, j2, j3 = [], [], [] | |
| 26 | + for _ in range(reps): | |
| 27 | + d = ex.simulate(rng, m * K, 1.0); db = gather_blocks(d, np.arange(m * K).reshape(K, m)) | |
| 28 | + fits = local_split_fits(ex, db, (1, 2, 4)) | |
| 29 | + raw.append(pilot_estimator(fits[1][0][:, 0])) | |
| 30 | + j2.append(pilot_estimator(worker_corrected_estimates(fits, (1, 2))[0])) | |
| 31 | + j3.append(pilot_estimator(worker_corrected_estimates(fits, (1, 2, 4))[0])) | |
| 32 | + raw, j2, j3 = np.array(raw).ravel(), np.array(j2).ravel(), np.array(j3).ravel() | |
| 33 | + se = 1 / np.sqrt(m * K * reps) | |
| 34 | + assert abs(raw.mean() - m / (m - 1)) < 4 * se | |
| 35 | + assert abs(j2.mean() - (1 - 2 / ((m - 1) * (m - 2)))) < 4 * se | |
| 36 | + assert abs(j3.mean() - (1 + 8 / ((m - 1) * (m - 2) * (m - 4)))) < 4 * se | |
| 37 | + | |
| 38 | + | |
| 39 | +def test_estimate_all_agrees_with_central_when_predicted(): | |
| 40 | + rng = np.random.default_rng(2); m = Logistic(3); N = 2 ** 15; K = 64 | |
| 41 | + d = m.simulate(rng, N, [0.5, 1.0, -1.0]); db = gather_blocks(d, np.arange(N).reshape(K, N // K)) | |
| 42 | + out, diag = estimate_all(m, db, splits_max=3, q_max=2) | |
| 43 | + c = out['central'] | |
| 44 | + gap = {k: np.sqrt(N) * np.linalg.norm(v - c) for k, v in out.items()} | |
| 45 | + assert gap['r1_q2'] < gap['r1_q1'] < gap['r1_q0'] | |
| 46 | + assert gap['r2_q1'] < gap['r1_q1'] | |
| 47 | + assert gap['r2_q2'] < 1e-6 and gap['r3_q2'] < 1e-5 | |
| 48 | + assert gap['csl2'] < gap['csl1'] | |
| 49 | + assert diag['central_status'] == 0 | |
| 50 | + | |
| 51 | + | |
| 52 | +def test_aee_and_csl_basic(): | |
| 53 | + ex, d, db, N, K = _exp_design(seed=3) | |
| 54 | + fits = local_split_fits(ex, db, (1,)) | |
| 55 | + th = fits[1][0][:, 0]; H = fits[1][2][:, 0] | |
| 56 | + a = aee_lin_xi(th, H) | |
| 57 | + assert np.isfinite(a).all() | |
| 58 | + path = csl(ex, db, th.mean(axis=0), rounds=1) | |
| 59 | + # exponential model: L_0 and L_N differ by a linear term, so the full CSL surrogate maximizer IS the centralized MLE | |
| 60 | + assert abs(path[1][0] - 1 / d['x'].mean()) < 1e-10 | |
| 61 | + # the quadratic (one-step) CSL variant with machine 0's Hessian equals the exact Newton step here (data-free Hessian) | |
| 62 | + single, _ = global_newton(ex, db, th.mean(axis=0), q=1, hessian_mode="single") | |
| 63 | + exact, _ = global_newton(ex, db, th.mean(axis=0), q=1) | |
| 64 | + assert abs(single[1][0] - exact[1][0]) < 1e-12 | |
| 65 | + | |
| 66 | + | |
| 67 | +def test_analytic_bias_coefficient_exponential(): | |
| 68 | + ex, d, db, N, K = _exp_design(seed=4) | |
| 69 | + th = local_split_fits(ex, db, (1,))[1][0][:, 0] | |
| 70 | + thc, b = analytic_bias_correction(ex, db, th) | |
| 71 | + assert abs(b.mean() - 1.0) < 0.05 # b_1 = theta0 = 1 | |
| 72 | + | |
| 73 | + | |
| 74 | +def test_quantize_relative(): | |
| 75 | + A = np.array([[1.234567, -0.001234], [3.0, 1e-9]]) | |
| 76 | + Q = quantize_relative(A, 8) | |
| 77 | + assert np.all(np.abs(Q - A) <= np.abs(A) * 2.0 ** -8 + 1e-300) | |
| 78 | + assert np.allclose(quantize_relative(A, 52), A) | |
added
tests/test_jackknife_weights.py
+43 −0
@@ -0,0 +1,43 @@ | ||
| 1 | +import numpy as np | |
| 2 | +from fractions import Fraction | |
| 3 | +from pamle.jackknife import (vandermonde_weights_exact, vandermonde_weights, leading_bias_coefficient, | |
| 4 | + amplification_constant, vandermonde_condition_number, combine_levels, DEFAULT_SPLITS) | |
| 5 | + | |
| 6 | + | |
| 7 | +def test_moment_constraints_exact(): | |
| 8 | + for s in [(1, 2), (1, 2, 4), (1, 2, 4, 8), (1, 2, 3), (1, 3, 9), (1, 2, 4, 8, 16)]: | |
| 9 | + w = vandermonde_weights_exact(s) | |
| 10 | + assert sum(w) == 1 | |
| 11 | + for a in range(1, len(s)): | |
| 12 | + assert sum(wj * Fraction(sj) ** a for wj, sj in zip(w, s)) == 0 | |
| 13 | + | |
| 14 | + | |
| 15 | +def test_known_weights(): | |
| 16 | + assert vandermonde_weights_exact((1, 2)) == [2, -1] | |
| 17 | + assert vandermonde_weights_exact((1, 2, 4)) == [Fraction(8, 3), -2, Fraction(1, 3)] | |
| 18 | + assert np.allclose(vandermonde_weights((1, 2, 4)), [8 / 3, -2, 1 / 3]) | |
| 19 | + | |
| 20 | + | |
| 21 | +def test_float_moment_constraints(): | |
| 22 | + for r in range(1, 6): | |
| 23 | + s = DEFAULT_SPLITS[r]; w = vandermonde_weights(s); sv = np.array(s, float) | |
| 24 | + assert abs(w.sum() - 1) < 1e-12 | |
| 25 | + for a in range(1, r): | |
| 26 | + assert abs(np.sum(w * sv ** a)) < 1e-9 * np.sum(np.abs(w) * sv ** a) | |
| 27 | + | |
| 28 | + | |
| 29 | +def test_leading_coefficients_and_amplification(): | |
| 30 | + assert leading_bias_coefficient((1, 2)) == -2 | |
| 31 | + assert leading_bias_coefficient((1, 2, 4)) == 8 | |
| 32 | + assert leading_bias_coefficient((1, 2, 4, 8)) == -64 | |
| 33 | + assert amplification_constant((1, 2)) == 4 and amplification_constant((1, 2, 4)) == 8 | |
| 34 | + assert vandermonde_condition_number((1, 2, 4, 8)) > vandermonde_condition_number((1, 2)) | |
| 35 | + | |
| 36 | + | |
| 37 | +def test_combine_levels_preserves_mean(): | |
| 38 | + rng = np.random.default_rng(0) | |
| 39 | + K, p = 5, 3 | |
| 40 | + common = rng.standard_normal((K, p)) | |
| 41 | + levels = np.stack([common, common, common]) # identical level means -> J equals them (sum w = 1) | |
| 42 | + J = combine_levels(levels, vandermonde_weights((1, 2, 4))) | |
| 43 | + assert np.allclose(J, common) | |
added
tests/test_math_identities.py
+81 −0
@@ -0,0 +1,81 @@ | ||
| 1 | +"""Exact exponential-model identities (Theorem 3 ingredients) checked numerically on random inputs.""" | |
| 2 | +import numpy as np | |
| 3 | +from scipy import stats | |
| 4 | +from pamle.jackknife import vandermonde_weights, leading_bias_coefficient | |
| 5 | + | |
| 6 | + | |
| 7 | +def test_gamma_inverse_moments(): | |
| 8 | + rng = np.random.default_rng(1) | |
| 9 | + for n in (5, 17, 64): | |
| 10 | + theta0 = 2.5 | |
| 11 | + S = rng.gamma(n, 1 / theta0, size=400000) | |
| 12 | + th = n / S | |
| 13 | + assert abs(th.mean() - theta0 * n / (n - 1)) < 0.01 * theta0 | |
| 14 | + assert abs(th.var() - theta0 ** 2 * n ** 2 / ((n - 1) ** 2 * (n - 2))) < 0.05 * theta0 ** 2 * n ** 2 / ((n - 1) ** 2 * (n - 2)) | |
| 15 | + | |
| 16 | + | |
| 17 | +def test_newton_identity(): | |
| 18 | + rng = np.random.default_rng(2) | |
| 19 | + for _ in range(1000): | |
| 20 | + x = rng.uniform(0.1, 5); xbar = rng.uniform(0.1, 5); th = 1 / xbar | |
| 21 | + g = 1 / x - xbar; H = -1 / x ** 2 | |
| 22 | + newton = x - g / H | |
| 23 | + assert abs(newton - (2 * x - x * x * xbar)) < 1e-12 | |
| 24 | + assert abs((newton - th) + (x - th) ** 2 / th) < 1e-12 | |
| 25 | + | |
| 26 | + | |
| 27 | +def test_newton_basin(): | |
| 28 | + th = 1.3 | |
| 29 | + for x in (0.01, 0.5, 1.0, 2.0, 2.59): | |
| 30 | + e = x - th | |
| 31 | + for _ in range(60): | |
| 32 | + e = -e * e / th | |
| 33 | + assert abs(e) < 1e-6 | |
| 34 | + for x in (-0.1, 2.61, 3.0): | |
| 35 | + e = x - th | |
| 36 | + for _ in range(12): | |
| 37 | + e = -e * e / th | |
| 38 | + assert abs(e) > 1 | |
| 39 | + e = 0.0 + (2 * th - th) # x = 2 theta_hat -> maps to 0 -> stays at 0 (error -theta_hat) | |
| 40 | + e = -e * e / th | |
| 41 | + assert abs(e + th) < 1e-12 | |
| 42 | + | |
| 43 | + | |
| 44 | +def test_exact_bias_of_split_estimators(): | |
| 45 | + """E[J_r] - theta0 = theta0 * sum_j w_j s_j/(m - s_j) exactly, leading term beta_r theta0 / m^r.""" | |
| 46 | + theta0 = 1.0 | |
| 47 | + for s in [(1, 2), (1, 2, 4), (1, 2, 4, 8)]: | |
| 48 | + w = vandermonde_weights(s); r = len(s) | |
| 49 | + for m in (64, 256, 1024): | |
| 50 | + bias = sum(wj * theta0 * m / (m - sj) for wj, sj in zip(w, s)) - theta0 | |
| 51 | + if s == (1, 2): | |
| 52 | + assert abs(bias + 2 * theta0 / ((m - 1) * (m - 2))) < 1e-14 | |
| 53 | + if s == (1, 2, 4): | |
| 54 | + assert abs(bias - 8 * theta0 / ((m - 1) * (m - 2) * (m - 4))) < 1e-14 | |
| 55 | + assert abs(bias * m ** r / (float(leading_bias_coefficient(s)) * theta0) - 1) < 20 / m | |
| 56 | + | |
| 57 | + | |
| 58 | +def test_boundary_constant_numeric(): | |
| 59 | + """Deterministic version of Theorem 3: plug e0 = beta_r theta0 c^r N^{-1/2^(q+1)} into the q-fold Newton identity.""" | |
| 60 | + theta0, c = 1.7, 0.8 | |
| 61 | + for s, q in [((1, 2), 1), ((1, 2, 4), 1), ((1, 2), 2), ((1, 2, 4, 8), 1)]: | |
| 62 | + r = len(s); beta = float(leading_bias_coefficient(s)) | |
| 63 | + for N in (1e12, 1e16): | |
| 64 | + alpha = 1 - 1 / (2 ** (q + 1) * r); m = N ** (1 - alpha) / c | |
| 65 | + e = beta * theta0 * m ** (-r) | |
| 66 | + for _ in range(q): | |
| 67 | + e = -e * e / theta0 | |
| 68 | + lim = -theta0 * (beta * c ** r) ** (2 ** q) | |
| 69 | + assert abs(np.sqrt(N) * e - lim) < 1e-8 * abs(lim) | |
| 70 | + assert abs((-theta0 * ((-2) * c ** 2) ** 2) - (-4 * theta0 * c ** 4)) < 1e-14 | |
| 71 | + | |
| 72 | + | |
| 73 | +def test_exact_variance_of_half_block_jackknife(): | |
| 74 | + """Var(J2) = theta0^2 m^2 (m^2 - 6m + 17)/((m-4)(m-2)^2(m-1)^2) (SymPy, theory/symbolic_checks.py S6) vs simulation.""" | |
| 75 | + rng = np.random.default_rng(5); m = 32; B = 400000 | |
| 76 | + halves = rng.gamma(m // 2, 1.0, size=(B, 2)) | |
| 77 | + full = m / halves.sum(axis=1); half = (m / 2) / halves | |
| 78 | + J2 = 2 * full - half.mean(axis=1) | |
| 79 | + exact = m ** 2 * (m ** 2 - 6 * m + 17) / ((m - 4) * (m - 2) ** 2 * (m - 1) ** 2) | |
| 80 | + assert abs(J2.var() / exact - 1) < 0.02 | |
| 81 | + assert abs(J2.mean() - (1 - 2 / ((m - 1) * (m - 2)))) < 0.002 | |
added
tests/test_models_and_gradients.py
+67 −0
@@ -0,0 +1,67 @@ | ||
| 1 | +import numpy as np | |
| 2 | +import pytest | |
| 3 | +from pamle.likelihoods import Exponential, Logistic, Poisson, gather_blocks, STATUS_OK, STATUS_DIVERGED | |
| 4 | + | |
| 5 | + | |
| 6 | +def fd_check(model, theta, db, eps=1e-6): | |
| 7 | + p = theta.shape[1] | |
| 8 | + g, H = model.score_hessian(theta, db) | |
| 9 | + gfd = np.array([(model.loglik(theta + eps * np.eye(p)[i], db) - model.loglik(theta - eps * np.eye(p)[i], db))[0] / (2 * eps) for i in range(p)]) | |
| 10 | + Hfd = np.array([(model.score_hessian(theta + eps * np.eye(p)[i], db)[0] - model.score_hessian(theta - eps * np.eye(p)[i], db)[0])[0] / (2 * eps) for i in range(p)]) | |
| 11 | + S = np.linalg.inv(-H) | |
| 12 | + T = model.third_contraction(theta, db, S) | |
| 13 | + Tfd = np.array([np.sum((model.score_hessian(theta + eps * np.eye(p)[i], db)[1] - model.score_hessian(theta - eps * np.eye(p)[i], db)[1])[0] / (2 * eps) * S[0]) for i in range(p)]) | |
| 14 | + return np.max(np.abs(g[0] - gfd)), np.max(np.abs(H[0] - Hfd)), np.max(np.abs(T[0] - Tfd)) | |
| 15 | + | |
| 16 | + | |
| 17 | +@pytest.mark.parametrize("model,theta0", [(Logistic(4), [0.3, 0.5, -0.7, 0.2]), (Poisson(3), [0.4, 0.3, -0.2]), (Exponential(), [1.7])]) | |
| 18 | +def test_derivatives_against_finite_differences(model, theta0): | |
| 19 | + rng = np.random.default_rng(11) | |
| 20 | + d = model.simulate(rng, 4000, theta0) | |
| 21 | + db = {k: v[None] for k, v in d.items()} | |
| 22 | + th = np.asarray(theta0, float)[None] * 0.8 + 0.05 | |
| 23 | + eg, eh, et = fd_check(model, th, db) | |
| 24 | + assert eg < 1e-8 and eh < 1e-7 and et < 1e-6 | |
| 25 | + | |
| 26 | + | |
| 27 | +def test_fit_matches_closed_form_exponential(): | |
| 28 | + rng = np.random.default_rng(3); ex = Exponential() | |
| 29 | + d = ex.simulate(rng, 500, 2.0) | |
| 30 | + th, info = ex.fit(d) | |
| 31 | + assert abs(th[0] - 1 / d['x'].mean()) < 1e-12 and info['status'] == STATUS_OK | |
| 32 | + | |
| 33 | + | |
| 34 | +def test_batched_fit_equals_individual_fits(): | |
| 35 | + rng = np.random.default_rng(4); m = Logistic(3) | |
| 36 | + d = m.simulate(rng, 2048, [0.2, 0.8, -0.5]) | |
| 37 | + db = gather_blocks(d, np.arange(2048).reshape(8, 256)) | |
| 38 | + th, info = m.fit_batched(db, tol=1e-11) | |
| 39 | + assert np.all(info['status'] == STATUS_OK) | |
| 40 | + for k in range(8): | |
| 41 | + thk, _ = m.fit({kk: v[k] for kk, v in db.items()}, tol=1e-11) | |
| 42 | + assert np.allclose(th[k], thk, atol=1e-8) | |
| 43 | + | |
| 44 | + | |
| 45 | +def test_score_zero_at_mle_tight(): | |
| 46 | + rng = np.random.default_rng(5); m = Poisson(3) | |
| 47 | + d = m.simulate(rng, 5000, [0.5, 0.3, -0.2]) | |
| 48 | + th, info = m.fit(d, tol=1e-12) | |
| 49 | + assert info['grad_max'] < 1e-12 | |
| 50 | + | |
| 51 | + | |
| 52 | +def test_separation_is_flagged(): | |
| 53 | + """Perfectly separated logistic sample: the MLE does not exist; status must not be OK.""" | |
| 54 | + m = Logistic(2) | |
| 55 | + X = np.hstack([np.ones((40, 1)), np.linspace(-1, 1, 40)[:, None]]) | |
| 56 | + y = (X[:, 1] > 0).astype(float) | |
| 57 | + th, info = m.fit({'X': X, 'y': y}, maxit=200) | |
| 58 | + assert info['status'] != STATUS_OK | |
| 59 | + | |
| 60 | + | |
| 61 | +def test_random_blocks_cover_everything(): | |
| 62 | + from pamle.partitioning import random_blocks, equal_blocks, skewed_sizes, unequal_blocks | |
| 63 | + rng = np.random.default_rng(0) | |
| 64 | + b = random_blocks(1000, 10, rng); assert sorted(b.ravel().tolist()) == list(range(1000)) | |
| 65 | + assert equal_blocks(12, 3).shape == (3, 4) | |
| 66 | + s = skewed_sizes(1000, 7, 4.0); assert s.sum() == 1000 and s.max() / s.min() > 3 | |
| 67 | + ub = unequal_blocks(s); assert sum(len(u) for u in ub) == 1000 | |
added
tests/test_reproducibility.py
+46 −0
@@ -0,0 +1,46 @@ | ||
| 1 | +"""Seeded runs are bit-reproducible, and the library reproduces the prior draft's exact-Gamma Monte Carlo numbers.""" | |
| 2 | +import csv, pathlib | |
| 3 | +import numpy as np | |
| 4 | +from pamle.jackknife import vandermonde_weights | |
| 5 | + | |
| 6 | +ROOT = pathlib.Path(__file__).resolve().parents[1] | |
| 7 | + | |
| 8 | + | |
| 9 | +def _prior_summary(): | |
| 10 | + f = ROOT / "references" / "initial_pamle" / "unpacked" / "simulation_summary.csv" | |
| 11 | + rows = list(csv.DictReader(open(f))) | |
| 12 | + return {(int(r['m']), r['method']): float(r['mean_scaled_error']) for r in rows} | |
| 13 | + | |
| 14 | + | |
| 15 | +def _run_prior_design(reps=500, seed=20261028, N=2 ** 20, m_list=(256, 64, 16)): | |
| 16 | + rng = np.random.default_rng(seed); w2 = vandermonde_weights((1, 2)); w3 = vandermonde_weights((1, 2, 4)) | |
| 17 | + out = {} | |
| 18 | + for m in m_list: | |
| 19 | + K = N // m; acc = {k: [] for k in ('raw', 'one_step_raw', 'two_step_raw', 'jack2', 'one_step_jack2', 'jack3', 'one_step_jack3')} | |
| 20 | + for _ in range(reps): | |
| 21 | + quarters = rng.gamma(shape=m // 4, scale=1.0, size=(K, 4)) | |
| 22 | + halves = quarters.reshape(K, 2, 2).sum(axis=2); totals = halves.sum(axis=1) | |
| 23 | + full_mle = N / totals.sum() | |
| 24 | + lf, lh, lq = m / totals, (m / 2) / halves, (m / 4) / quarters | |
| 25 | + raw = lf.mean() | |
| 26 | + J2 = (w2[0] * lf + w2[1] * lh.mean(axis=1)).mean() | |
| 27 | + J3 = (w3[0] * lf + w3[1] * lh.mean(axis=1) + w3[2] * lq.mean(axis=1)).mean() | |
| 28 | + newton = lambda x: 2 * x - x * x / full_mle | |
| 29 | + vals = {'raw': raw, 'one_step_raw': newton(raw), 'two_step_raw': newton(newton(raw)), 'jack2': J2, | |
| 30 | + 'one_step_jack2': newton(J2), 'jack3': J3, 'one_step_jack3': newton(J3)} | |
| 31 | + for k, v in vals.items(): | |
| 32 | + acc[k].append(np.sqrt(N) * (v - full_mle)) | |
| 33 | + for k, v in acc.items(): | |
| 34 | + out[(m, k)] = float(np.mean(v)) | |
| 35 | + return out | |
| 36 | + | |
| 37 | + | |
| 38 | +def test_prior_draft_numbers_are_reproduced(): | |
| 39 | + prior = _prior_summary(); mine = _run_prior_design() | |
| 40 | + for key, val in prior.items(): | |
| 41 | + assert abs(mine[key] - val) <= 1e-9 * max(1.0, abs(val)), (key, mine[key], val) | |
| 42 | + | |
| 43 | + | |
| 44 | +def test_seeded_determinism(): | |
| 45 | + a = _run_prior_design(reps=5, m_list=(64,)); b = _run_prior_design(reps=5, m_list=(64,)) | |
| 46 | + assert a == b | |
added
theory/symbolic_checks.py
+168 −0
@@ -0,0 +1,168 @@ | ||
| 1 | +#!/usr/bin/env python3 | |
| 2 | +"""Symbolic verification of the algebraic identities used in WP28 (SymPy). | |
| 3 | + | |
| 4 | +Every check raises AssertionError on failure. Run: python3 theory/symbolic_checks.py | |
| 5 | +Items (numbering follows theory/THEOREM_AUDIT.md): | |
| 6 | + S1 Vandermonde extrapolation weights for several split patterns; moment constraints. | |
| 7 | + S2 Exponential(rate) model: Gamma inverse moments, E[1/Xbar_n], Var, remainder R_n moments. | |
| 8 | + S3 Exact bias of the split estimators J_r in the exponential model and its leading coefficient beta_r. | |
| 9 | + S4 Exact Newton identity N(x) = 2x - x^2 Xbar and the error identity N(x) - theta_hat = -(x - theta_hat)^2 / theta_hat. | |
| 10 | + S5 q-fold Newton error identity and the boundary limit -theta0 (beta_r c^r)^(2^q). | |
| 11 | + S6 Exact variance of the half-block jackknife J_2 in the exponential model (nested dependence handled exactly). | |
| 12 | + S7 Newton basin of attraction in the exponential model: convergence iff 0 < x < 2 theta_hat. | |
| 13 | + S8 Remainder-amplification constant A_r = sum |w_j| s_j and the sufficient exponent table alpha(r,q). | |
| 14 | +""" | |
| 15 | +import sympy as sp | |
| 16 | +from fractions import Fraction | |
| 17 | + | |
| 18 | +def vandermonde_weights(s): | |
| 19 | + """Weights w solving sum w_j = 1, sum w_j s_j^a = 0 (a=1..r-1); exact rationals.""" | |
| 20 | + r = len(s) | |
| 21 | + M = sp.Matrix([[sp.Integer(sj)**a for sj in s] for a in range(r)]) | |
| 22 | + rhs = sp.Matrix([1] + [0]*(r-1)) | |
| 23 | + return list(M.LUsolve(rhs)) | |
| 24 | + | |
| 25 | +# ---------- S1: weights ---------- | |
| 26 | +checks = {} | |
| 27 | +for s in [(1,2), (1,2,4), (1,2,4,8), (1,2,3), (1,3,9), (1,2,4,8,16)]: | |
| 28 | + w = vandermonde_weights(s) | |
| 29 | + for a in range(len(s)): | |
| 30 | + tot = sum(wj * sp.Integer(sj)**a for wj, sj in zip(w, s)) | |
| 31 | + assert tot == (1 if a == 0 else 0), (s, a, tot) | |
| 32 | + checks[s] = w | |
| 33 | +assert checks[(1,2)] == [2, -1] | |
| 34 | +assert checks[(1,2,4)] == [sp.Rational(8,3), -2, sp.Rational(1,3)] | |
| 35 | +print("S1 OK weights:", {k: [str(x) for x in v] for k, v in checks.items()}) | |
| 36 | + | |
| 37 | +# ---------- S2: exponential model moments ---------- | |
| 38 | +n, th = sp.symbols('n theta0', positive=True) | |
| 39 | +# S_n ~ Gamma(n, rate theta0): E[S^-k] = theta0^k * Gamma(n-k)/Gamma(n) | |
| 40 | +def Einv(k): # E[(theta0*Xbar_n)^(-k)] = E[(S_n/n)^-k] * theta0^-k ... compute E[Y^-k] with Y = theta0*Xbar_n, nY ~ Gamma(n,1) | |
| 41 | + # E[(nY)^-k] = Gamma(n-k)/Gamma(n) => E[Y^-k] = n^k Gamma(n-k)/Gamma(n) | |
| 42 | + return sp.simplify(n**k * sp.gamma(n-k)/sp.gamma(n)) | |
| 43 | +EY1 = sp.simplify(sp.gammasimp(Einv(1))) | |
| 44 | +EY2 = sp.simplify(sp.gammasimp(Einv(2))) | |
| 45 | +assert sp.simplify(EY1 - n/(n-1)) == 0 | |
| 46 | +assert sp.simplify(EY2 - n**2/((n-1)*(n-2))) == 0 | |
| 47 | +# E[theta_hat_n] = theta0 * n/(n-1); Var = theta0^2 n^2/((n-1)^2 (n-2)) | |
| 48 | +var = sp.simplify(th**2*(EY2 - EY1**2)) | |
| 49 | +assert sp.simplify(var - th**2*n**2/((n-1)**2*(n-2))) == 0 | |
| 50 | +# remainder R_n = theta0 (Y^-1 - 2 + Y), E Y = 1, E Y^2 = 1 + 1/n | |
| 51 | +ER = sp.simplify(th*(EY1 - 2 + 1)) | |
| 52 | +assert sp.simplify(ER - th/(n-1)) == 0 | |
| 53 | +ER2 = sp.simplify(th**2*(EY2 - 4*EY1 + 6 - 4*1 + (1 + 1/n))) | |
| 54 | +assert sp.simplify(ER2 - th**2*(3*n+2)/(n*(n-1)*(n-2))) == 0 | |
| 55 | +# series: E R_n = theta0 (1/n + 1/n^2 + ...) => b_a = theta0 for all a | |
| 56 | +ser = sp.series(th/(n-1), n, sp.oo, 4).removeO() | |
| 57 | +assert sp.simplify(ser - th*(1/n + 1/n**2 + 1/n**3)) == 0 | |
| 58 | +# uniform constant: sup_{n>=3} n^2 E R_n^2 / theta0^2 | |
| 59 | +f = sp.lambdify(n, n**2*(3*n+2)/(n*(n-1)*(n-2))) | |
| 60 | +Csup = max(f(k) for k in range(3, 2000)) | |
| 61 | +assert abs(Csup - 16.5) < 1e-12 | |
| 62 | +print("S2 OK E[theta_hat_n]=theta0 n/(n-1), Var=theta0^2 n^2/((n-1)^2(n-2)), E R_n=theta0/(n-1), E R_n^2=theta0^2(3n+2)/(n(n-1)(n-2)), sup n^2 E R_n^2/theta0^2 = 16.5 (n=3)") | |
| 63 | + | |
| 64 | +# ---------- S3: exact bias of J_r and leading coefficient ---------- | |
| 65 | +m = sp.symbols('m', positive=True) | |
| 66 | +def exact_bias(s): | |
| 67 | + w = vandermonde_weights(s) | |
| 68 | + # sub-block size m/s_j, E[theta_hat] = theta0 (m/s_j)/((m/s_j)-1) = theta0 m/(m - s_j) | |
| 69 | + return sp.simplify(sum(wj*th*m/(m-sj) for wj, sj in zip(w, s)) - th) | |
| 70 | +b2 = exact_bias((1,2)); b3 = exact_bias((1,2,4)) | |
| 71 | +assert sp.simplify(b2 + 2*th/((m-1)*(m-2))) == 0 | |
| 72 | +assert sp.simplify(b3 - 8*th/((m-1)*(m-2)*(m-4))) == 0 | |
| 73 | +betas = {} | |
| 74 | +for s in [(1,2),(1,2,4),(1,2,4,8),(1,2,3),(1,3,9)]: | |
| 75 | + w = vandermonde_weights(s); r = len(s) | |
| 76 | + beta = sum(wj*sp.Integer(sj)**r for wj, sj in zip(w, s)) | |
| 77 | + lead = sp.limit(exact_bias(s)*m**r, m, sp.oo) | |
| 78 | + assert sp.simplify(lead - th*beta) == 0, (s, lead, beta) | |
| 79 | + betas[s] = beta | |
| 80 | +assert betas[(1,2)] == -2 and betas[(1,2,4)] == 8 and betas[(1,2,4,8)] == -64 | |
| 81 | +print("S3 OK exact biases: J2 = -2 theta0/((m-1)(m-2)), J3 = 8 theta0/((m-1)(m-2)(m-4)); beta_r:", {k: str(v) for k, v in betas.items()}) | |
| 82 | + | |
| 83 | +# ---------- S4: Newton identity ---------- | |
| 84 | +x, xbar = sp.symbols('x xbar', positive=True) | |
| 85 | +L = sp.log(x) - x*xbar # average log-likelihood of Exp(rate x) | |
| 86 | +g = sp.diff(L, x); H = sp.diff(L, x, 2) | |
| 87 | +Nx = sp.simplify(x - g/H) | |
| 88 | +assert sp.simplify(Nx - (2*x - x**2*xbar)) == 0 | |
| 89 | +that = 1/xbar | |
| 90 | +assert sp.simplify((Nx - that) + (x - that)**2/that) == 0 | |
| 91 | +print("S4 OK N(x) = 2x - x^2 Xbar, N(x) - theta_hat = -(x - theta_hat)^2/theta_hat") | |
| 92 | + | |
| 93 | +# ---------- S5: q-fold identity and boundary limit ---------- | |
| 94 | +e0, thh, c, N = sp.symbols('e0 thetahat c N', positive=True) | |
| 95 | +def newton_err(e, q): | |
| 96 | + for _ in range(q): | |
| 97 | + e = -e**2/thh | |
| 98 | + return sp.simplify(e) | |
| 99 | +for q in range(1, 4): | |
| 100 | + assert sp.simplify(newton_err(e0, q) + e0**(2**q)/thh**(2**q-1)) == 0 | |
| 101 | +# boundary: K = c N^alpha, alpha = 1 - 1/(2^(q+1) r) => m = N/K = N^{1/(2^(q+1) r)}/c ; e0 ~ beta_r theta0 m^{-r} | |
| 102 | +for s in [(1,2),(1,2,4),(1,2,4,8)]: | |
| 103 | + r = len(s); beta = betas[s] | |
| 104 | + for q in range(0, 3): | |
| 105 | + alpha = 1 - sp.Rational(1, 2**(q+1)*r) | |
| 106 | + mm = N**(1-alpha)/c | |
| 107 | + e_pilot = beta*th*mm**(-r) | |
| 108 | + if q == 0: | |
| 109 | + lim = sp.limit(sp.sqrt(N)*e_pilot * N**(sp.Rational(1,2)) / N**(sp.Rational(1,2)), N, sp.oo) # sqrt(N) e0 -> infinity: use N^{1/2} scaling only for q>=1 | |
| 110 | + err = newton_err(e_pilot, q).subs(thh, th) | |
| 111 | + scaled = sp.simplify(sp.sqrt(N)*err) | |
| 112 | + lim = sp.limit(scaled, N, sp.oo) | |
| 113 | + expected = -th*(beta*c**r)**(2**q) if q >= 1 else sp.nan | |
| 114 | + if q >= 1: | |
| 115 | + assert sp.simplify(lim - expected) == 0, (s, q, lim, expected) | |
| 116 | + else: | |
| 117 | + # q = 0: the scaled bias N^{1/2} m^{-r} = c^r N^{1/2 - 1/2} -> beta theta0 c^r (alpha(r,0) = 1 - 1/(2r) makes sqrt(N) m^{-r} = c^r) | |
| 118 | + assert sp.simplify(lim - beta*th*c**r) == 0, (s, q, lim) | |
| 119 | +assert sp.simplify(-th*((-2)*c**2)**2 + 4*th*c**4) == 0 | |
| 120 | +print("S5 OK N^q(x) - theta_hat = -(x-theta_hat)^(2^q)/theta_hat^(2^q-1); boundary limit sqrt(N)(theta^(q) - theta_hat) -> -theta0 (beta_r c^r)^(2^q) for q>=1 (r=2,q=1: -4 theta0 c^4) and beta_r theta0 c^r for q=0") | |
| 121 | + | |
| 122 | +# ---------- S6: exact variance of J2 in the exponential model ---------- | |
| 123 | +# S_A, S_B ~ Gamma(h, 1), h = m/2 (theta0 = 1 w.l.o.g.); S = S_A+S_B ~ Gamma(2h,1), W = S_A/S ~ Beta(h,h) indep. of S. | |
| 124 | +# theta_full = m/S = 2h/S ; theta_A = h/S_A = h/(S W) ; theta_B = h/(S (1-W)). | |
| 125 | +# J2 = 2*theta_full - (theta_A+theta_B)/2 = 4h/S - (h/2)(1/W + 1/(1-W))/S | |
| 126 | +h = sp.symbols('h', positive=True) | |
| 127 | +def EgammaInv(a, k): # E[S^-k], S ~ Gamma(a,1) | |
| 128 | + return sp.gamma(a-k)/sp.gamma(a) | |
| 129 | +def EbetaInv(a, b, k, l): # E[W^-k (1-W)^-l], W ~ Beta(a,b) | |
| 130 | + return sp.gamma(a-k)*sp.gamma(b-l)*sp.gamma(a+b)/(sp.gamma(a)*sp.gamma(b)*sp.gamma(a+b-k-l)) | |
| 131 | +# J2 = (1/S) * [4h - (h/2) (1/W + 1/(1-W))] =: (1/S) * G(W) | |
| 132 | +# E J2 = E[S^-1] * E[G]; E J2^2 = E[S^-2] * E[G^2] | |
| 133 | +EG = 4*h - (h/2)*(EbetaInv(h,h,1,0) + EbetaInv(h,h,0,1)) | |
| 134 | +EG2 = 16*h**2 - 2*4*h*(h/2)*(EbetaInv(h,h,1,0)+EbetaInv(h,h,0,1)) + (h**2/4)*(EbetaInv(h,h,2,0) + EbetaInv(h,h,0,2) + 2*EbetaInv(h,h,1,1)) | |
| 135 | +EJ2 = sp.gammasimp(EgammaInv(2*h,1)*EG) | |
| 136 | +EJ22 = sp.gammasimp(EgammaInv(2*h,2)*EG2) | |
| 137 | +EJ2_m = sp.simplify(EJ2.subs(h, m/2)) | |
| 138 | +VarJ2_m = sp.simplify((EJ22 - EJ2**2).subs(h, m/2)) | |
| 139 | +assert sp.simplify(EJ2_m - (1 - 2/((m-1)*(m-2)))) == 0 # matches S3 (theta0 = 1) | |
| 140 | +# asymptotic expansion of the variance: 1/m + v2/m^2 + ... | |
| 141 | +serV = sp.series(VarJ2_m, m, sp.oo, 4).removeO() | |
| 142 | +v1 = sp.limit(VarJ2_m*m, m, sp.oo); v2 = sp.limit((VarJ2_m - 1/m)*m**2, m, sp.oo) | |
| 143 | +assert v1 == 1 | |
| 144 | +# compare with the full-block MLE variance m^2/((m-1)^2 (m-2)) = 1/m + 3/m^2 + ... | |
| 145 | +v2_full = sp.limit((m**2/((m-1)**2*(m-2)) - 1/m)*m**2, m, sp.oo) | |
| 146 | +print("S6 OK Var(J2) exact =", sp.factor(VarJ2_m), " ~ 1/m + (%s)/m^2 ; Var(theta_hat_m) ~ 1/m + (%s)/m^2" % (v2, v2_full)) | |
| 147 | +VARJ2_EXACT = VarJ2_m | |
| 148 | + | |
| 149 | +# ---------- S7: Newton basin ---------- | |
| 150 | +# e_{j+1} = -e_j^2/theta_hat ; |e_{j+1}|/theta_hat = (|e_j|/theta_hat)^2 -> 0 iff |e_0| < theta_hat iff 0 < x < 2 theta_hat | |
| 151 | +u = sp.symbols('u', positive=True) | |
| 152 | +assert sp.simplify(sp.Abs(-(u*thh)**2/thh)/thh - u**2) == 0 | |
| 153 | +print("S7 OK Newton basin: |x - theta_hat| < theta_hat <=> 0 < x < 2 theta_hat ; x <= 0 or x >= 2 theta_hat diverges (x = 2 theta_hat maps to 0 then stays at 0)") | |
| 154 | + | |
| 155 | +# ---------- S8: amplification constants and exponent table ---------- | |
| 156 | +def A_const(s): | |
| 157 | + w = vandermonde_weights(s); return sum(abs(wj)*sj for wj, sj in zip(w, s)) | |
| 158 | +def W_sum_abs(s): | |
| 159 | + w = vandermonde_weights(s); return sum(abs(wj) for wj in w) | |
| 160 | +tab = [] | |
| 161 | +for s in [(1,), (1,2), (1,2,4), (1,2,4,8), (1,2,3), (1,3,9)]: | |
| 162 | + tab.append((s, A_const(s), W_sum_abs(s))) | |
| 163 | +assert A_const((1,2)) == 4 and A_const((1,2,4)) == 8 | |
| 164 | +alpha = {(r,q): sp.Rational(1) - sp.Rational(1, 2**(q+1)*r) for r in (1,2,3,4) for q in (0,1,2,3)} | |
| 165 | +assert alpha[(1,0)] == sp.Rational(1,2) and alpha[(1,1)] == sp.Rational(3,4) and alpha[(2,1)] == sp.Rational(7,8) and alpha[(3,1)] == sp.Rational(11,12) and alpha[(1,2)] == sp.Rational(7,8) and alpha[(2,0)] == sp.Rational(3,4) and alpha[(3,0)] == sp.Rational(5,6) | |
| 166 | +print("S8 OK A_r = sum|w_j| s_j:", [(s, str(a), str(b)) for s, a, b in tab]) | |
| 167 | +print(" alpha(r,q) = 1 - 1/(2^(q+1) r):", {k: str(v) for k, v in alpha.items()}) | |
| 168 | +print("ALL SYMBOLIC CHECKS PASSED") | |
| 169 | ||