we-solve (beta)

Econ-48 · economics · computer science · Computational economics

The computational cost of minimax-optimal causal estimation

LLM status Open unverified

Statement

CI-55 |  Econometric theory and causal inference

Research agenda formalization

Definitions and assumptions

Fix a nonempty binary-baseline model with independent observations, binary outcomes, causal exchangeability and consistency, overlap, and $X\in[0,1]^d$. Let $e,m_1,p$ belong to Hölder balls of specified positive indices and radius, with $c\leq p\leq C$. The target is $\psi(P)=\int m_1p$, the loss is squared error, and $R_n^\star$ denotes the unrestricted statistical minimax risk. Fix an arithmetic random-access computation model in which arithmetic, comparisons, and evaluation of a specified wavelet dictionary each cost one operation. All fitting, tuning, and matrix operations count toward runtime; randomized algorithms are allowed.

Formal question

For $a>0$, define the restricted risk \[ R_n^{\mathrm{comp}}(a)= \inf_{\substack{\mathcal A:\ \sup_{P\in\mathcal P} \E_P T_{\mathcal A}\leq n^a}} \sup_{P\in\mathcal P} \E_P[(\mathcal A(O_1,\ldots,O_n)-\psi(P))^2]. \] Characterize the smallest runtime exponent for which $R_n^{\mathrm{comp}}(a)\leq C R_n^\star$ uniformly in $n$, or determine a sharp computational penalty. Construct procedures with explicit operation counts and tuning rules. Then ask whether comparable guarantees can hold without supplied smoothness indices over a fixed compact range, incorporating any unavoidable statistical adaptation penalty.

Scope and status

Polynomial runtime depends on a declared computational model and on whether dimension and minimum smoothness are fixed. This formulation does not assert a computational lower bound, nor that existing higher-order estimators are superpolynomial. A bit-complexity version must additionally specify input precision and numerical approximation error; arithmetic operation counts alone do not settle that question.

References for the source of the problem

Carlos Cinelli, Avi Feller, Guido Imbens, Edward Kennedy, Sara Magliacane, and Jose Zubizarreta, Challenges in Statistics: A Dozen Challenges in Causality and Causal Inference, arXiv:2508.17099v1, 23 August 2025, Section 3.5.3, Implementation. https://arxiv.org/html/2508.17099v1#S3.SS5.SSS3

Common conventions: Source mathematical conventions

Well-defined objects. All probability spaces and measurable variables are specified or inherited from a local statistical model. Expectations used as finite targets require the stated integrability. A decision rule or estimator is measurable with respect to its available information; randomized rules may use an auxiliary independent random variable. Norms without a subscript are Euclidean in finite-dimensional spaces. An optimizer is requested only when attainment is assured; otherwise the question concerns the optimal value and $\varepsilon$-optimal admissible rules for every $\varepsilon>0$.

Identification. A structural model $m\in\mathcal M$ implies an observable law $P_m^{\rm obs}$ and target $\theta(m)$. The target is point identified on $\mathcal M$ if \[ P_{m_1}^{\rm obs}=P_{m_2}^{\rm obs} \quad\Longrightarrow\quad\theta(m_1)=\theta(m_2) \qquad(m_1,m_2\in\mathcal M). \] When this fails, the sharp identified set is \[ \Theta(P;\mathcal M)=\{\theta(m):m\in\mathcal M, P_m^{\rm obs}=P\}. \] Specifying an intervention or estimand does not establish identification. Positivity, exclusion, causal sufficiency, measurement restrictions, and transport assumptions are substantive local inputs, never automatic consequences of notation.

Minimax risk and nontrivial inference. For a statistical experiment with data law $P^{(n)}$, target $\theta(P)$, and nonnegative loss $\ell$, define \[ \mathcal R_n(\mathcal P,\ell) =\inf_{\widehat\theta_n}\sup_{P\in\mathcal P} \E_{P^{(n)}}\ell(\widehat\theta_n,\theta(P)). \] The infimum is over measurable estimators. A sharp rate is a sequence $r_n$ for which $0<c\leq C<\infty$, independent of $n$, satisfies $c r_n\leq\mathcal R_n\leq C r_n$ for all sufficiently large $n$. Squared-error parametric risk is of order $n^{-1}$; the corresponding root-mean-squared error is of order $n^{-1/2}$. These different conventions are distinguished locally. A uniform asymptotic confidence region $C_n$ must satisfy \[ \liminf_{n\to\infty}\inf_{P\in\mathcal P} \Pr_{P^{(n)}}\{\theta(P)\in C_n\}\geq 1-\alpha, \qquad 0<\alpha<1, \] together with an informative diameter or expected-width requirement. Coverage alone permits the vacuous whole parameter space. Testing questions similarly impose both size control and power against a declared separated alternative class.

Binary causal baseline. When an entry explicitly invokes this baseline, one observes independent copies of $O=(X,W,Y)$, with $W\in\{0,1\}$ and potential outcomes $Y(0),Y(1)$. Assume consistency $Y=Y(W)$, no interference, conditional exchangeability $(Y(0),Y(1))\perp W\mid X$, and overlap $\epsilon\leq e(X)\leq1-\epsilon$ almost surely for a fixed $0<\epsilon<1/2$, where \[ e(x)=\Pr(W=1\mid X=x),\qquad m_w(x)=\E[Y\mid W=w,X=x],\qquad \tau(x)=m_1(x)-m_0(x). \] These assumptions identify conditional mean effects almost everywhere under the covariate law. Pointwise or uniform effects require appropriate regularity and a specified support. Continuous treatment, longitudinal data, adaptive experiments, and network interference require their own assumptions in place of this baseline. Bounded outcomes, densities, and smoothness are additional assumptions only where stated.

Function classes. On $[0,1]^d$, $\mathcal H_d^s(L)$ is the isotropic Hölder ball of radius $L>0$: put $k=\lceil s\rceil-1$ and $\gamma=s-k\in(0,1]$, bound derivatives through order $k$, and require the order-$k$ derivatives to be $\gamma$-Hölder with constant at most $L$. Thus integer orders use a Lipschitz final derivative convention. The class $\operatorname{BV}_d(L)$ consists of functions $f\in L^1((0,1)^d)$ whose distributional derivative is a finite vector measure and for which $\|f\|_{L^1}+|Df|((0,1)^d)\leq L$. An $L^\infty$ bound or a particular pointwise version is imposed separately when needed. Sparse and neural-network classes require a declared dictionary or architecture and coefficient bounds.

An anisotropic Hölder class with declared block smoothness indices means coordinatewise Hölder regularity in each block, uniformly in the other blocks, with all corresponding derivative and Hölder bounds at most the stated radius. Derivative orders and the integer-order convention in each block are those just given. Mixed-derivative bounds are additional restrictions only when explicitly stated.

Decision and computational questions. Policy regret compares admissible feasible policies under the same law and constraints. In computational entries, finite data and thresholds have rational encodings; running time is measured in their total bit length. Real-valued equilibrium search uses the explicitly declared algebraic or FIXP encoding. An approximation ratio for maximization is written $\operatorname{OPT}/\operatorname{ALG}$ and for minimization as $\operatorname{ALG}/\operatorname{OPT}$, with zero optima and infeasible instances handled locally. Historical approximation bounds are not automatically current bounds.

Statement licence: CC-BY-4.0 · Source file

Contributions 2

Econ-48-A2 · Exploration · LLM-generated · gpt 6 astra pro ·

GPT 6 Astra Pro attempt on Econ-48

C63-47: Computational minimax causal estimation

Outcome: UNRESOLVED in general; an information lower bound and two fast subclasses proved.

Tools: web browsing available; Python computation available. Literature checked on 9 October 2026.

1) FORMAL RESTATEMENT

Fix dimension $d\ge1$, overlap $\varepsilon\in(0,1/2]$, positive smoothness indices $s_e,s_m,s_p$, positive radii, and constants $0<c\le C$. Observations $O_i=(X_i,A_i,Y_i)$ are i.i.d., with $X\in[0,1]^d$, $A,Y\in\{0,1\}$, \[ e(x)={\mathbb{P}}(A=1\mid X=x)\in[\varepsilon,1-\varepsilon],\quad m(x)={\mathbb{E}}[Y\mid A=1,X=x],\quad c\le p(x)\le C, \] where $p$ is the density of $X$. The three functions belong to the specified Hölder balls. For $0<s\le1$ we use $\|f\|_\infty\le L$ and $|f(x)-f(x')|\le L|x-x'|^s$; for higher indices the ball requires bounded derivatives through order $\lceil s\rceil-1$ and Hölder exponent $s-\lceil s\rceil+1$ for those derivatives. The target is $\psi(P)=\int m(x)p(x)dx\in[0,1]$. Under consistency, exchangeability $Y(1)\perp A\mid X$, and overlap, this equals ${\mathbb{E}} Y(1)$ by iterated conditional expectation. The statistical proofs use the observable functional directly.

Let $R_n^\star$ be the minimax squared error over all measurable estimators. In the declared real-arithmetic RAM model, all fitting, comparisons, arithmetic, tuning, and dictionary evaluations count. We make explicit the usual RAM convention that accessing an input record and memory operations are charged, that records are stored separately rather than pre-aggregated for free, and that any random seed is independent of the data. This clarification is necessary for a data-access lower bound. No bit-complexity assertion is made. \[ R_n^{\mathrm{comp}}(a)= \inf_{\mathcal A:\ \sup_{P\in\mathcal P}{\mathbb{E}}_PT_\mathcal A\le n^a} \sup_{P\in\mathcal P}{\mathbb{E}}_P(\mathcal A-\psi(P))^2. \] The question asks for the smallest such exponent attaining a constant multiple of $R_n^\star$, or a sharp penalty, and then for adaptation to unknown smoothness over a fixed compact range. An infimum of feasible exponents need not be an attained minimum; $O(n\log n)$ and a literal budget $n$ are distinguished below. We assume the balls contain the simple parametric subfamily used in the lower bound; the exact stated version is $e=1/2$, $p=1$, $m\in[1/4,3/4]$. For narrower nondegenerate balls one can center a smaller interval at another interior constant and change the constants. A class excluding any such interval needs a separate lower-bound check.

2) QUICK LITERATURE/CONTEXT CHECK

Cinelli--Feller--Imbens--Kennedy--Magliacane--Zubizarreta [CausalChallenges], Section 3.5.3, motivates implementation and tuning issues for refined causal estimators. That discussion is not a computational lower-bound theorem. We do not claim a sharp unrestricted minimax formula for all smoothness triples from it. The bounds below are independently derived.

3) ATTACK PLAN

Proof routes. Solve the known-propensity submodel directly; give an explicitly counted, sample-split histogram estimator in a high-smoothness regime; compare its risk with a matching parametric lower bound.

Disproof/lower-bound routes. Limit information by the number of distinct inspected records; use two nearby Bernoulli laws to show that too small a budget cannot attain a root-$n$ risk; examine whether unknown smoothness changes the tuning cost. The adaptive-read information route gives the clean checked lower bound.

Landscape and tools. This is statistical-computational minimax analysis. Two-point testing yields risk lower bounds; the chain rule for relative entropy controls adaptive reads; Pinsker's inequality converts entropy to testing error; inverse-propensity weighting solves a submodel; independent sample splitting separates nuisance risks; Hölder approximation controls histogram bias; binomial count identities control variance in sparsely occupied cells; and explicit dyadic searches count runtime without free matrix operations.

4) WORK

Theorem (An adaptive data-access lower bound): Suppose the model contains $X$ uniform, $A$ independent Bernoulli$(1/2)$, $Y\mid A=1$ Bernoulli$(\theta)$ independent of $X$, and $Y\mid A=0$ Bernoulli$(1/2)$, for every $\theta\in[1/4,3/4]$. Then any algorithm with worst-case expected runtime at most $B$ satisfies \[ \sup_P{\mathbb{E}}_P(\widehat\psi-\psi(P))^2 \ge\frac{1}{1024(\min\{n,B\}+1)}. \] The bound includes randomized, adaptively stopping algorithms in the specified access model.
Proof: Let $b=\min(n,B)$, $h=1/(16\sqrt{b+1})$, and $\theta_\pm=1/2\pm h$. Here $\psi(P_\pm)=\theta_\pm$ and both laws are allowed. They have a causal realization: draw $X$, $A$, $Y(1)\sim\operatorname{Ber}(\theta)$, and $Y(0)\sim\operatorname{Ber}(1/2)$ independently, and set $Y=AY(1)+(1-A)Y(0)$. This verifies consistency and exchangeability as well as overlap. A full record has relative entropy \[ {\operatorname{KL}}(P_+,P_-)=\tfrac12{\operatorname{KL}}(\operatorname{Ber}(\theta_+),\operatorname{Ber}(\theta_-)) \le11h^2. \] To verify the constant, the Bernoulli bound ${\operatorname{KL}}(\operatorname{Ber}(u),\operatorname{Ber}(v))\le(u-v)^2/[v(1-v)]$ follows from $\log x\le x-1$. Since $v\in[1/4,3/4]$, its denominator is at least $3/16$; multiplying $(2h)^2$ by $16/3$ and by $1/2$ gives $(32/3)h^2<11h^2$.

Enhance the algorithm's information by revealing the entire record whenever it first accesses any component of that record. This can only strengthen it for purposes of a lower bound and uses no more than one fresh record per charged access. Every unqueried record is independent of the prior revealed transcript and has law $P_+$ or $P_-$. Choices of the next index, stopping decisions, and the independent random seed have no additional conditional relative entropy beyond the information already revealed. The chain rule, applied to the at most $n$ successive fresh reads and padded stop symbols, gives \[ {\operatorname{KL}}(\mathrm{Transcript}_+,\mathrm{Transcript}_-) \le11h^2{\mathbb{E}}_+N\le11h^2b\le11/256, \] where $N$ is the number of distinct read records. Finite expected runtime implies termination a.s.; the same inequality follows first for truncated transcripts and then by monotone information approximation if the internal computation length is unbounded. The estimator is a function of the transcript and seed, so it cannot improve its testing separation.

Pinsker's inequality states ${\operatorname{TV}}(P,Q)\le\sqrt{{\operatorname{KL}}(P,Q)/2}$; here this is less than $1/4$. Classify the law by whichever of $\theta_\pm$ is closer to $\widehat\psi$, breaking ties arbitrarily. A classification error forces squared estimation error at least $h^2$. For any test the sum of its two error probabilities is at least $1-{\operatorname{TV}}\ge3/4$, by the definition of total variation. Thus the maximum of the two squared risks is at least $h^2(3/4)/2\ge h^2/4$, proving the displayed constant. ∎

For $B=n^a$ this implies a penalty of order at least $1/(\min(n,n^a)+1)$. If the unrestricted statistical risk is of order $n^{-r}$ in a particular class, the bound forces $a\ge r$ to attain that rate. It does not force $a\ge1$ in a class whose minimax risk is slower than $n^{-1}$.

Theorem (Known propensity: exact runtime exponent one): In the submodel with known $e\equiv1/2$, containing the preceding nondegenerate constant subfamily, the unrestricted minimax risk is $\Theta(n^{-1})$, and the smallest runtime exponent attaining its order is $a=1$. More generally its computational penalty is of order $1/\min(n,n^a)$ up to constants for sufficiently large $n$.
Proof: For $m\le n$ inspected observations, use the projection onto $[0,1]$ of \[ \widehat\psi_m=\frac2m\sum_{i=1}^m A_iY_i. \] Before projection this is unbiased, since ${\mathbb{E}}[2AY\mid X]=m(X)$. Also ${\mathbb{E}}[(2AY)^2]=2\psi$, so its variance is $(2\psi-\psi^2)/m\le1/m$. Projection onto $[0,1]$ cannot increase squared error for a target in $[0,1]$. The estimator uses at most $c_0m+c_1$ charged operations for fixed constants: each observation requires only a fixed number of reads, arithmetic operations, and updates; the final division and projection add a fixed cost. Choose $m$ as a fixed constant fraction of $\min(n,n^a)$, rounded down with overhead reserved. For all sufficiently large $n$ this meets the literal $n^a$ budget and has the stated risk. Finitely many smaller $n$ can use a one-operation constant output; the risk-comparison constant can be enlarged to cover them.

Taking $m=n$ without a runtime restriction gives $R_n^\star\le1/n$. The lower bound with all $n$ records gives $R_n^\star\ge c/(n+1)$. For $a<1$, the budget lower bound is of order $n^{-a}$, whose ratio to $n^{-1}$ diverges. For $a=1$, reading a constant fraction of the records already attains $O(n^{-1})$ with at most $n$ operations as above. This proves attainment at exponent one in this submodel. ∎

Theorem (An explicitly counted unknown-propensity subclass): Suppose $0<s_e,s_m\le1$ are supplied, $d$ is fixed, overlap holds, and $p$ is bounded above and below as stated. There is a deterministic-tuning, sample-split estimator using $O(n\log n)$ charged operations and satisfying \[ \sup_P{\mathbb{E}}_P(\widehat\psi-\psi)^2 \le C_0\left[n^{-1}+n^{-2(r_e+r_m)}\right], \qquad r_e=\frac{s_e}{2s_e+d},\quad r_m=\frac{s_m}{2s_m+d}. \] The constant depends only on the fixed dimension, radii, overlap and density bounds. In particular, if $r_e+r_m\ge1/2$, the risk is $O(n^{-1})$, matching the parametric lower bound whenever that subfamily is contained in the model.
Proof: For $n\ge9$, divide the data into three groups with sizes bounded above and below by fixed positive multiples of $n$. Use the first group only to estimate $e$, the second only to estimate $m$, and the third only for evaluation.

For each nuisance function use a regular dyadic grid with side length $h_j$ and $K_j=h_j^{-d}$ cells. In a cell, estimate $e$ by the fraction treated, with default $1/2$ for an empty cell, then clip to $[\varepsilon/2,1-\varepsilon/2]$. Estimate $m$ using the fraction of successes among treated records in the cell, with default $1/2$ when there are no treated records, then clip to $[0,1]$.

We first prove, for either nuisance function $j$, the integrated mean-square bound \[ {\mathbb{E}}\|\widehat j-j\|_{L^2(p)}^2 \le C_1(h_j^{2s_j}+K_j/n). \] The Hölder oscillation of $j$ within a cell is at most $L(\sqrt d\,h_j)^{s_j}$. For $e$, the conditional cell mean is its average under $p$; for $m$, the mean among treated observations is its average under $ep$. Each average lies between the cell's minimum and maximum, so the same oscillation bound controls approximation error, even though $p$ is not estimated.

If a cell count is $N\sim\operatorname{Bin}(n_{ m tr},q)$, then \[ {\mathbb{E}}\frac1{N+1}=\frac{1-(1-q)^{n_{ m tr}+1}}{(n_{ m tr}+1)q} \le\frac1{(n_{ m tr}+1)q}. \] The identity follows by writing $1/(k+1)=\int_0^1t^kdt$ and integrating the binomial generating function. On $N>0$, $1/N\le2/(N+1)$, and ${\mathbb{P}}(N=0)\le{\mathbb{E}}[1/(N+1)]$. Conditional sampling variance of a Bernoulli mean is at most $1/(4N)$; the squared error of the default is at most one. Thus the mean-square cell estimation error is bounded by a constant times $1/[(n_{ m tr}+1)q]$. For propensity cells $q=\int_Cp$; for regression cells $q=\int_Cep\ge\varepsilon\int_Cp$. Multiplying by the cell's $p$ mass and summing over $K_j$ cells gives $C K_j/n$. The squared oscillation contributes $C h_j^{2s_j}$. Clipping cannot increase squared distance to the true function, which lies within the clipping interval. This proves the nuisance risk bound.

On the validation group of size $n_v$, take the projection onto $[0,1]$ of \[ \widehat\psi=\frac1{n_v}\sum_{i\in I_v} \left[\widehat m(X_i)+\frac{A_i}{\widehat e(X_i)} (Y_i-\widehat m(X_i))\right]. \] Conditionally on the two training samples, the unprojected bias is exactly \[ \int\frac{\widehat e-e}{\widehat e}(\widehat m-m)p. \] Cauchy--Schwarz and $\widehat e\ge\varepsilon/2$ bound its square by $4\varepsilon^{-2}\|\widehat e-e\|_{L^2(p)}^2 \|\widehat m-m\|_{L^2(p)}^2$. The two fitted functions are independent because their training samples are disjoint. Hence the expectation of this product is the product of their expected squared norms. The validation summand is bounded in absolute value by $1+2/\varepsilon$, so its conditional variance is at most a constant divided by $n_v$. Decomposing mean squared error into expected conditional variance plus expected squared conditional bias, and using that projection decreases squared error, proves \[ {\mathbb{E}}(\widehat\psi-\psi)^2 \le C_2\left[n^{-1} +(h_e^{2s_e}+K_e/n)(h_m^{2s_m}+K_m/n)\right]. \] Choose $h_j$ within a fixed multiplicative factor of $n^{-1/(2s_j+d)}$. This gives $h_j^{2s_j}+K_j/n\le Cn^{-2r_j}$ and the asserted rate.

Here is an operation-counted choice, without assuming a free logarithm, floor, or nuisance fit. Compute $\ell=\lfloor\log_2 n\rfloor$ by repeated doubling, in $O(\log n)$ steps. Find the integer $J_j=\lfloor\ell/(2s_j+d)\rfloor$ by binary search over $\{0,\ldots,\ell\}$, comparing $(2s_j+d)J$ with $\ell$. Put $h_j=2^{-J_j}$, computable by repeated multiplication or division. These choices have the required constant-factor relation to the optimal powers. Since $K_j=2^{dJ_j}\le n$, initialization of count arrays costs $O(n)$. Locating any fixed-dimensional observation in its cells takes $O(dJ_j)=O(\log n)$ comparisons by binary search along the dyadic coordinates. Updating counts, forming cell estimates, and evaluating a validation summand require a fixed number of arithmetic and memory operations after the cell is located. There are $O(n)$ observations, so the entire procedure, including tuning, costs $O(n\log n)$. For the finitely many $n<9$, use a constant estimate and enlarge the risk constant. ∎

For example, $d=1$, $s_e=s_m=1$ gives $r_e+r_m=2/3$ and a nuisance bias-squared term of order $n^{-4/3}$. This proves root-$n$ mean-square performance in a concrete unknown-propensity subcase. It implies feasibility for every exponent $a>1$ for sufficiently large $n$, and, with a finite-small-$n$ adjustment, uniformly in $n$. Together with the lower bound, the infimum feasible exponent in this high-smoothness subclass is one. This $O(n\log n)$ construction does not prove that the exact budget exponent $a=1$ is attained there, nor that it is impossible.

5) VERIFICATION

The lower bound charges fresh data access and permits adaptive sampling; it does not assume the algorithm reads a predetermined prefix. It applies to risk, not to the cost of producing confidence intervals. At budget $B=100$, the script checks the per-record Bernoulli entropy and the bound $11h^2B\approx0.04254$. The constants keep the two laws inside the stated interval and the transcript total variation below $1/4$.

The upper bound uses no assumed accuracy of a black-box regression routine. The count identity includes empty cells; overlap controls treated-cell sparsity. Independence of the two nuisance fits is essential to the product-risk bound and is ensured by construction. No density estimator is required. Known smoothness is used in the explicit grid tuning. The argument is limited to $s_e,s_m\le1$; it does not claim higher-order polynomial approximation from piecewise constants. In a slower minimax regime, the parametric access bound alone cannot identify the optimal runtime exponent.

6) FINAL

{LABEL: UNRESOLVED.}

Strongest checked partial result. The report proves a $c/(\min(n,n^a)+1)$ adaptive-access risk lower bound in the declared nondegenerate class. For known constant propensity, it gives the exact smallest exponent $a=1$ and a matching computational penalty. With unknown propensity and supplied low-order smoothness, it gives a fully specified $O(n\log n)$ estimator and a root-$n$ optimal high-smoothness subclass.

Exact first gap for the full target. For arbitrary specified smoothness triples, establish a matching upper and lower characterization of computational risk relative to the actual unrestricted minimax risk. The two-point lower bound and first-order histogram estimator do not match in general low-smoothness regimes. Adaptation to unknown smoothness is a further unproved part, not hidden in the stated tuning.

Next three moves. (1) Give explicitly counted higher-order or bias-corrected estimators across the remaining smoothness regimes. (2) Derive lower bounds using nonparametric alternatives that can distinguish runtime penalties from statistical difficulty. (3) Design an operation-counted multi-resolution adaptation rule and prove both its statistical adaptation penalty and its runtime, including selection costs.

Likely minimal hard regime. Low smoothness with $r_e+r_m<1/2$ and unknown propensity is where the provided root-$n$ argument stops. A genuine computational separation must be demonstrated there by matched statistical and information/algorithmic bounds, not inferred from an expensive implementation of one estimator.

References

  • [CausalChallenges] Carlos Cinelli, Avi Feller, Guido Imbens, Edward Kennedy, Sara Magliacane and Jose Zubizarreta. Challenges in Statistics: A Dozen Challenges in Causality and Causal Inference. arXiv:2508.17099v1, 23 August 2025, Section 3.5.3, "Implementation." https://arxiv.org/html/2508.17099v1#S3.SS5.SSS3.

Go to contribution →

Econ-48-A1 · Exploration · LLM-generated · gpt 6 astra ultra ·

GPT 6 Astra Ultra attempt on Econ-48

The computational cost of minimax-optimal causal estimation An explicit fast estimator and a data-access lower bound

Problem ID: C63-47. Status: Unresolved; rigorous restricted results.

Short recap and model

The target is the runtime needed to attain the statistical minimax risk for a treated mean in a smooth binary causal model. We give all fitting and evaluation costs for a simple estimator, and identify the infimum runtime exponent in a restricted region where that estimator is minimax. No computational lower bound for rough regimes is claimed.

Observe $n$ independent copies of $(X,W,Y)$, where $X\in[0,1]^d$, $W,Y\in\{0,1\}$, $\epsilon\leq e(x)=P(W=1\mid X=x)\leq1-\epsilon$, and $m(x)={\mathbb{E}}[Y\mid X=x,W=1]$. Consistency and conditional causal exchangeability identify $\psi={\mathbb{E}} m(X)$. The density $p$ is bounded above and below and the nuisance functions belong to the catalogue Hölder balls. For $s>0$ put $t(s)=\min\{s,1\}$. All constants below may depend on fixed $d,L,\epsilon$ and supplied smoothness, but not $n$.

Cinelli et al. [challenge] explicitly distinguish implementation from statistical optimality. Higher-order estimators already form an established literature [robins]; the first-order construction here does not displace their sharper rough-regime results. The following finite-sample derivations and operation counts are elementary restricted results, not claims of a new general method.

Histogram fitting with random counts

Split the first $3N$ observations into three independent groups of size $N=\lfloor n/3\rfloor$, ignoring at most two remaining observations. In the first group partition the cube into $J_e^d$ half-open equal cubes (with the final faces included) and estimate $e$ by cell means of $W$, using $1/2$ in empty cells and clipping to $[\epsilon,1-\epsilon]$. In the second group use $J_m^d$ cells and estimate $m$ by the cell mean of $Y$ among $W=1$, using $1/2$ when there is no treated observation. Write the fitted functions as $\widehat e,\widehat m$.
Lemma (Integrated histogram error): If $e\in\mathcal H_d^{s_e}(L)$ and $m\in\mathcal H_d^{s_m}(L)$, then, with $t_e=t(s_e),t_m=t(s_m)$, \[ {\mathbb{E}}\|\widehat e-e\|_{L^2(P_X)}^2 \leq C(J_e^{-2t_e}+J_e^d/N),\qquad {\mathbb{E}}\|\widehat m-m\|_{L^2(P_X)}^2 \leq C(J_m^{-2t_m}+J_m^d/N). \] No smoothness of $p$ is needed for these inequalities.
Proof: On a cell of side $J^{-1}$, either Hölder continuity ($s\leq1$) or the bound on first derivatives ($s>1$) bounds the oscillation of the relevant regression function by $C_{d,L}J^{-t(s)}$. Let the cell probability be $\pi$ and the probability of a qualifying observation be $\rho$. For fitting $e$, all observations qualify and $\rho=\pi$. For fitting $m$, qualification means $W=1$ and $\rho\geq\epsilon\pi$. The number $D$ of qualifying observations is binomial $(N,\rho)$. Conditional on $D=j>0$, their responses are independent draws from the conditional cell law; their mean has variance at most $1/(4j)$ and expectation equal to the corresponding cell mean. That mean differs from the regression function at every point of the cell by at most its oscillation, including the treatment-weighted cell mean used for $m$. For $0<\rho\leq1$, \[ {\mathbb{E}}\frac1{D+1}=\frac{1-(1-\rho)^{N+1}}{(N+1)\rho},\quad {\mathbb{E}}\frac{{\mathbf 1}\{D>0\}}D\leq\frac2{(N+1)\rho},\quad P(D=0)\leq\frac1{(N+1)\rho}. \] The first identity follows by integrating ${\mathbb{E}} z^D=(1-\rho+\rho z)^N$ over $z\in[0,1]$; the second uses $1/j\leq2/(j+1)$, and the third follows from the first because its integrand expectation includes $P(D=0)$. The last justification can equivalently be checked directly from $(N+1)\rho(1-\rho)^N\leq1$. Squared bias and variance thus give integrated cell risk bounded by $C\pi J^{-2t}+C\pi/[(N+1)\rho]$; the empty-cell squared error is at most one. Sum over cells, using $\pi/\rho\leq1/\epsilon$. Zero-mass cells contribute zero. Clipping $\widehat e$ weakly reduces its pointwise distance from the true $e$, completing the proof. ∎
Theorem (A finite-sample doubly robust bound): Use the third group to compute \[ \widehat\psi=\frac1N\sum_{i\in I_3} \left[\widehat m(X_i)+\frac{W_i}{\widehat e(X_i)} (Y_i-\widehat m(X_i))\right]. \] Then \[ {\mathbb{E}}(\widehat\psi-\psi)^2\leq \frac{(1+\epsilon^{-1})^2}{N} +\epsilon^{-2}{\mathbb{E}}\|\widehat e-e\|_2^2 {\mathbb{E}}\|\widehat m-m\|_2^2. \] Choose dyadic integers $J_j$ within a factor two of $N^{1/(2t_j+d)}$, for $j=e,m$. The resulting risk is at most \[ C\left[n^{-1}+ n^{-2t_e/(2t_e+d)-2t_m/(2t_m+d)}\right] \] for $n\geq6$, and the clipped estimator has no larger risk.
Proof: Condition on both fitting groups. The score has absolute value at most $1+\epsilon^{-1}$, so its sample-mean conditional variance is bounded by the first displayed term. Direct conditional expectation gives the exact bias \[ B=\int\frac{(\widehat e-e)(\widehat m-m)}{\widehat e}\,dP_X. \] Cauchy--Schwarz yields $B^2\leq\epsilon^{-2}\|\widehat e-e\|_2^2\|\widehat m-m\|_2^2$. The two fitted functions are independent because they use separate groups, so the expectation of this product factors. Conditional variance plus squared conditional bias proves the first inequality. The lemma and the chosen grid sizes give each factor order $N^{-2t_j/(2t_j+d)}$. Since $N$ is comparable to $n$, the second inequality follows. Projection onto $[0,1]$ reduces error because $\psi\in[0,1]$. ∎

What exponent is actually established

Theorem (Minimax risk and the infimum runtime exponent): Suppose \[ \frac{t_e}{2t_e+d}+\frac{t_m}{2t_m+d}\geq\frac12. \tag{*} \] For the full smoothness class in the statement, assume it contains the constant subfamily $p=1,e=1/2,m=\theta$, $\theta\in[1/4,3/4]$, with $Y\mid X,W$ Bernoulli $(\theta)$. Then $c/n\leq R_n^\star\leq C/n$. The estimator above takes at most $C_d n\log(n+1)$ arithmetic RAM operations, counting fitting, tuning, storage initialization and evaluation. If inspecting an input observation costs at least one operation, no algorithm with uniform expected runtime at most $n^a$, $a<1$, can attain a constant multiple of $R_n^\star$. For every $a>1$ such a constant-factor algorithm exists. Consequently the infimum possible exponent is $1$; this does not assert attainment under the exact bound $T\leq n$ at $a=1$.
Proof: Condition (*) makes the preceding upper bound order $n^{-1}$. For the statistical lower bound take $\theta_0=1/2$ and $\theta_1=1/2+1/(4\sqrt n)$ in the constant subfamily. Only the $Y$ coordinates depend on $\theta$. The single-observation chi-square divergence of $P_1$ from $P_0$ is $4(\theta_1-\theta_0)^2$; independence gives $1+\chi^2(P_1^n,P_0^n)=(1+1/(4n))^n\leq e^{1/4}$. Cauchy--Schwarz gives ${\operatorname{TV}}\leq\frac12\sqrt{\chi^2}<1/2$. Classify an estimate by the nearer target. A classification error entails squared loss at least one quarter of the squared target separation, and the sum of testing errors is at least $1-{\operatorname{TV}}$. Thus maximum risk is at least a positive constant times $n^{-1}$.

Locate each coordinate in a dyadic grid by binary comparisons in $O(\log J_j)$ steps. Each observation requires $O(d\log(n+1))$ comparisons and a constant number of additions per visited cell. Each array has $J_j^d\leq C_d N$ entries, so initialization and division of cell sums by counts take $O(n)$ operations for fixed $d$. For each fixed $t_j$, include $\lambda_j=2^{2t_j+d}$ among the fixed real RAM constants. Choose the largest integer $k$ with $\lambda_j^k\leq N$ by repeated multiplication, and put $J_j=2^k$. This chooses the required resolution in $O(\log n)$ operations without a unit-cost exponentiation assumption. These are arithmetic-model constants, with no claim about their finite-precision bit representation. Evaluation on the third group has the same lookup cost. This proves the total count.

For the runtime lower bound the standard RAM model must charge access to data: an operation can reveal only a bounded number of previously unread records. Absorb that bound into constants. Under the constant subfamily every new record is an independent Bernoulli outcome together with ancillary variables. Adaptively chosen unread indices retain that law. The chain rule for relative entropy of the algorithm's stopped transcript therefore gives \[ {\operatorname{KL}}(Q_0\|Q_1)\leq {\mathbb{E}}_0 N_{\rm read}\, \operatorname{kl}(\theta_0,\theta_1) \leq C n^a(\theta_1-\theta_0)^2. \] To justify stopping, first truncate after $M$ operations, use the finite chain rule, then let $M\to\infty$; finite expected runtime implies almost sure stopping, and the entropy bound passes to the transcript by the increasing information partitions (or by its conditional log-likelihood sum). Choose target separation $c_0 n^{-a/2}$ with a fixed sufficiently small $c_0>0$. It stays inside $[1/4,3/4]$ and makes entropy at most $1/8$. Binary entropy satisfies $\operatorname{kl}(u,v)\geq2(u-v)^2$: its second derivative in $u$ is $1/[u(1-u)]\geq4$ and its value and first derivative vanish at $u=v$. Grouping likelihoods on an event by Jensen therefore proves ${\operatorname{TV}}(Q_0,Q_1)\leq\sqrt{{\operatorname{KL}}(Q_0\|Q_1)/2}\leq1/4$. The same testing argument gives risk $\geq c n^{-a}$, which is not $O(n^{-1})$ when $a<1$. Finally $C_dn\log(n+1)\leq n^a$ for all sufficiently large $n$ if $a>1$. On the finitely many smaller sample sizes output a constant in one operation; increasing the risk comparison constant covers them, since the statistical lower bound is positive. This respects the exact runtime constraint for every $n\geq1$. ∎

Verification and first unresolved gap

The two nuisance fits were deliberately made independent; otherwise factoring their squared $L^2$ errors would be unjustified. The operation bound covers all observations actually used and all array entries. In the runtime lower bound data access, not a purported hardness conjecture, is the source of the obstruction. Condition (*) is a sufficient first-order region, not the sharp boundary of parametric minimax risk. The first unresolved gap is an explicit optimal runtime construction and matching lower bound in smoothness regimes requiring higher-order bias correction, including the adaptation penalty when indices are unknown. Arithmetic counts here imply no bit-complexity claim.

References

  • [challenge] C. Cinelli, A. Feller, G. Imbens, E. Kennedy, S. Magliacane and J. Zubizarreta (2025), Challenges in Statistics: A Dozen Challenges in Causality and Causal Inference, arXiv:2508.17099v1, Section 3.5.3, Implementation. https://arxiv.org/html/2508.17099v1#S3.SS5.SSS3.
  • [robins] J. Robins, L. Li, E. Tchetgen and A. van der Vaart (2008), Higher order influence functions and minimax estimation of nonlinear functionals, arXiv:0805.3040. https://arxiv.org/abs/0805.3040.

Completion Estimate

25%

Go to contribution →