1\documentclass[11pt,letterpaper]{article}2\usepackage[T1]{fontenc}3\usepackage[utf8]{inputenc}4\usepackage{lmodern}5\usepackage{microtype}6\usepackage[margin=1.08in,headheight=23pt]{geometry}7\usepackage{amsmath,amssymb,amsthm,mathtools,bm}8\usepackage{booktabs,tabularx,array,multirow}9\usepackage{graphicx}10\usepackage{xcolor}11\usepackage{fancyhdr}12\usepackage{titlesec}13\usepackage{enumitem}14\usepackage{natbib}15\usepackage{hyperref}16\usepackage[nameinlink,noabbrev,capitalise]{cleveref}17\usepackage{listings}18\usepackage{float}19\usepackage{tikz}20\usetikzlibrary{arrows.meta,positioning,calc}21\definecolor{navy}{HTML}{17304A}22\definecolor{teal}{HTML}{216C78}23\definecolor{charcoal}{HTML}{242E35}24\definecolor{grayblue}{HTML}{728291}25\definecolor{pale}{HTML}{E8F0F2}26\definecolor{rulegray}{HTML}{D6E0E5}27\hypersetup{colorlinks=true,linkcolor=navy,citecolor=teal,urlcolor=teal,pdftitle={Local Bias Correction versus Global Synchronization in Distributed Maximum Likelihood},pdfauthor={Simon-Pierre Boucher and GPT-6.1}}28\titleformat{\section}{\large\bfseries\color{navy}}{\thesection.}{0.6em}{}29\titleformat{\subsection}{\normalsize\bfseries\color{navy}}{\thesubsection}{0.6em}{}30\titleformat{\subsubsection}{\normalsize\bfseries\color{charcoal}}{\thesubsubsection}{0.6em}{}31\titlespacing*{\section}{0pt}{19pt}{8pt}32\titlespacing*{\subsection}{0pt}{12pt}{5pt}33\setlength{\parskip}{4.2pt}34\setcounter{tocdepth}{1}35\setlength{\parindent}{0pt}36\setlength{\emergencystretch}{3em}37\setlist[itemize]{topsep=3pt,itemsep=2pt,leftmargin=1.4em}38\setlist[enumerate]{topsep=3pt,itemsep=3pt,leftmargin=1.6em}39\pagestyle{fancy}40\fancyhf{}41\fancyhead[L]{\small\color{grayblue}\textsc{Boucher \& GPT-6.1}}42\fancyhead[R]{\small\color{grayblue}\textsc{Distributed Maximum Likelihood}}43\fancyfoot[C]{\color{grayblue}\thepage}44\renewcommand{\headrulewidth}{0.35pt}45\renewcommand{\headrule}{\hbox to\headwidth{\color{rulegray}\leaders\hrule height \headrulewidth\hfill}}46\newtheorem{theorem}{Theorem}[section]47\newtheorem{proposition}[theorem]{Proposition}48\newtheorem{corollary}[theorem]{Corollary}49\newtheorem{lemma}[theorem]{Lemma}50\theoremstyle{definition}51\newtheorem{assumption}[theorem]{Assumption}52\newtheorem{remark}[theorem]{Remark}53\newtheorem{example}[theorem]{Example}54\DeclareMathOperator{\Var}{Var}55\DeclareMathOperator{\Bias}{Bias}56\DeclareMathOperator{\MSE}{MSE}57\DeclareMathOperator{\RMSE}{RMSE}58\newcommand{\EE}{\mathbb E}59\newcommand{\PP}{\mathbb P}60\newcommand{\RR}{\mathbb R}61\newcommand{\op}{O_{\!p}}62\newcommand{\smallop}{o_{\!p}}63\newcommand{\thet}{\theta_0}64\newcommand{\thn}{\widehat\theta_N}65\newcommand{\score}{s}66\newcommand{\I}{\mathcal I}67\newcommand{\D}{\mathcal D}68\lstdefinestyle{clean}{basicstyle=\ttfamily\footnotesize,backgroundcolor=\color{pale!25},frame=single,rulecolor=\color{rulegray},keywordstyle=\color{navy}\bfseries,commentstyle=\color{grayblue},showstringspaces=false,breaklines=true,tabsize=2,xleftmargin=1em,xrightmargin=1em}69\lstset{style=clean}70\begin{document}71\thispagestyle{empty}72\vspace*{0.23in}73{\color{teal}\rule{1.35in}{3pt}}\\[22pt]74{\LARGE\bfseries\color{navy} Local Bias Correction versus Global Synchronization in Distributed Maximum Likelihood\par}75\vspace{11pt}76{\Large\color{charcoal} Rate Frontiers, Exact Exponential Analysis, and Reproducible Monte Carlo Evidence\par}77\vspace{17pt}78{\large\bfseries Simon-Pierre Boucher\textsuperscript{1} \quad and \quad GPT-6.1\textsuperscript{2}\par}79\vspace{7pt}80{\small\color{grayblue}\textsuperscript{1} Universit\'e du Qu\'ebec en Outaouais \qquad \textsuperscript{2} AI co-author designation\par}81\vspace{10pt}82{\small\color{grayblue}Working paper \;\textbullet\; October 28, 2026 \;\textbullet\; Theoretical and numerical research note\par}83\vspace{18pt}84{\color{rulegray}\hrule}85\vspace{12pt}86{\bfseries\color{navy} Abstract}\par87Distributed maximum likelihood estimation (MLE) trades statistical efficiency against network communication and local computation. We analyze a two-phase construction that combines a within-worker split-sample jackknife with a single globally aggregated Newton step. Under fixed-dimensional asymptotic linearity, a second-order bias expansion and local Hessian regularity, we prove first-order equivalence to centralized MLE for $K=o(N^{7/8})$ partitions; the proof explicitly retains the jackknife's centered random remainder. A higher-order local extrapolation yields $K=o(N^{11/12})$ under stronger assumptions. For exponential observations, exact bias formulas and a Newton identity establish a nonzero limiting shift at the $N^{7/8}$ boundary of the specified algorithm. Reproducible exact-Gamma Monte Carlo experiments confirm the qualitative pattern. We account for derivative-communication costs and compare extra local estimation with an additional Newton synchronization. These results provide a transparent theoretical and numerical working note, not an established claim of methodological priority or improved wall-clock performance.\par88\vspace{13pt}89{\small\textbf{Keywords:} maximum likelihood; distributed econometrics; split-sample jackknife; one-step Newton; asymptotic equivalence; communication complexity; Monte Carlo.}\par90{\small\textbf{JEL codes:} C13, C15, C55, C63.}\par91\vspace{10pt}92{\color{rulegray}\hrule}93\vspace{10pt}94{\small\color{grayblue}\textbf{Disclosure.} The second author label is an AI designation supplied for this draft. The mathematics is conditional on explicit assumptions, and the numerical experiment is executable and independently reproducible. Researcher review, literature-priority checks, and empirical performance benchmarking remain necessary prior to submission.}\par95\clearpage96\tableofcontents97\clearpage9899\section{Introduction}100Massive datasets distributed across geographically or operationally separate devices create a tension between three objectives: statistically efficient estimation, low network communication, and computational scalability. For a regular parametric model the centralized maximum likelihood estimator can attain the inverse Fisher-information asymptotic variance, but collecting all observations on one machine may be expensive or prohibited. An alternative is to fit the model independently on $K$ data partitions and average the local estimators. This does not automatically preserve centralized statistical efficiency: nonlinear estimation induces a local bias that does not disappear merely by averaging independent estimates.101102Early systematic analyses of divide-and-conquer statistical optimization include \citet{zhang2013}, who study averaging and communication-efficient alternatives, while \citet{huang2019} develop a distributed one-step estimator. Multi-round and adaptive distributed procedures are represented by \citet{fan2023}. Separately, jackknife extrapolation has a long history in bias correction; \citet{dhaene2015} treat split-panel corrections, including higher-order variants, in nonlinear fixed-effects models. Our starting point is to combine these familiar elements and quantify the \emph{exchange rate} between extra local statistical work and an additional global synchronization.103104We focus on a sharply defined question: for fixed parameter dimension $p$, how fast may the partition count $K$ grow with the global sample size $N$ while a locally bias-corrected, globally one-step estimator remains first-order equivalent to the centralized MLE? For the half-block jackknife, we establish $K=o(N^{7/8})$ under explicit regularity conditions. Our proof takes special care to retain a random term of order $1/(m\sqrt K)$, where $m=N/K$ is the local sample size. Omitting that term would produce an incomplete argument, even if the final exponent happened to be correct.105106Our contributions are deliberately bounded. We give a transparent sufficient-rate theorem, an exact verification and boundary distribution for an exponential model, a higher-order generalization, and a reproducible computational illustration. We do \emph{not} establish that this procedure dominates the best available distributed estimators, nor do we prove a universal communication lower bound. The asymptotic regime is fixed-$p$, i.i.d., correctly specified, and locally regular. The exact exponential model actually admits trivial sufficient-statistic aggregation; it is therefore a proof-of-concept and counterexample laboratory, not a favorable application of distributed approximate MLE.107108\section{Model, estimators, and information exchange}109\subsection{Likelihood and partitioning}110Let $X_1,\ldots,X_N$ be i.i.d. observations with density $f(x;\theta)$, $\theta\in\Theta\subseteq\RR^p$, and true value $\thet$ in the interior of $\Theta$. Let $K$ divide $N$, set $m=N/K$, and partition the index set into disjoint blocks $\D_1,\ldots,\D_K$, each of size $m$. The average log-likelihood and centralized MLE are111\begin{align}112 L_N(\theta)&=\frac1N\sum_{i=1}^N\log f(X_i;\theta),113 &\thn&\in\arg\max_{\theta\in\Theta}L_N(\theta).\label{eq:global}114\end{align}115Worker $k$ computes a local MLE $\widehat\theta_{k,m}$ based on $\D_k$. In the baseline construction $m$ is even and each block has two disjoint halves $A_k,B_k$, with half-sample MLEs $\widehat\theta_{k,A}$ and $\widehat\theta_{k,B}$. Write116\begin{align}117 \bar\theta_{\mathrm{raw}}&=K^{-1}\sum_{k=1}^K\widehat\theta_{k,m},\\118 J_k&=2\widehat\theta_{k,m}-\tfrac12\bigl(\widehat\theta_{k,A}+\widehat\theta_{k,B}\bigr),\label{eq:jack}\\119 \theta^{(0)}&=K^{-1}\sum_{k=1}^KJ_k.\label{eq:init}120\end{align}121This local operation needs three likelihood fits per block but, because all three fits access locally stored observations, introduces no within-cluster data transmission beyond the final parameter vector.122123Define the full score and Hessian by $g_N(\theta)=\nabla_\theta L_N(\theta)$ and $H_N(\theta)=\nabla_\theta^2L_N(\theta)$. After the coordinator broadcasts $\theta^{(0)}$, each worker sends its contributions to $g_N(\theta^{(0)})$ and $H_N(\theta^{(0)})$. The coordinator computes124\begin{equation}125 \boxed{\quad\theta^{(1)}=\theta^{(0)}-H_N(\theta^{(0)})^{-1}g_N(\theta^{(0)})\quad.}\label{eq:newton}126\end{equation}127An invertible Hessian and well-defined, consistent local solutions are assumed throughout; a production implementation may require regularization and step damping, neither of which is covered by the present theory.128129\begin{figure}[ht]130\centering131\begin{tikzpicture}[node distance=8mm and 7mm,>=Latex,scale=.95,transform shape,132 box/.style={draw=rulegray,rounded corners=3pt,fill=white,align=center,minimum height=0.65cm,font=\small,inner sep=5pt},133 cent/.style={box,fill=pale!75},arr/.style={->,thick,draw=teal}]134\node[box,text width=10.5cm] (data) {Partitioned observations: $N=Km$ and fixed $p$};135\node[box,below left=of data,xshift=17mm,text width=3.1cm] (w1) {Worker $1$\\full / half-sample MLEs};136\node[box,right=of w1,text width=3.1cm] (w2) {Worker $2$\\full / half-sample MLEs};137\node[box,right=of w2,text width=3.1cm] (wk) {Worker $K$\\full / half-sample MLEs};138\node[cent,below=13mm of w2,text width=8.2cm] (avg) {Coordinator: average jackknifed estimates $\theta^{(0)}$};139\node[cent,below=10mm of avg,text width=8.2cm] (agg) {Broadcast $\theta^{(0)}$; aggregate scores and Hessians};140\node[box,below=9mm of agg,text width=8.2cm] (out) {One Newton correction: $\theta^{(1)}$};141\draw[arr] (data.south) -- (w1.north);\draw[arr] (data.south) -- (w2.north);\draw[arr] (data.south) -- (wk.north);142\draw[arr] (w1.south) -- (avg.north west);\draw[arr] (w2.south) -- (avg.north);\draw[arr] (wk.south) -- (avg.north east);143\draw[arr] (avg) -- (agg);\draw[arr] (agg) -- (out);144\end{tikzpicture}145\caption{Local bias correction followed by one globally aggregated Newton step. Arrows indicate logical phases, not measured network latency.}\label{fig:architecture}146\end{figure}147148\subsection{Regularity and the benchmark MLE}149Our central statistical result is conditional upon the following strong but explicit assumptions, stated for a \emph{sequence} with $N\to\infty$ and $K=K_N\to\infty$ permitted.150151\begin{assumption}[Uniform asymptotic linear expansion]\label{ass:lin}152For every relevant sub-sample size $n\in\{m,m/2,N\}$, its MLE has the representation153\begin{equation}154 \widehat\theta_n-\thet=\frac1n\sum_{i=1}^n\psi(X_i)+R_n,\label{eq:lin}155\end{equation}156where $\EE\psi(X_i)=0$, $\EE\|\psi(X_i)\|^2<\infty$, $\Var(\psi(X_i))=\I(\thet)^{-1}$ for a positive-definite Fisher information matrix, and a common constant $C<\infty$ satisfies157\begin{equation}158 \EE\|R_n\|^2\le Cn^{-2}.\label{eq:r2}159\end{equation}160For correctly specified, regular parametric MLE, $\psi(x)=\I(\thet)^{-1}\nabla_\theta\log f(x;\thet)$. The expansion is required jointly with the same $\psi$ on all nested subsamples.161\end{assumption}162163\begin{assumption}[Second-order bias expansion]\label{ass:bias}164There exist vectors $b_1,b_2$ and a constant $C_b$, not depending on $n$, for which165\begin{equation}166 \EE R_n=\frac{b_1}{n}+\frac{b_2}{n^2}+\rho_n,167 \qquad\|\rho_n\|\le C_b n^{-3},\label{eq:biasexp}168\end{equation}169for the sample sizes involved. The existence of this expansion is not implied by mere consistency or asymptotic normality.170\end{assumption}171172\begin{assumption}[Local Newton regularity]\label{ass:hess}173With probability tending to one, $\widehat\theta_N$ is a consistent interior stationary maximizer, $g_N(\thn)=0$, and there is a fixed convex neighborhood $\mathcal U$ of $\thet$ such that174\begin{equation}175 \sup_{\theta\in\mathcal U}\|H_N(\theta)^{-1}\|=\op(1),\qquad176 \sup_{u\ne v\in\mathcal U}\frac{\|H_N(u)-H_N(v)\|}{\|u-v\|}=\op(1).\label{eq:localreg}177\end{equation}178The probability that both $\thn$ and the initializer lie in $\mathcal U$ tends to one. The centralized MLE obeys $\sqrt N(\thn-\thet)\Rightarrow\mathcal N(0,\I(\thet)^{-1})$.179\end{assumption}180181\begin{remark}182The assumptions deliberately fix $p$ and require i.i.d. partitions, finite moments of inverse-information-type quantities, and a well-behaved local optimizer. Weak identification, incidental parameters, dependent samples, and growing $p$ require different analysis. The theorem does not apply merely because code can execute on multiple machines.183\end{remark}184185\section{An exact decomposition of local bias correction}186\begin{proposition}[Half-block jackknife initializer]\label{prop:init}187Under \cref{ass:lin,ass:bias}, with independent blocks and $m=N/K$ even,188\begin{align}189 \EE[J_k]-\thet&=-\frac{2b_2}{m^2}+O(m^{-3}),\label{eq:jackbias}\\190 \theta^{(0)}-\thn&=191 \op\left(m^{-2}+\frac1{m\sqrt K}+N^{-1}\right)\\192 &=\op\left(\frac{K^2}{N^2}+\frac{\sqrt K}{N}+\frac1N\right).\label{eq:initrate}193\end{align}194\end{proposition}195\begin{proof}196Write $U_k=m^{-1}\sum_{i\in\D_k}\psi(X_i)$ and define $U_{k,A}$ and $U_{k,B}$ analogously. Because the halves partition $\D_k$ exactly,197\begin{equation}198 \tfrac12(U_{k,A}+U_{k,B})=U_k.\label{eq:influencecancel}199\end{equation}200Applying \cref{eq:lin} to the three local estimators yields201\begin{equation}202 J_k-\thet=U_k+R_k^J,\qquad203 R_k^J=2R_{k,m}-\tfrac12(R_{k,A,m/2}+R_{k,B,m/2}).\label{eq:rj}204\end{equation}205Using \cref{eq:biasexp} at $m$ and $m/2$ gives206\begin{align*}207 \EE R_k^J208 &=2\left(\frac{b_1}{m}+\frac{b_2}{m^2}\right)209 -\left(\frac{2b_1}{m}+\frac{4b_2}{m^2}\right)+O(m^{-3})\\210 &=-\frac{2b_2}{m^2}+O(m^{-3}).211\end{align*}212The $L^2$ triangle inequality and \cref{eq:r2} imply213\begin{align*}214 (\EE\|R_k^J\|^2)^{1/2}215 &\le 2\sqrt C/m+\frac12\bigl(2\sqrt C/m+2\sqrt C/m\bigr)\\216 &=4\sqrt C/m.217\end{align*}218Hence $\EE\|R_k^J\|^2\le16C/m^2$. Independence of blocks eliminates cross terms for centered remainders:219\begin{equation}220 \EE\left\|K^{-1}\sum_{k=1}^K(R_k^J-\EE R_k^J)\right\|^2221 \le\frac{16C}{Km^2}.222\end{equation}223Thus, by Markov's inequality,224\begin{equation}225 \theta^{(0)}-\thet226 =\frac1N\sum_{i=1}^N\psi(X_i)227 -\frac{2b_2}{m^2}+O(m^{-3})228 +\op\left(\frac1{m\sqrt K}\right).229\end{equation}230The full-sample expansion gives $\thn-\thet=N^{-1}\sum_i\psi(X_i)+\op(N^{-1})$. Subtraction cancels the leading influence term \emph{exactly}, proving the second assertion.\end{proof}231232\begin{remark}[Why the random remainder matters]\label{rem:noise}233A small \emph{mean} remainder does not imply a small realized remainder. The term $1/(m\sqrt K)=\sqrt K/N$ in \cref{eq:initrate} describes variation of centered local second-order terms. Bias cancellation alone would not justify removing it. Averaging over workers does remove part of its magnitude, but only after an explicit second-moment argument.234\end{remark}235236\section{One Newton step and the $7/8$ frontier}237\begin{lemma}[Quadratic local error contraction]\label{lem:newton}238Under \cref{ass:hess}, for any initializer $\vartheta_N$ satisfying $\vartheta_N\to_p\thet$ and lying with $\thn$ in the regularity neighborhood, its undamped Newton update $\mathcal N_N(\vartheta_N)=\vartheta_N-H_N(\vartheta_N)^{-1}g_N(\vartheta_N)$ obeys239\begin{equation}240 \|\mathcal N_N(\vartheta_N)-\thn\|241 =\op(\|\vartheta_N-\thn\|^2).\label{eq:quadratic}242\end{equation}243\end{lemma}244\begin{proof}245Let $e=\vartheta_N-\thn$. By the fundamental theorem of calculus, $g_N(\vartheta_N)=\int_0^1 H_N(\thn+te)e\,dt$. Since $g_N(\thn)=0$, one obtains246\begin{align*}247 \mathcal N_N(\vartheta_N)-\thn248 &=H_N(\vartheta_N)^{-1}\int_0^1249 [H_N(\vartheta_N)-H_N(\thn+te)]e\,dt.250\end{align*}251The local Lipschitz property bounds the integrand by $\Lambda_N(1-t)\|e\|^2$. Taking norms and integrating yields252\begin{equation*}253 \|\mathcal N_N(\vartheta_N)-\thn\|254 \le\tfrac12\|H_N(\vartheta_N)^{-1}\|\Lambda_N\|e\|^2,255\end{equation*}256where the multiplicative factor is $\op(1)$ by \cref{ass:hess}.\end{proof}257258\begin{theorem}[Jackknife--Newton equivalence]\label{thm:seven}259Under \cref{ass:lin,ass:bias,ass:hess}, with $K=o(N^{7/8})$ and $m=N/K$ even, the estimator in \cref{eq:jack,eq:init,eq:newton} satisfies260\begin{equation}261 \boxed{\quad\sqrt N\bigl(\theta^{(1)}-\thn\bigr)\xrightarrow{p}0.\quad}\label{eq:equiv}262\end{equation}263Consequently,264\begin{equation}265 \sqrt N(\theta^{(1)}-\thet)\Rightarrow266 \mathcal N\left(0,\I(\thet)^{-1}\right).267\end{equation}268\end{theorem}269\begin{proof}270By \cref{prop:init}, the starting discrepancy is271\begin{equation*}272 e_0=\theta^{(0)}-\thn273 =\op(m^{-2}+m^{-1}K^{-1/2}+N^{-1}).274\end{equation*}275The assumed $K$ regime implies $m\to\infty$ and hence consistency of the initializer, permitting \cref{lem:newton}. Then276\begin{align}277 \|\theta^{(1)}-\thn\|278 &=\op\left(m^{-4}+\frac1{Km^2}+N^{-2}\right)\nonumber\\279 &=\op\left(\frac{K^4}{N^4}+\frac K{N^2}+\frac1{N^2}\right).\label{eq:newtonrate}280\end{align}281Multiplication by $\sqrt N$ gives282\begin{equation}283 \sqrt N\|\theta^{(1)}-\thn\|284 =\op\left(\frac{K^4}{N^{7/2}}+\frac K{N^{3/2}}+N^{-3/2}\right).285\end{equation}286All three deterministic orders vanish whenever $K=o(N^{7/8})$. The first result follows; the limiting distribution then follows from the centralized MLE central limit theorem and Slutsky's lemma.\end{proof}287288\begin{corollary}[Comparison with an uncorrected start]\label{cor:raw}289If the ordinary local average instead has bias $\op(m^{-1})$ and a linear expansion of the form used above, then one global Newton step has a sufficient statistical-equivalence regime $K=o(N^{3/4})$. Two global Newton steps, started from the same raw average, have a sufficient regime $K=o(N^{7/8})$ under analogous local regularity at successive steps.290\end{corollary}291\begin{proof}292The raw starting discrepancy is $\op(m^{-1}+m^{-1}K^{-1/2}+N^{-1})$ and is dominated by $m^{-1}$ in the relevant regime. One squaring requires $\sqrt N m^{-2}\to0$, or $K=o(N^{3/4})$; two squarings require $\sqrt Nm^{-4}\to0$, or $K=o(N^{7/8})$. All stochastic remainders are negligible under these regimes.\end{proof}293294\begin{remark}[Interpretation rather than priority claim]295\cref{thm:seven} does not by itself establish a new statistical estimator relative to all prior literature. The well-studied ingredients are local estimator combination, one-step distributed optimization, and split-sample bias correction \citep{zhang2013,huang2019,dhaene2015}. Its role here is to establish a self-contained rate calculation and a specific communication--computation comparison suitable for further research.296\end{remark}297298\section{Exact exponential-model analysis}299\subsection{An analytically tractable likelihood}300Suppose $X_i\stackrel{\mathrm{iid}}\sim\operatorname{Exp}(\theta_0)$ using the \emph{rate} parameterization. The average log-likelihood is301\begin{equation}302 L_n(\theta)=\log\theta-\theta\bar X_n,303 \qquad\widehat\theta_n=\bar X_n^{-1}.304\end{equation}305Let $S_n=\sum_{i=1}^n X_i$, so $S_n\sim\operatorname{Gamma}(n,\text{rate}=\theta_0)$. The Gamma inverse moments give, for $n>2$,306\begin{align}307 \EE[\widehat\theta_n]&=\theta_0\frac n{n-1},\label{eq:expmom1}\\308 \Var(\widehat\theta_n)&=\theta_0^2\frac{n^2}{(n-1)^2(n-2)}.\label{eq:expmom2}309\end{align}310Here $\psi(X)=\theta_0-\theta_0^2X$. The remainder has the \emph{exact} form311\begin{equation}312 R_n=\widehat\theta_n-\theta_0-\frac1n\sum_i\psi(X_i)313 =\frac{(1-\theta_0\bar X_n)^2}{\bar X_n},\label{eq:expremainder}314\end{equation}315whose moments are316\begin{align}317 \EE R_n&=\frac{\theta_0}{n-1}=\theta_0(n^{-1}+n^{-2}+O(n^{-3})),\\318 \EE R_n^2&=\theta_0^2\frac{3n+2}{n(n-1)(n-2)}=O(n^{-2}).\label{eq:expR2}319\end{align}320These identities verify the key strong mean-square and bias assumptions for this example, with $b_1=b_2=\theta_0$. The score and Hessian are $g_N(\theta)=1/\theta-\bar X_N$ and $H_N(\theta)=-1/\theta^2$; local Newton regularity holds around every fixed $\theta_0>0$.321322\subsection{The exact jackknife bias}323The local full and half estimators have means $\theta_0m/(m-1)$ and $\theta_0m/(m-2)$ respectively. Therefore324\begin{align}325 \EE J_k&=2\theta_0\frac m{m-1}-\theta_0\frac m{m-2},\nonumber\\326 \boxed{\;\EE J_k-\theta_0=-\frac{2\theta_0}{(m-1)(m-2)}.\;}\label{eq:exactbias2}327\end{align}328The identity confirms the $m^{-2}$ order without relying on an asymptotic expansion. The split estimates are dependent because they share observations with the full estimate; accounting for this dependence is essential when evaluating variance. Nevertheless, \cref{eq:expremainder} and the decomposition in \cref{eq:rj} ensure $\Var(J_k)=\theta_0^2/m+O(m^{-2})$ as $m\to\infty$.329330\subsection{Exact Newton identity and boundary distribution}331A Newton step started from $x$ simplifies to332\begin{equation}333 \mathcal N_N(x)=2x-x^2\bar X_N.334\end{equation}335Because $\thn=1/\bar X_N$, the following error formula holds \emph{exactly}:336\begin{equation}337 \boxed{\quad\mathcal N_N(x)-\thn338 =-\frac{(x-\thn)^2}{\thn}.\quad}\label{eq:exactnewton}339\end{equation}340341\begin{proposition}[Sharp boundary for the specified exponential construction]\label{prop:boundary}342Suppose $N\to\infty$, $m=N/K$ is even and $K/N^{7/8}\to c\in(0,\infty)$. For the half-block jackknife average and one Newton update,343\begin{equation}344 \sqrt N\bigl(\theta^{(1)}-\thn\bigr)345 \xrightarrow{p}-4\theta_0c^4,\label{eq:frontiergap}346\end{equation}347and hence348\begin{equation}349 \sqrt N(\theta^{(1)}-\theta_0)350 \Rightarrow\mathcal N(-4\theta_0c^4,\theta_0^2).\label{eq:frontierdist}351\end{equation}352\end{proposition}353\begin{proof}354From \cref{eq:exactbias2}, the jackknife average has bias asymptotic to $-2\theta_0m^{-2}$. Since $m=N/K\sim N^{1/8}/c$,355\begin{equation}356 N^{1/4}(\EE\theta^{(0)}-\theta_0)\to-2\theta_0c^2.357\end{equation}358The centered average fluctuates on scale $N^{-1/2}$, since its leading influence function averages all $N$ observations and the remaining centered term has standard deviation $O((m\sqrt K)^{-1})=o(N^{-1/4})$. Likewise $\thn-\theta_0=O_p(N^{-1/2})$. Therefore359\begin{equation}360 N^{1/4}(\theta^{(0)}-\thn)\xrightarrow{p}-2\theta_0c^2.361\end{equation}362Use the exact identity \cref{eq:exactnewton}, $\thn\to_p\theta_0$, and the continuous mapping theorem to obtain \cref{eq:frontiergap}. The centralized MLE satisfies $\sqrt N(\thn-\theta_0)\Rightarrow\mathcal N(0,\theta_0^2)$, proving \cref{eq:frontierdist}.\end{proof}363364\begin{remark}[Not a lower bound for all distributed estimators]365In the exponential model, each worker can transmit $\sum_{i\in\D_k}X_i$. Summing these $K$ scalars recovers $\thn=N/\sum_i X_i$ \emph{exactly}. Thus the sharp frontier in \cref{prop:boundary} concerns this particular jackknife-plus-one-step algorithm; it is not an information-theoretic limitation of distributed likelihood estimation.366\end{remark}367368\section{Higher-order local extrapolation}369\subsection{Cancelling two bias terms with four quarter-block fits}370Assume $m$ is divisible by four. Let $A_{k,2}$ denote the average of MLEs on the two disjoint halves and $A_{k,4}$ the corresponding average on four disjoint quarter-blocks. Define371\begin{equation}372 \boxed{\quad J_k^{[3]}=\frac83\widehat\theta_{k,m}-2A_{k,2}+\frac13 A_{k,4}.\quad}\label{eq:jack3}373\end{equation}374The coefficients solve375\begin{align}376 \tfrac83-2+\tfrac13&=1,\\377 \tfrac83-2(2)+\tfrac13(4)&=0,\\378 \tfrac83-2(2^2)+\tfrac13(4^2)&=0.379\end{align}380The first identity preserves the leading influence term; the others cancel $n^{-1}$ and $n^{-2}$ terms of a common bias expansion. With a third-order bias expansion $\EE R_n=b_1n^{-1}+b_2n^{-2}+b_3n^{-3}+O(n^{-4})$, the local bias is381\begin{equation}382 \EE[J_k^{[3]}]-\theta_0=8b_3m^{-3}+O(m^{-4}).\label{eq:bias3}383\end{equation}384The analog of the $L^2$ triangle inequality establishes a centered remainder of order $\op((m\sqrt K)^{-1})$ after averaging workers, despite sharing nested data. Consequently385\begin{equation}386 \theta^{(0,[3])}-\thn387 =\op\left(m^{-3}+\frac1{m\sqrt K}+N^{-1}\right).\label{eq:init3}388\end{equation}389390\begin{corollary}[A sufficient $11/12$ regime]\label{cor:eleven}391Suppose \cref{ass:lin,ass:hess} hold for $n\in\{m,m/2,m/4,N\}$, and the bias expansion through $n^{-3}$ holds with $O(n^{-4})$ remainder. For $m$ divisible by four and $K=o(N^{11/12})$, one global Newton step starting from $K^{-1}\sum_kJ_k^{[3]}$ is first-order equivalent to centralized MLE.392\end{corollary}393\begin{proof}394From \cref{eq:init3} and \cref{lem:newton},395\begin{equation*}396 \sqrt N\|\theta^{(1,[3])}-\thn\|397 =\op\left(\frac{K^6}{N^{11/2}}+\frac K{N^{3/2}}+N^{-3/2}\right)=o_p(1).398\end{equation*}399The first expression tends to zero if $K=o(N^{11/12})$; the others do so automatically.\end{proof}400401\subsection{Exact higher-order bias in the exponential model}402Using the exact mean in \cref{eq:expmom1}, for $m>4$,403\begin{equation}404 \boxed{\quad\EE[J_k^{[3]}]-\theta_0405 =\frac{8\theta_0}{(m-1)(m-2)(m-4)}.\quad}\label{eq:exactbias3}406\end{equation}407The bias magnitude is indeed $O(m^{-3})$. However, higher-order jackknife coefficients can amplify the finite-sample variance; the correction also requires local fits on quarter-blocks, which can have unstable or nonexistent moments for very small $m$. This is a reason to compare bias, variance, wall-clock time, and convergence failures rather than optimizing the asymptotic exponent alone.408409\subsection{A general bias--synchronization tradeoff}410Consider a fixed bias-removal order $r\ge1$. Suppose a local combination retains the same influence function and achieves411\begin{equation}412 \theta_r^{(0)}-\thn413 =\op\left(m^{-r}+\frac1{m\sqrt K}+N^{-1}\right),\label{eq:generalinit}414\end{equation}415with the requisite uniform moment bounds. After $q\ge0$ locally regular global Newton updates, repeated quadratic error contraction yields416\begin{equation}417 \|\theta_r^{(q)}-\thn\|418 =\op\left[419 \left(m^{-r}+\frac1{m\sqrt K}+N^{-1}\right)^{2^q}420 \right].\label{eq:qsteps}421\end{equation}422\begin{proposition}[Sufficient rate envelope]\label{prop:envelope}423For fixed integers $r\ge1$ and $q\ge0$, under \cref{eq:generalinit} and the required local quadratic contraction at each update, a sufficient condition for $\sqrt N(\theta_r^{(q)}-\thn)\to_p0$ is424\begin{equation}425 \boxed{\quad K=o\left(N^{\,1-1/(2^{q+1}r)}\right).\quad}\label{eq:envelope}426\end{equation}427\end{proposition}428\begin{proof}429With $m=N/K$, the bias component of the scaled remainder is $\sqrt N m^{-r2^q}=N^{1/2}(K/N)^{r2^q}$, which vanishes precisely under the stated sufficient regime. The centered-remainder component is $\sqrt N(\sqrt K/N)^{2^q}$. For $q=0$ it equals $\sqrt{K/N}$, which vanishes since the displayed exponent is less than one; for $q\ge1$ it is at most of order $N^{1/2-2^{q-1}}$ when $K\le N$, hence vanishes. The full-sample term $\sqrt N(N^{-1})^{2^q}$ also vanishes. The sum of finitely many terms may be bounded by the corresponding sum of powers, finishing the proof.\end{proof}430431\begin{table}[ht]432\centering\small433\caption{Sufficient partition-growth exponents under fixed-dimensional regularity. These are \emph{not} hardware speedups or universal lower bounds.}\label{tab:frontiers}434\begin{tabularx}{\linewidth}{@{}Xccc@{}}435\toprule436Initialization & Local bias order $r$ & Global Newton steps $q$ & $K=o(N^\alpha)$\\437\midrule438Raw local mean & 1 & 0 & $\alpha=1/2$\\439Raw local mean & 1 & 1 & $\alpha=3/4$\\440Raw local mean & 1 & 2 & $\alpha=7/8$\\441Half-block jackknife & 2 & 0 & $\alpha=3/4$\\442Half-block jackknife & 2 & 1 & $\alpha=7/8$\\443Quarter-block extrapolation & 3 & 0 & $\alpha=5/6$\\444Quarter-block extrapolation & 3 & 1 & $\alpha=11/12$\\445\bottomrule446\end{tabularx}447\end{table}448449\begin{remark}450\Cref{prop:envelope} is a rate-accounting result, conditional on an initializer satisfying \cref{eq:generalinit}; it is \emph{not} a proof that such an estimator exists for every $r$ or every econometric model. Constructing stable higher-order initializers in heterogeneous or high-dimensional settings remains an open implementation and theory problem.451\end{remark}452453\section{Numerical experiment: exact Gamma simulation}454\subsection{Design and reproducibility}455To isolate estimator behavior from optimization and simulation discretization, we draw the sufficient statistics of exponential observations exactly. For each worker $k$, we generate four independent Gamma variates, each representing the sum of $m/4$ independent exponential observations; their sum is the block sufficient statistic. From these sums we evaluate local full, half, and quarter MLEs; the centralized MLE follows from the sum over all workers.456457The common sample size is $N=2^{20}=1{,}048{,}576$, the true rate is $\theta_0=1$, and the partition sizes are $m\in\{256,64,16\}$, corresponding to $K\in\{4{,}096;16{,}384;65{,}536\}$. We run $B=500$ independent Monte Carlo replications in each design with NumPy's default random generator initialized at seed \texttt{20261028}. These are \emph{simulated data partitions, not physical computers}. We compare raw averaging with one or two Newton steps and both jackknife corrections with one Newton step. The estimand for each replication is the gap to the \emph{same-sample centralized MLE}, scaled by $\sqrt N$:458\begin{equation}459 \Delta_b=\sqrt N(\widehat\theta_b^{\mathrm{dist}}-\widehat\theta_{N,b}).\label{eq:mcgap}460\end{equation}461A gap to centralized MLE directly measures the numerical-statistical approximation, not the sampling error $\widehat\theta_{N,b}-\theta_0$. The exact code and the full 500-replication outputs are included in the accompanying reproducibility bundle.462463\subsection{Results}464\begin{table}[ht]465\centering\small466\caption{Mean scaled gap $\overline\Delta$ across 500 independent Monte Carlo replications. More negative numbers indicate a larger downward difference relative to centralized MLE. Rounding to six decimal places.}\label{tab:simulation}467\begin{tabular}{@{}lrrr@{}}468\toprule469Estimator & $m=256$ & $m=64$ & $m=16$\\470\midrule471Raw mean + 1 Newton & $-0.015728$ & $-0.258367$ & $-4.549613$\\472Raw mean + 2 Newton & $\approx 0$ & $-0.000065$ & $-0.020216$\\473Half-block jackknife + 1 Newton & $-0.000016$ & $-0.000331$ & $-0.093720$\\474Quarter-block jackknife + 1 Newton & $-0.000030$ & $-0.000147$ & $-0.011034$\\475\bottomrule476\end{tabular}477\end{table}478479\begin{figure}[ht]480\centering481\includegraphics[width=\linewidth]{figures/mc_comparison.pdf}482\caption{Absolute mean and root mean squared scaled gaps to the centralized MLE, calculated from the reproducibility script. Logarithmic axes make both small and large discrepancies visible. For the quarter-block method, its increased finite-sample variability can sometimes offset reduced bias.}\label{fig:mc}483\end{figure}484485The half-block correction substantially reduces the one-step gap for large local samples, consistent with \cref{thm:seven}. At $m=16$, higher-order correction improves the mean one-step gap still further, but also increases local fitting requirements and possible sensitivity to extreme quarter-sample MLEs. With two Newton steps, uncorrected averaging performs competitively: this is expected from \cref{cor:raw}, not evidence against the theory. The simulation is not a runtime evaluation and provides no evidence that one strategy is faster on actual distributed nodes.486487\subsection{What the simulation does not establish}488Three caveats are material. First, $N$ is fixed, whereas \cref{thm:seven,cor:eleven} describe triangular-array limits. Second, the Monte Carlo experiment uses exact aggregate statistics, hence omits network latency, local iterative optimizer cost, load imbalance, failure recovery, and floating-point differences across devices. Third, because the Gamma model has an elementary sufficient statistic, any empirical advantage here is a property of the estimator \emph{approximation}, not an argument against exact aggregation.489490\section{Communication and computation accounting}491\subsection{Payload and synchronization rounds}492For $p$ parameters and dense symmetric Hessians, each worker's initial jackknife estimate is a $p$-vector. One subsequent derivative evaluation requires one $p$-vector and $p(p+1)/2$ unique Hessian entries per worker. Excluding protocol headers, numerical precision and broadcasts, the total amount of worker-to-coordinator data is on the order of493\begin{equation}494 \boxed{\quad K\left[2p+\frac{p(p+1)}2\right]\quad\text{numbers}.}\label{eq:comm}495\end{equation}496The two communication \emph{phases} are not reducible to one if the second-phase score and Hessian are evaluated at the common, data-dependent $\theta^{(0)}$. One can transmit asynchronous or compressed derivatives, but then the current mathematical guarantee does not automatically apply. Dense-Hessian payload and coordinator matrix inversion may become prohibitive when $p$ is large.497498\begin{table}[ht]499\centering\small500\caption{Logical comparison. A ``derivative phase'' is one globally synchronized derivative evaluation after the initial communication of local estimates. Local-fit counts assume direct refitting on nested partitions.}\label{tab:comm}501\begin{tabularx}{\linewidth}{@{}Xccc@{}}502\toprule503Method & Local MLE fits / worker & Derivative phases & Sufficient exponent\\504\midrule505Raw mean + 1 Newton & 1 & 1 & $3/4$\\506Raw mean + 2 Newton & 1 & 2 & $7/8$\\507Half-block jackknife + 1 Newton & 3 & 1 & $7/8$\\508Quarter-block extrapolation + 1 Newton & 7 & 1 & $11/12$\\509\bottomrule510\end{tabularx}511\end{table}512513The meaningful systems-level claim is conditional: if an additional full-cluster synchronization has greater cost than two extra local likelihood fits, half-block jackknifing may be favorable; if local fits dominate and network latency is low, extra Newton steps may be faster. The cut-off depends on model dimension, optimizer complexity, scheduling, bandwidth, and numerical conditioning.514515\subsection{A cost model for an empirical follow-up}516Let $T_{\mathrm{fit}}(m,p)$ denote the elapsed time for one local MLE fit, $T_{\mathrm{der}}(m,p)$ for one local derivative evaluation, $L$ the per-phase synchronization latency, and $B$ the effective worker-to-coordinator bandwidth. A first-order approximation is517\begin{align}518 T_{\mathrm{half+N}}&\approx 3T_{\mathrm{fit}}(m,p)+T_{\mathrm{der}}(m,p)+2L+\frac{bK[2p+p(p+1)/2]}{B},\\519 T_{\mathrm{raw+2N}}&\approx T_{\mathrm{fit}}(m,p)+2T_{\mathrm{der}}(m,p)+3L+\frac{bK[3p+p(p+1)]}{B},520\end{align}521where $b$ is bytes per transmitted scalar and the expression deliberately suppresses overlapping parallelism, coordinator inversion, broadcasts, network contention, and data loading. These are accounting illustrations, \emph{not} validated time-complexity theorems. An actual benchmark must record energy, cold starts, memory use, fit failures, $p$, hardware topology, and communication bytes independently of statistical accuracy.522523\section{Limitations and research agenda}524\subsection{Identification, misspecification, and dependence}525The Fisher-information variance in the main theorem presumes a correctly specified regular model. Under misspecification, the asymptotic covariance is generally the sandwich form and the common influence function must be modified. If observations exhibit cross-node dependence or unequal distributions, the independent-block variance calculation in \cref{prop:init} changes. Fixed effects with incidental parameters can have nonstandard bias expansions; indeed, that setting is one of the motivations for \citet{dhaene2015}. Neither validity nor exact exponents should be assumed to transfer without new proofs.526527\subsection{High-dimensional Hessians and nonconvex likelihoods}528Our proof assumes fixed $p$. When $p$ grows with $N$, the norm of the Fisher inverse and the constants in the high-moment expansion can depend on $p$, and the $K$ exponent may weaken. In nonconvex likelihoods, local workers may converge to different stationary points, while a full Newton step can escape the basin in which the quadratic bound holds. Regularization, trust regions and generalized linear-system solvers may mitigate these issues but require separate theory. Stronger generality without control of such effects would overstate the current result.529530\subsection{Statistical versus computational optimality}531We have not proved that the number of communications can be independent of $N$ in a minimax sense, nor that this particular scheme achieves a lower bound over all algorithms. Indeed, exact sufficient-statistic aggregation defeats any supposed general lower bound in the exponential model. The appropriate open question is narrower: \emph{over a specified class of nonlinear likelihoods, and under a model of allowable local computation and communication, what Pareto frontier links local refits, derivative synchronizations, transferred bytes, and statistical precision?}532533A useful research program would (i) characterize the finite-sample risk of high-order jackknife starts; (ii) derive feasible model-specific analytical bias corrections that use fewer local optimizations; (iii) compare with modern multi-round techniques such as CEASE \citep{fan2023}; (iv) extend the expansions to clustered or time-dependent panels; and (v) implement a wall-clock benchmark across physical nodes, varying the cost ratio of fitting to synchronization.534535\section{Conclusion}536This note develops and verifies a local-bias-correction-plus-global-Newton construction for distributed maximum likelihood estimation. Under explicit uniform mean-square remainders, second-order bias expansions, and local Hessian regularity, a half-block jackknife initializer followed by one global Newton step is first-order equivalent to centralized MLE when $K=o(N^{7/8})$. The stochastic component of the jackknife remainder has been kept throughout the proof. A third-order local extrapolation yields the sufficient exponent $11/12$ under stronger assumptions. Exact Gamma calculations establish the rate and a nonzero boundary shift in the exponential model, while a reproducible Monte Carlo experiment confirms the expected qualitative pattern.537538The scientific interpretation is appropriately limited. Related distributed one-step estimation and higher-order jackknife ideas are established, and the method has not yet been benchmarked in a real networked environment. The potentially productive direction is the \emph{quantitative exchange between local correction and global synchronization}, rather than an unsupported claim to have surpassed all existing maximum likelihood algorithms.539540\section*{Statements and reproducibility}541\addcontentsline{toc}{section}{Statements and reproducibility}542\textbf{Code and data.} The accompanying package contains \texttt{monte\_carlo.py}, the simulation CSV outputs, and figure sources. The simulation requires Python, NumPy, and Matplotlib. All generated observations are synthetic; no proprietary data or user data were used.\\543\textbf{Authorship.} The designation ``GPT-6.1'' is included at the request of the user preparing this draft. It refers to an AI-assistance authorship label, rather than a verified human contributor. Responsibility for scientific validation and compliance with any target journal's authorship policy rests with the submitting researcher.\\544\textbf{Research status.} Unsubmitted methodological working paper. Theoretical claims are explicitly assumption-dependent. Priority and novelty have not been exhaustively established.545546\begin{thebibliography}{99}547\bibitem[Zhang et~al.(2013)Zhang, Duchi, and Wainwright]{zhang2013}548Zhang, Y., Duchi, J.~C., and Wainwright, M.~J. (2013).549Communication-efficient algorithms for statistical optimization.550\emph{Journal of Machine Learning Research}, 14:3321--3363.551\url{https://www.jmlr.org/papers/v14/zhang13b.html}.552553\bibitem[Huang and Huo(2019)]{huang2019}554Huang, C. and Huo, X. (2019).555A distributed one-step estimator.556\emph{Mathematical Programming}, 174(1--2):41--76.557\url{https://doi.org/10.1007/s10107-019-01369-0}.558559\bibitem[Dhaene and Jochmans(2015)]{dhaene2015}560Dhaene, G. and Jochmans, K. (2015).561Split-panel jackknife estimation of fixed-effect models.562\emph{The Review of Economic Studies}, 82(3):991--1030.563\url{https://doi.org/10.1093/restud/rdv007}.564565\bibitem[Fan et~al.(2023)Fan, Guo, and Wang]{fan2023}566Fan, J., Guo, Y., and Wang, K. (2023).567Communication-efficient accurate statistical estimation.568\emph{Journal of the American Statistical Association}, 118(542):1000--1010.569\url{https://doi.org/10.1080/01621459.2021.1969238}.570571\bibitem[Efron(1982)]{efron1982}572Efron, B. (1982). \emph{The Jackknife, the Bootstrap and Other Resampling Plans}.573SIAM, Philadelphia.574\url{https://doi.org/10.1137/1.9781611970319}.575\end{thebibliography}576577\appendix578\section{Proof details for the exponential remainder}579Let $Y=\theta_0\bar X_n$, so $nY\sim\operatorname{Gamma}(n,\text{rate}=1)$. We have $\widehat\theta_n=\theta_0/Y$ and $R_n=\theta_0(1-Y)^2/Y=\theta_0(Y^{-1}-2+Y)$. The requisite moments for $n>2$ are580\begin{equation}581 \EE Y=1,\quad\EE Y^2=1+1/n,\quad582 \EE Y^{-1}=\frac n{n-1},\quad583 \EE Y^{-2}=\frac{n^2}{(n-1)(n-2)}.584\end{equation}585Therefore586\begin{align*}587 \EE R_n&=\theta_0\left(\frac n{n-1}-2+1\right)=\frac{\theta_0}{n-1},\\588 \EE R_n^2&=\theta_0^2\EE(Y^{-2}-4Y^{-1}+6-4Y+Y^2)\\589 &=\theta_0^2\left[\frac{n^2}{(n-1)(n-2)}-\frac{4n}{n-1}+3+\frac1n\right]\\590 &=\theta_0^2\frac{3n+2}{n(n-1)(n-2)}.591\end{align*}592This supplies an elementary verification of \cref{eq:expR2}. For nested subsamples, the same $\psi$ is used for each sample, which is why the linear term cancels in the jackknife identity.593594\section{Exact bias of the quarter-block extrapolation}595At full, half and quarter-block sizes, the rate MLE has expectations $\theta_0m/(m-1)$, $\theta_0m/(m-2)$ and $\theta_0m/(m-4)$ respectively. Inserting these values into \cref{eq:jack3} yields596\begin{align*}597 \EE J_k^{[3]}-\theta_0598 &=\theta_0\left[\frac{8/3}{m-1}-\frac4{m-2}+\frac{4/3}{m-4}\right]\\599 &=\frac{8\theta_0}{(m-1)(m-2)(m-4)},600\end{align*}601proving \cref{eq:exactbias3}. The high-order bias coefficients are positive in this particular model, whereas the half-block bias in \cref{eq:exactbias2} is negative. Such alternating signs are one reason direct extrapolation is not necessarily monotone in finite-sample accuracy.602603\end{document}604