{'fraction_non_alphanumeric': 0.03647771605271201, 'fraction_numerical': 0.012512809449328515, 'mean_word_length': 5.103807747174367, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Empirical Gateaux Derivatives for Causal Inference

We study a constructive algorithm that approximates Gateaux derivatives for statistical functionals by finite-differencing, with a focus on causal inference functionals. We consider the case where probability distributions are not known a priori but also need to be estimated from data. These estimated distributions lead to empirical Gateaux derivatives, and we study the relationships between empirical, numerical, and analytical Gateaux derivatives. Starting with a case study of estimating the mean potential outcome (hence average treatment effect), we instantiate the exact relationship between finite-differences and the analytical Gateaux derivative. We then derive requirements on the rates of numerical approximation in perturbation and smoothing that preserve the statistical benefits of one-step adjustments, such as rate-double-robustness. We then study more complicated functionals such as dynamic treatment regimes and the linear-programming formulation for policy optimization in infinite-horizon Markov decision processes. The newfound ability to approximate bias adjustments in the presence of arbitrary constraints illustrates the usefulness of constructive approaches for Gateaux derivatives. We also find that the statistical structure of the functional (rate-double robustness) can permit less conservative rates of finite-difference approximation. This property, however, can be specific to particular functionals, e.g. it occurs for the mean potential outcome (hence average treatment effect) but not the infinite-horizon MDP policy value.


Introduction
Inferential targets in causal machine learning often take the form of statistical functionals of the data distribution. Examples include average treatment effects or average policy values for infinite-horizon offline reinforcement learning. Estimation of these statistical functionals can be vulnerable to the first-stage bias introduced by the estimation of nuisance functions such as outcome regressions or transition probabilities. However, by leveraging the causal structure, it is possible to derive bias-adjusted estimators that weaken the need for accurate estimation of these nuisance functions. The celebrated doubly robust estimator is one such bias adjustment, and there are general frameworks for deriving such estimators, including semiparametric efficiency, double robustness, Neyman orthogonality, and "debiased/double" machine learning [Bickel et al., 1993, Chernozhukov et al., 2016, Hines et al., 2022, Kennedy, 2016, Newey, 1994a, Robins et al., 1994. These general frameworks can also be applied to more complicated functionals in longitudinal causal inference and offline reinforcement learning (also known as offline policy evaluation/learning) [Bibaut et al., 2019, Chakraborty and Moodie, Jiang and Li, 2016, Kallus and Uehara, 2020, Kallus and Zhou, 2020, Tang et al., 2019, Thomas and Brunskill, 2016, Van der Laan and Robins, 2003.
A drawback of these general frameworks is that significant analytic effort is often necessary to obtain concrete estimators that are appropriate for particular situations. Indeed, a practitioner may be interested in novel variants of an estimand, or be working within a constrained class of functionals that are appropriate for a particular class of problems. For example, in constrained Markov Decision Processes (MDPs) [Altman, 1999], or optimization-based estimators more broadly, a practitioner can easily avail themselves of custom convex-optimization solvers that are computationally efficient and easy to deploy to add additional constraints. Unfortunately, such choices may require case-specific re-analysis to establish the desired statistical properties, and indeed deriving the actual estimator that yields bias adjustment may require significant analytical effort. It is therefore important to develop a suite of constructive, numerical or algorithmic methods that yield desired forms of bias adjustment. Such a suite would be complementary to analytic derivations.
To this end, we develop off-the-shelf procedures for estimating the statistical functionals used for bias adjustment in causal inference. We do so via numerical approximation of the Gateaux derivatives that underlie the general analytic frameworks, building on prior work in this area. We focus on interventional effects. Newey [1994b] suggested numerical differentiation to estimate Gateaux derivatives of functionals unavailable in closed form but approximated by computational procedures, such as the solution to a system of differential equations; this approach is later used in Hausman and Newey [1995]. Recent pedagogical surveys of Hines et al. [2022], Kennedy [2022] also emphasize the usefulness of Gateaux derivatives for deriving influence functions, an approach which is directly applicable for discrete data and generalizes to continuous distributions via a smoothing argument as discussed in Carone et al. [2018], Ichimura and Newey [2015]. In this spirit, we build on the line of work of Carone et al. [2018], Frangakis et al. [2015], Ichimura and Newey [2015], Newey [1994b], which proposes a constructive procedure by numerical approximation of a Gateaux derivative (in the sense of Hampel [1974], Huber [2004]) by finite differences.
Numerical differentiation can be a useful paradigm because evaluating the one-step estimator reduces to multiple evaluations of the black-box statistical functional. In other areas such as optimization, numerical differentiation underpins common subroutines that enable analysts to use optimization algorithms without specialized training in a particular modeling language or paradigm. 1 Although more advanced paradigms can compute exact gradients, they can require specialized training to use or implement.
In addition, the numerical differentiation viewpoint connects the conceptual interpretation of influence functions as qualitative sensitivity analysis (i.e. sensitivity of a functional to perturbations in the distribution) to the use of influence functions in causal inference. Typical presentations of the latter framework are not concretely connected to the former notion, where influence functions are well-known in robust statistics. For example, this former notion was applied in the machine learning setting to probe machine learning models [Koh and Liang, 2017]. We provide a concrete connection from the viewpoint of exactly quantifying the finite difference approximation error, which can be interesting independent of computerized estimation for algorithmic, pedagogical, or conceptual reasons.
In this paper, we build on the previously mentioned line of work on the use of numerical derivatives to approximate influence functions, focusing on the implications for statistical estimation. In more detail, we provide exact characterizations of empirical Gateaux derivatives, in contrast to numerical derivatives computed based on oracle knowledge of probability distributions. This work is focused on exact characterizations that can inform statistical and computational trade-offs. Although we expect that some of these characterizations may be apparent to experts, we also hope that this level of concreteness can be helpful to non-experts. We use these characterizations to study statistical and computational trade-offs: how statistical properties of functionals in causal inference can weaken computational approximation requirements, as well as more complicated functionals. These trade-offs could eventually be practically relevant because the finerscaled discretization required for lower approximation error can incur fundamental issues with floating-point computation that prohibit further gains in accuracy [Fornberg, 1981, Sauer, 2017.
Our contributions are as follows. We begin with a case study of the average potential outcome. The average treatment effect is the difference of the average potential outcomes, and our analysis applies also to the average treatment effect. We discuss the treated mean for simplicity; by symmetry the same results hold for the control mean. In this setting, we characterize the exact relationship between finite-differences and the analytical Gateaux derivative. This characterization helps concretize the proposal of empirical Gateaux derivatives by allowing us to study the rates of numerical approximation that can preserve statistical performance. Finally, with these concretizations in hand, we turn to more complicated motivating examples including a finite-horizon dynamic treatment regime and infinite-horizon reinforcement learning. We illustrate with the case of policy learning in finite-state, finite-action infinite-horizon Markov decision processes how this constructive approach to Gateaux derivatives can directly extend also to the case of custom constraints.

Problem Setup: From Numerical to Empirical Gateaux Derivatives
We begin by introducing our first example, that of the mean potential outcome, and the canonical doublyrobust estimator. We then briefly introduce the key objects in the more general framework of orthogonalizing or debiasing statistical functionals by Gateaux derivative adjustments via the one-step estimator. We also describe the numerical approximation of these adjustments via perturbed black-box evaluations of the statistical functional and plug-in estimation. After introducing these frameworks, we clarify what we mean by empirical rather than numerical Gateaux derivatives and delineate our specific research questions.
We let O ∼ P denote a draw of observations following the distribution P belonging to a statistical model M. We will focus on the estimation of statistical functionals Ψ(P ), where Ψ : M → R q is pathwise differentiable. Throughout, we take as given the assumptions that grant causal identifiability (i.e., overlap/positivity, unconfoundedness/ignorability).
The mean potential outcome and augmented inverse propensity weighting. We very briefly overview the celebrated doubly robust estimate of the mean potential outcome. The average treatment effect is simply the difference of these doubly robust estimators for the treated and control means; for simplicity, we discuss the treated mean only.
Example 1 (Mean potential outcome). The observation O = (X, A, Y ) includes covariate X, treatment A ∈ {0, 1}, and outcome Y ∈ R. The statistical functional corresponding to the mean potential outcome is: Plug-in estimation of the conditional expectation is termed regression adjustment. The overall functional is also referred to as inverse propensity weighting, e(X) , where e(x) = P (A = 1 | X = x) is the propensity score. (We omit explicit discussion of assumptions such as consistency and ignorability, under which these identification results hold). Given two identifying functionals for the mean potential outcome, a natural question is: can these identification approaches be combined in some way to improve estimation? The doubly-robust estimator of Robins et al. [1994] achieves this as follows: Appealing properties such as the mixed-bias property and the rate double robustness properoty can be verified readily from its functional form. This provides a canonical example of one-step adjustment by influence functions that applies more broadly.
Numerical Gateaux derivatives. Although this definition of influence function via perturbation is not quite enough for semiparametric inference with continuous distributions, Carone et al. [2018] and Ichimura and Newey [2015] show how to obtain absolute continuity in semiparametric models by taking an additional limit as a smoothed perturbation converges to the point mass. That is, they replace the Dirac delta measure δ o with a smoothed distributionδ λ o (o), a function satisfying g(u)dδ λ o (u) = 1 and dominated by P . A common choice forδ λ o has the formδ λ where K λ is a bandwidth-normalized kernel function and notation is overloaded to also designate product kernels. We then define a perturbation in the direction of a generic observation o i as follows: Carone et al. [2018], Frangakis et al. [2015] discuss numerical approximation by finite differences: for a fixed , λ with λ, evaluate the numerical derivativeφ(o; P ) = −1 (Ψ(P ,λ ) − Ψ(P )). Further, Hines et al. [2022], Kennedy [2022] highlight influence function derivations based on the analytic Gateaux derivative.
Carone et al. [2018] establish the limit of P o i ,λ when and λ tend to zero is precisely the Gateaux derivative. and they provided generic approximation rates in ( , λ) for a multi-point finite difference scheme. These rates establish sufficient conditions for approximation error that ensure validity of the numerical approximation, but can be conservative. For a two-point approximation with a uniform kernel, the rate from Carone et al. [2018, Thm. 5] is λ −d + λ −2 . We will strengthen this result later in the paper.
Empirical Gateaux derivatives. In the context of statistical estimation, the underlying probability distributions also need to be estimated. We consider plug-in estimation of the statistical functional via density estimation of the collection of joint probability estimatesp(o) that compriseP . For example, a plug-in representation of the mean potential outcome is as follows: It is important to distinguish between the empirical Gateaux derivative computed by numerical approximation with estimated probability densities (e.g.,P (Y, A = 1, X)), the numerical Gateaux derivative obtained by finite differencing based on the true probability densities, and the analytic Gateaux derivative: empirical derivative at smoothed and estimated distributions, numerical derivative at smoothed and true distributions, analytical Gateaux derivative.
It is the first of these representations that is appropriate for serving as a statistical estimator of influence functions. In particular, the one-step estimator is justified by the following asymptotically linear representation [Pfanzagl, 1990]: See Appendix A for further discussion of this expansion. Given these definitions, we now turn to the major mathematical and statistical questions that we address in this paper: Q1: How exactly does the empirical Gateaux derivative approximate the analytical derivative?
Q2: What are the required rates of numerical approximation in perturbation and smoothing based on ( , λ) that preserve the beneficial statistical properties of the constructed estimator?
In Section 4, we study these questions in the simple case of our running example: estimation of the mean potential outcome. In that setting we can provide a fine-grain comparison of the finite-difference estimator and the classical doubly robust estimator. In Section 5, we study more complicated functionals to highlight more general applicability.
Example 2 (Example 1, continued). The influence function [Hahn, 1998, Robins et al., 1994 for a mean potential outcome is: Hence the canonical doubly-robust estimator (augmented inverse propensity weighting, or AIPW) can be viewed as a one-step adjustment.
We conclude our introductory remarks by noting some limitations of our study. We omit discussion of semiparametric efficiency since we focus on computing influence functions without the consideration of model restrictions and tangent spaces. That is, we focus on Gateaux derivatives under a nonparametric model. Note that Carone et al. [2018] suggest that one can handle the semiparametric case by evaluating a Gateaux derivative and projecting onto the model space directly; such an investigation is out of scope of the current paper.
Throughout, we assume the functional is pathwise differentiable. This assumption does warrant caution, because this condition can pose a fundamental barrier to efforts for "automation"; hence empirical Gateaux derivatives are intended to augment but not supplant or replace analytical expertise. However, under the (stronger) assumption of Fréchet rather than Gateaux differentiability, previous work gives sufficient conditions that imply the sufficiency of our assumption. Ichimura and Newey [2015, sec. 4] provides a generic sufficient condition for asymptotic linearity assuming additionally that the Fréchet derivative is Lipschitz, while Kandasamy et al. [2015] show the remainder term is second order; both under convergence rate conditions on the input densities. See Appendix D.2 for more discussion, including a class of problems where general guarantees are possible, i.e., strongly convex stochastic programs.

Related Work
In this section, we discuss the most closely related directions but do not attempt to provide a survey or overview of influence functions or semiparametric statistics. See Bickel et al. [1993], Fisher and Kennedy [2021], Kennedy [2022], Tsiatis [2006] or the appendix for a broader overview and further details. In Appendix D.1, for additional context, we also overview work on numerical derivatives in optimization and statistical machine learning more broadly.
Numerical approximation of Gateaux derivatives. Recent pedagogical reviews [Fisher and Kennedy, 2021, Hines et al., 2022, Kennedy, 2022 emphasize the Gateaux derivative and its use as a heuristic tool for deriving influence functions. While our point of departure is a finite-difference approach to the former, we also consider more complicated optimization functionals, with a focus on how dual variables characterize the Gateaux derivatives therein. Carone et al. [2018] derive high-level sufficient conditions for approximation error via general results for finite-difference approximations. Relative to those general sufficient conditions, we focus on specific improvements in approximation error for various functionals. For the specific case of the mean potential outcome/average treatment effect, Frangakis et al. [2015] exploit additional structural knowledge. In particular, they interchange integration and differentiation as a result of an analytical insight and obtain functional evaluations in terms of the conditional expectation. For more complex functionals, making these determinations is similar to analytically deriving the influence function; i.e., it is less "automatic." In general, it may be difficult to deduce the form of the nuisance functions without using knowledge of the estimand or working with specific graphical or algorithmic structure. The approach in Appendix B of Bravo et al. [2020] considers the specific case of debiasing moment conditions with respect to first stages but also implicitly leverages this change in order of integration.
Other work on automating causal inference. Recent work in causal inference aims at "automatic" or "computerized" semiparametric inference via a variety of methods [Bhattacharya et al., 2020, Chernozhukov et al., 2020, Farrell et al., 2021, Jung et al., 2021b, Kandasamy et al., 2015. To summarize briefly the difference with our approach, our approach focuses on the use of numerical derivatives, and requires only a zeroth-order computation oracle of perturbed functional evaluations, without additional nuisance estimation. This approach can be particularly helpful for incorporating debiased causal inference in the design and analysis of algorithms. On the other hand, it may not lead to improvements in estimation performance relative to other approaches that modify the mode of estimation or estimate additional nuisance functions.
In more detail, Kandasamy et al. [2015] outline a general approach for Gateaux derivative estimation based on leave-one-out techniques. Chernozhukov et al. [2021] build on the variational characterizations of Ichimura and Newey [2022] via Riesz representation for generalized regression residuals. Chernozhukov et al. [2020,2022] and other works further develop "automatic" estimation based on this approach. Our references throughout to Ichimura and Newey [2015] are primarily for the smoothing arguments required in a semiparametric setting. Farrell et al.
[2021] develop a framework for modeling heterogeneity with deep networks, similar to a functional coefficient model that admits an influence function computable by automatic differentiation (or numerical differentiation of Hessians/Jacobians for a given parametrization). Bhattacharya et al. [2020] focus on graphical algorithms and characterizations of the efficient influence function for identified latent variable models, for ATE-like functionals, given a causal graph. Alaa and Van Der Schaar [2019] develop approaches for validating causal machine learning methods with influence functions.
Influnce functions or Gateaux derivatives in other settings. More broadly in machine learning, influence functions have been used for qualitative sensitivity analysis or bias correction [Giordano et al., 2019a,b, Koh and Liang, 2017, Wilson et al., 2020; further connections to numerical methods could be interesting. However, the use of influence functions in the causal setting is different because one computes influence functions of the statistical functional with respect to perturbations in the distribution, rather than perturbing the moment function that identifies the estimator. Recently, in data-driven optimization, finite-difference approaches for Gateaux derivatives have been used to debias in-sample bias of the stochastic optimization functional in general; see, e.g., Gupta et al. [2021]. Guo et al. [2022] build on the finite-difference method of Ito et al. [2018] for debiasing optimization with causal outcome coefficients but without characterization of the adjustment.

Empirical Gateaux Derivative of the Mean Potential Outcome
To illustrate the numerical approach to Gateaux derivatives, we first work through the running example of the mean potential outcome and indicate precisely how the numerical differentiation incurs error relative to the analytical derivative, which in this case is the doubly robust or augmented IPW (AIPW) estimator.

Relating the empirical and analytical Gateaux derivatives
Setup. We target the estimation of a mean potential outcome as in Example 1. For the smoothed perturbation, we specialize to kernel functions K λ (u) from kernel density estimation satisfying ∞ −∞ K(u)du = 1 (normalization to a probability density) and K(−u) = K(u), for all u (symmetry). 3 The mean potential outcome is identified as follows, We obtain Ψ(P ) by plugging in thep probability density estimates. Note this is distinct from the plug-in estimator in empirical process theory, which plugs in the empirical CDF; we plug in statistical estimates. See Algorithm 1 for a summary of the procedure.
Comparison of empirical Gateaux derivative and augmented inverse propensity weighting. We have argued that, in the context of statistical estimation, approximating Gateaux derivatives requires plug-in estimation of the statistical functional with estimated probabilities. Towards resolving Q1, we observe that the perturbed probability densities induce the functional form of the nuisance functions in the analytical derivative. This is in contrast to the analytical derivative, which describes the estimand in terms of, e.g., conditional expectations, and allows the analyst to choose the functional form for estimating the nuisance function.
We will denote theP -induced conditional expectation as For example, the density-induced conditional expectation of kernel density estimates is exactly Nadaraya-Watson regression. 4 Given this, we study empirical derivatives with kernel density estimation. We study KDEs because they provide a classic textbook example of density estimates that induce conditional expectation estimators; our analysis is modular in that other density estimation approaches can be used and their rates of convergence would simply be substituted in Lemma 1. In Appendix D we discuss other choices of density estimates: any "modern" generative modeling approach/density estimation approach can be used for plug-in evaluation of the statistical functional. Of course, all such nonparametric estimation approaches suffer from the curse of dimensionality in general so additional structure is required in general to satisfy the product-rate conditions. Our characterization will show exactly how plug-in evaluation with finite differences is nearly equivalent to evaluating the doubly-robust estimator with a Nadaraya-Watson regressor as the conditional outcome regression for E[Y | A = 1, X] and a kernel density estimate of the propensity score. Indeed, it is equivalent up to an additive bias that arises due to the smoothed perturbation distribution, and what we will term a smoothed nuisance evaluation. We next discuss what this smoothed nuisance function is, and how we study its approximation error relative to estimation error of the induced nuisance.
Integrating the smoothed perturbation induces smoothed nuisance functions. Integrating with respect to the smoothed perturbation distribution can be interpreted as a smoothed evaluation of the resulting nuisance functions. We define theδ λ x 0 -smoothed conditional expectation which smooths the evaluation point, i.e., smooths evaluation around x 0 rather than precisely at x 0 .
A similar definition appears in the literature on nonregular inference or the localization of a global functional [Newey, 1994b]. For example, van der Laan et al. [2018] proposes smoothing a nonregular target estimand, such as a dose-response curve or density evaluated at a point x 0 , by conducting a smoothed evaluation locally around x 0 . In that work, a function evaluated at x 0 , g(x 0 ), is approximated with a kernel smooth over x in a neighborhood of [2021a] use a similar smoothing/localization device but for a different problem in general. The following 4 For example, Algorithm 1 Empirical Gateaux derivatives 1: Inputs:P probability density estimates, data {O i } 1:n , perturbation and smoothing ( , λ) parameters, functional Ψ 2: for i = 1, . . . n do 3: ) 5: end for 6: Output: one-step adjusted estimate from empirical Gateaux derivatives, Ψ(P ) + 1 n iφ (o i ) lemma, similar to van der Laan et al. [2018, Lemma 25.1], summarizes convergence rate implications of these smoothed nuisances; these are direct consequences of typical analyses of kernel density estimation.
Assume thatδ λ is a product kernel in a Hölder class of degree β, and |u| β |δ λ u | du < ∞ and u sδλ u du = 0 for s ≤ β. Assume Y is bounded. Then the perturbed nuisances satisfy the following rates: The perturbed nuisances are asymptotically consistent, so long asp induces asymptotically consistent nuisances and we fix as a sequence vanishing in n at a rate appropriate to counteract the growth in λ due to "roughness"; i.e., the variability ofδ λ . The difference between the function evaluation and its smoothed nuisance arises from the bias analysis of kernel density estimators and the dimension-dependence of ensuring absolute continuity, but it is a non-stochastic, entirely deterministic argument. While the analysis improves mildly in the dimension dependence on λ compared to the general analysis of Carone et al. [2018, Thm. 5], this still highlights the generally unfavorable dependence of the numerical approach on the dimension due to smoothing.
Our first proposition establishes the exact approximation error induced by finite differencing.
Proposition 1 (Gateaux derivative of probability density representation). Consider perturbations in the direction of o i = (x i , a i , y i ). Letp(x),p(A = 1, x),p(y, A = 1, x) denote kernel density estimates. Let Ψ(P ) be as in Equation (2): We compare the expansion of Proposition 1 to the canonical AIPW estimator in Example 2. We specialize to the case of a uniform kernel forδ λ y i (y) so that yδ λ y i (y)dy = y i . We next observe that the convolution with respect to the smoothing kernel introduces additive bias relative to evaluating terms in the expression at the observation o i .

Corollary 1. When the perturbation observation is the observation datapoint
Lemma 1 allows us to further simplify and conclude that the decomposition of Proposition 1 is close to the canonical AIPW estimator up to the O(λ β ) bias induced by smoothed evaluation.

Implications for estimation
Next, we use our characterization in Proposition 1, Lemma 1, and Corollary 1 to study the statistical properties of our computational/numerical approximation. Corollary 1 arises from deterministic equivalences and allows us to study the relationship to the typical "oracle" AIPW estimator; albeit with induced -perturbed nuisances.
An appealing property of bias-adjusted treatment effect estimation is rate double robustness: because we incur the product of convergence rates of the nuisances, we may enjoy parametric n − 1 2 rate convergence of the target functional while nuisances converge at slower rates, for example at n − 1 4 ; see Chernozhukov et al. [2018], Rotnitzky et al. [2021]. We use the rate double robustness property of the target and Lemma 1 to infer the required rate conditions on ( , λ) that retain the beneficial statistical properties of the canonical AIPW estimator. The next result combines our previous characterizations with the standard analysis of AIPW.
Assumption 1 (Regularity conditions). Assume the following regularity conditions hold: (iii) The estimates p(y, a, x), p(a, x), p(x) belong to a Donsker function class.
Theorem 1 (Rate double robustness). Consider the one step-estimator with the empirical Gateaux derivative. Under Assumption 1, when λ −d/2 = o(n − max(rµ,re) ), and λ β = o(n − 1 2 ): We now interpret the qualitative estimation implications of the analysis. Theorem 1 states the rate conditions required of numerical approximation in , λ to preserve the statistical property of rate double robustness. Assuming that the nuisances induced by the probability density estimates satisfy a product-rate condition that the product of their RMSEs is faster than the parametric O(n − 1 2 ) rate, the form of Proposition 1 suggests that the rate of RMSE convergence of perturbed nuisances is required to be faster than the rates of the unperturbed nuisance functions, max(r µ , r e ), which can be a slower rate than the = o(n − 1 2 ) implied by the generic analysis of finite differences. For a direct comparison, the rate of Carone et al. [2018, Thm. 5] suggests a rate requirement of λ −d 1 s + λ −2 = o(n − 1 2 ). In general, the dimension dependence on λ in Lemma 1 will prevent setting a slower rate for for nontrivial dimensionality. Note that without the smoothing requirement, the potential improvement in the rate of could be on the order of generic rate improvements implied by a central difference scheme. That is, we would have attained the approximation benefits of generically approximation-improved finite-difference schemes via the specialized statistical structure of the target adjustments. Our analysis does make use of our characterization in Corollary 1; i.e., a deterministic approximation error from AIPW with induced nuisance functions.

Applications: Mean Potential Outcome and Dynamic Settings
As empirical illustrations, we first show a plot of numerical approximation error ranging over ( , λ) in the simple case of the mean potential outcome. We then apply the tools that we derived in the previous section to analyze more complicated settings. We study longitudinal data (dynamic treatment regime) and the policy linear program for an infinite-horizon Markov decision with linear constraints, with attention to how different functionals have different interactions between approximation error and statistical properties.

Mean potential outcome
We first conduct a validation of our decomposition and theoretical characterization. Our data-generating process includes a piecewise linear outcome model (such that kernel regression is misspecified). We conduct plug-in evaluation of the mean potential outcome with kernel density estimates. In Figure 1a we include an ( , λ)-plot for the simple case of AIPW (see Carone et al. [2018] for more discussion and examples). We consider a one-dimensional case with uniformly distributed X, piecewise-linear Y , and smooth propensity scores that are logistic in sin(X). We use n = 500 and fix the bandwidth h = 0.05. We use colors to denote magnitude of the mean absolute error, included in text on the heatmap. Without loss of generality, we study the estimation of a mean under missingness, E[Y (1)]. Figure 1b illustrates the estimation error of various strategies with the comparable kernel-based estimates (DM is regression adjustment). Note that even in this simple setting AIPW offers benefits and sample efficiency relative to DM. We include further numerical experiments in the Appendix.

Dynamic treatment regimes
We study the estimation properties of the empirical Gateaux derivative for multi-stage DTR from the probability density representation. In T -stage dynamic treatment regimes, the causal quantity of interest is the mean potential outcome E [Y (ā)], whereā = (a 0 , . . . , a T ) is the (deterministic) treatment strategy. Assuming Y (ā) is sequentially ignorable given the treatment and covariate history, (Ā t ,X t ), at each time t, this causal quantity can be identified by the g-formula, To derive its influence function, we take the empirical Gateaux derivative by considering the perturbation in the direction of Below we characterize how this empirical Gateaux derivative at the smoothed distribution differs from the one at the (unsmoothed) estimated distribution. The following result is analogous to Proposition 1 and Corollary 1 but is extended to the dynamic treatment regime.
The expansion verifies that the requirements for numerical approximation are similar in as in the case of a single-timestep mean potential outcome, despite the additional complexity of the dynamic treatment regime function as a nested expectation. Again, smoothing incurs overall unfavorable dependence in the dimension, as the extension of Lemma 1 to EP [Y |Ā,X] also incurs unfavorable dependence on the dimension. The history-dependent nuisances imply that the dimension additionally grows in the time horizon.

Perturbations for (constrained) infinite-horizon off-policy evaluation
We show how this framework can be applied to derive generic bias corrections for infinite-horizon offpolicy optimization in offline reinforcement learning, via the canonical linear programming characterization [De Farias andVan Roy, 2003, Puterman, 2014]. Our application is relevant to a recent line of work on offpolicy evaluation and learning (OPE/L) [Bibaut et al., 2019, Jiang and Li, 2016, Kallus and Uehara, 2022, Kallus and Zhou, 2020, Tang et al., 2019, Thomas and Brunskill, 2016.
Setup. We recall the classical linear programming formulation of finding an optimal policy in a tabular (finite-state, finite-action) Markov decision process, due to De Farias and Van Roy [2003], Puterman [2014]. -Yadkori et al., 2019, Lee and He, 2018, Serrano and Neu, 2020, Wang, 2020. The dynamics follow an infinite-horizon Markov decision process with discount factor γ. The offline dataset is comprised of trajectories of state, action, next state observations: {(s 0 i , a 0 i , . . . , s t i , . . . , s T i )} 1:N . We derive our bias adjustment under a statistical model where the stationary distribution factorizes as p(s, a, s ) = d(s, a)P (s | s, a). From the joint state-action-state occupancy distribution p(s, a, s ), we estimate the marginalization d(s, a), the stationary state-action occupancy distribution, which we use to estimate the transition probability P (s | a, s) = p(s,a,s ) d(s,a) . Let µ 0 (s) denote the initial state distribution, estimated from offline data and assumed to be equivalent in the original confounded and evaluation setting. Because we focus on the tabular case with discrete distributions, we do not require smoothing.

Recent work has revisited this formulation with interest in developing primal-dual algorithms [Abbasi
The optimal policy linear program. The optimal value function solves the following linear program in a finite-state and finite-action setting; the objective value at optimality is the policy value. Here, P a ∈ R |S|×|S| is transition matrix at action a, with P a (s, s ) = P (s | s, a), µ 0 (s) is the marginal occupancy distribution of the initial state. We assume, as is fairly common in theoretical analyses of reinforcement learning, that the reward function is a known function of state and action, r(s, a). Let r a = {r(s, a)} s∈S be the vector of reward values per state, for fixed action a. We have: The dual formulation is well known to parametrize the stationary occupancy probabilities of the optimal policy (at optimality): We verify the bias correction from the Gateaux derivative d d Ψ D (P )

=0
. Note that evaluating this functional at the perturbed distribution corresponds to evaluating the linear program with perturbed left-hand side coefficients; i.e., perturbing the constraints of the coefficient matrix as well as the objective (marginalizing over a perturbed initial state distribution. We evaluate the Gateaux derivative in the direction of (s,ã,s ). We assume nondegeneracy (unique optimal solution); this can be relaxed by using a more sophisticated perturbation analysis.
Assumption 2. The optimal solution is unique; there is no degeneracy.
Proposition 3. Perturb in the direction of a generic observation (s,ã,s ) : Sketch of argument. The statement of Proposition 3 is not surprising per se, since we can immediately verify the mean-zero sample version of the Bellman residual, and so one might plausibly derive this form from double robustness. Rather, we provide an alternative argument that applies perturbation analysis of optimization programs with the specific perturbations from the finite difference. We build on the sensitivity analysis characterization of Freund [1985], which studies perturbations of linear programs with respect to constraint matrix perturbations. For small enough, the active basis remains the same. When the perturbation matrix is P + G for some matrix G, Freund [1985] notes that the derivative of the optimal value is d ( ) = −πGx (where (π, x) are the dual and primal optimal solutions, respectively) can be obtained from the derivative of a matrix inverse, which yields higher-order expansions. We generalize this argument for the specific form of the perturbed transition probabilities. We next compute the approximation error that arises from evaluating finite differences, so long as is small enough to maintain the same active basis, which is empirically verifiable.
Proposition 4 (Error analysis from finite differencing). Perturb in the direction of a generic observation (s,ã,s ). Then: In this setting, unlike the earlier case of the treatment mean in Corollary 1 and Theorem 1, even in a favorable discrete-state case without smoothing requirements, it does not appear that rate double robustness would admit weaker numerical requirements on . Ultimately, this arises because finite-difference approximation in this case involves perturbed matrix inversion, rather than composition of sum and product rules of differentiation as in the case of the previous functional. In the analysis, Equation (45) shows that V (s) − V (s) = O( ), so that Proposition 4 holds with the unperturbed nuisances and the same order of approximation error.
Interpreting the finite-difference approach. We take a moment to translate the benefits of a finitedifference approach to estimating a one-step adjustment in the context of reinforcement learning. In reinforcement learning, a common paradigm is that of "model-based" reinforcement learning; i.e., learning a dynamics model (estimating transition distribution) and plugging it into a planning problem or the linear programming formulation. A benefit of finite differences is that one can approximate doubly-robust adjustments by additional black-box function evaluations that they would already be computing, without estimating new nuisance functions.
Generalization to constraints. We argue the main benefit of this approach is that by computing these finite differences, we can obtain generic bias corrections via a more broadly applicable argument that does not require interpretation of the closed-form solution to an optimization problem. For example, consider a relevant subclass of constrained Markov decision processes comprised of additional constraints on the policy variables: the linear programming formulation is particularly appealing because linear constraints can be directly added [Altman, 1999]. Another benefit of the finite-differences approach is that one-step estimation can be obtained by evaluating the functional at feasible points. In our case, this corresponds to optimizing with respect to feasible probability distributions.

Example 3 (Constrained policy optimization).
Ψ P (P ) = max a∈A µ a r a : The additional constraint µ ∈ P reflects linearly representable constraints on state-occupancy under the optimal policy. For example, relevant constraints could require the occupancy in certain risky states to be bounded. Paternain et al. [2019] develops a framework with a "caution function," a convex function on the state-action occupancy measures; this analysis applies to linear caution functions.
To compare to related work, we note that double robustness has been studied (for policy evaluation) in Tang et al. [2019]. Our derivation via perturbation analysis is similar to Tang et al. [2019, Thm. 4.1]'s illustration via the Lagrangian. This is unsurprising since the Lagrangian provides the dual characterization that also underlies perturbation analysis. Relative to the doubly robust estimators derived elsewhere, this influence function may be higher variance.

Conclusions
We have presented a constructive algorithm that approximates Gateaux derivatives finite differencing, with a focus on the statistical functionals used in causal inference. There are several directions which are natural follow-ups to this work. First, the smoothing perturbation, although it restores absolute continuity, introduces technical barriers in the implied approximation rates with regards to dimension dependence. Other arguments that reduce the dependence on smoothing would be of interest. Regarding Q2 (the preservation of the beneficial statistical properties of estimators under numerical approximation), our analysis of different functionals surfaced fine-grained differences in opportunities for weakening requirements for computational approximation: Theorem 1 and Proposition 4 have different, estimand-dependent implications.
It would also be interesting to generalize our comparisons of finite versus analytical differences to functionals such as Z-and M-estimators. Another useful direction would be to explore methods that are fully adaptive, not requiring any prior knowledge of beneficial statistical properties. Finally, our work has been mainly theoretical, and it would be useful to conduct further empirical comparisons of our approach based on numerical differentiation to various alternatives. Appendix: Empirical Gateaux Derivatives for Causal Inference

A Overview and background for influence functions
For a broad overview of influence functions and their role in causal inference, see Hines et al. [2022], Kennedy [2022], Tsiatis [2006]. We focus on functionals that admit an asymptotically linear representation: where the influence function, φ(u;P ), captures the first-order bias in the functional. That is, we focus on functionals that are smooth enough to admit a von Mises expansion with a mean-zero, finite-variance influence function and a second-order remainder term R 2 (P 1 , P ).

B.0.2 Proofs
Proof of Proposition 1. We have: The first term can be expanded using the decomposition of Lemma 2. We recognize the second term as ẼP [Y | A = 1, X = x i ] − Ψ(P ) by definition of the perturbed nuisances. We obtain: Proof of Lemma 2. Verify by direct computation. For simplicity we consider a single summand in the sum (the claim follows by applying the same argument to every summand). Recall thatP (A = 1, X j ) = (1 − )p(A = 1, X j ) + Ĩ (o \y )ĩ and expanding the latter yields: Proof of Lemma 1. To bound E[(ẽ (X) − e(X)) 2 ], argue similarly and observe that sinceẽ =p (A = 1 | x) =p (A = 1, x)/p(x), we have: In the second-to-last inequality we lower boundp(x) ≥ ν and upper boundp(A = 1, x) ≤ 1. In the last inequality we apply a change of variables to identify the term that arises in analysis of variance of kernel density estimates. We incur the λ −d dependence since we assumed that multivariate kernels are product kernels. The roughness term R(δ λ ) = δ 2 (u)du depends on the particularδ kernel function and is a kernel-dependent constant. Hence the term scales overall as λ −d as in standard analyses of variance of kernel density estimates; see e.g. Ullah and Pagan [1999, pg. 23-26].
Next observe that: where the first inequality comes from the triangle inequality, and the second equality from arguments similar to those in the previous bound.
Proof of Theorem 1. Letẽ denote the propensity and conditional outcome nuisances arising from plug-in estimation ofP . Letẽ denote the nuisance arising from plug-in estimation at the perturbed distribution. First observe that by Corollary 1: This decomposition characterizes convergence of an AIPW estimator with possibly biased but asymptotically consistent nuisancesμ,ẽ in terms of the AIPW estimator with oracle nuisances-a standard analysis-and an additional term that vanishes on the order of (assuming strong overlap and boundedness). Define Γ * (O) as the doubly-robust score for AIPW with oracle nuisances e * , µ * , and let Γ * j denote its empirical counterpart at O j = (X k , A k , Y k ): We will study the decomposition relative to the AIPW estimand with oracle nuisance estimates, with Γ(O; e, µ) denoting the above score but with nuisance functions e, µ: where the second equality follows by Proposition 1, and the third by adding/subtracting E[Γ(O;ẽ ,μ)].
The first term, (E n − E) [Γ(O;ẽ ,μ)] is the standard empirical process term of the decomposition, obtained by the Donsker assumption and classical stability results for Donsker classes [Van der Vaart, 2000]: Hence, we have We obtain the product error conditions by standard decomposition of the remaining first term: The second equality follows by adding/subtracting I [A = 1] /ẽ , I [A = 1]μ/e * , I [A = 1] µ * /e * ). The third equality follows from iterated expectations. The last inequality follows from the Cauchy-Schwarz inequality and overlap. Therefore, must satisfy a rate so thatẽ converges at typical double robustness fast rates, such as n −1/4 , under a suitable sample splitting scheme. For the last term of Equation (10), note that by Lemma 2, the numerical error is multiplied by error rate of the outcome nuisance: by application of Cauchy-Schwarz inequality and strong overlap. The claim follows by application of Lemma 1 to deduce implications on from standard product-error rate requirements onμ,ẽ.

C Proofs for Section 5 C.1 Empirical Gateaux derivative for multi-stage dynamic treatment regimes
We study the estimation properties of the empirical Gateaux derivative for multi-stage DTR from the probability density representation. In T -stage dynamic treatment regimes, the causal quantity of interest is the mean potential outcome E [Yā], whereā = (a 0 , . . . , a T ) is the (deterministic) treatment strategy. Assuming Yā is sequentially ignorable given treatment and covariate history (Ā t ,X t ) at each time t, this causal quantity can be identified by the g-formula To derive its influence function, we take the empirical Gateaux derivative by considering the perturbation in the direction of o Below we characterize how this empirical Gateaux derivative at the smoothed distribution differs from the one at the (unsmoothed) estimated distribution. The following result is analogous to Proposition 1 but is extended to the dynamic treatment regimes.
For comparison, the exact efficient influence function under the non-Markovian dynamic treatment regime is This proposition performs a decomposition of the empirical derivative. It quantifies the additive bias due to the smoothed perturbation compared to the (raw) perturbation at the observation o i .
Proof. We first decompose the empirical Gateaux derivative.
Equation (13) is due to the definition of the smoothed perturbation. Equation (14) separates the effect of smoothed perturbation into two pieces: the evaluation of y and the evaluation of the density ratios. Equation (15)  Calculating the second term of Equation (15). We first decompose this term into two difference terms: The second difference term can be decomposed in a similar way as in Lemma 2.
y ·p We next analyze the first difference term by induction: Next we evaluate the key term in this inductionp (xs,ā s−1 ,x s−1 ) p (ā s−1 ,x s−1 )) −p (xs,ā s−1 ,x s−1 ) p(ā s−1 ,x s−1 )) : Proof of Proposition 4. We analyze the error introduced by finite differencing. Consider the dual formulation of Equation (6), and the perturbed formulation of Equation (40). Fix small enough such that the active basis remains the same under P and P ; let P π denote the submatrix of the transition matrix P (stacked for the optimization problem) corresponding to the active basis (equivalently, optimal policy). Define the entries of G such that P = P + G ,(analogous to Lemma 2): Let V π (s), V π (s) be the active subvector of the optimal solution achieving the optimal values of Ψ(P ), Ψ(P ), respectively. Then . (44) The optimal solutions V π , V π differ in the perturbation matrix; inverting the basic feasible solution, where invertibility is assured by nondegeneracy, gives that: We analyze V π − V π = ((I − γP π ) −1 − (I − γP π ) −1 )r π via a general result for perturbed matrices. Note that for two matrices A, B, we have (A + B) −1 = A −1 − A −1 B(A + B) −1 . Applying this recursively we obtain the power series expansion Applying this result to our setting, with A = (I − γP π ) and B = −γ G , we obtain V π − V π = γ(I − γP π ) −1 G π (I − γP π ) −1 r π + O( 2 ).
The statement follows by substituting in the above expression, evaluating based on definition of the entries of G π , and substituting the above and Equation (44) into Equation (42) of the proof of Proposition 3 and following that proof onwards.
D Additional discussion D.1 Additional related work on numerical differentiation Numerical derivatives, statistical machine learning, and optimization. Independently of causal inference, other works study the use of numerical gradients in the context of machine learning and stochastic gradients [Oktay et al., 2020], and statistical estimation and inference Li, 2018, Hong et al., 2015]. Analysis of finite-difference estimators has been of interest in (stochastic) zeroth-order optimization and derivative-free optimization [Duchi et al., 2015, Flaxman et al., 2005.

D.2 Discussion on pathwise differentiability
Assuming pathwise differentiability is a fundamental limitation of attempts to "automate" semiparametrics because it is a technical condition that requires a similar amount of specialized knowledge to understand and check as to derive an influence function itself. We can also provide an example of a broad class of problems where this assumption is satisfied in general: strongly convex stochastic programs. We find it useful to explicitly contrast these two different types of structural assumptions, because a reasonably informed analyst could check the condition of strong convexity, while checking Fréchet vs. Gateaux differentiability can be more difficult.
Regarding strongly convex stochastic programs, these provide an example where asymptotic linearity is established for the class as a whole but analysts can specify arbitrary functionals within this class. Duchi and Ruan [2021, Prop. 1] establish asymptotic linearity under primitive conditions of strong convexity (as they mention, these results are known in the literature on perturbation analysis of stochastic programs, Bonnans and Shapiro [1998]).
In the main text, we referenced sufficient conditions that Fréchet differentiability of the functional can imply asymptotic linearity and that the remainder is second-order. These are helpful, in contrast to current typical approaches that exhibit pathwise (non-)differentiability by checking the (non-)boundedness of the influence function. Kallus et al.
[2019] study a different problem, also make a stronger assumption of Fréchet differentiability, and hence provide some concrete examples (such as estimating equations with incomplete data) where this is satisfied. Therefore, we think that previously articulated sufficient conditions of Fréchet differentiability can be made even more useful in the context of additional results that establish Fréchet differentiability of functionals, based on the composition of simpler functionals, akin to identification of convex structure via composition rules and structure of known atoms Grant and Boyd [2008]. We think these are promising directions for future study, though beyond the scope of this paper.

D.3 Kernel density estimation practical considerations
Monte Carlo sampling for uniform kernel for smoothed perturbation. For computational purposes, it is helpful to evaluate the final integral over a uniform kernel in x by Monte Carlo sampling. 5 In particular, for perturbations in the direction of (x i , a i , y i ), we evaluate the integral using N M C Monte Carlo samples: (1− )p(A=1,x k )+ I[a i =1] Bandwidth for kernel density estimation of functionals. Bickel and Ritov [2003] discuss optimal bandwidth selection for estimation of functionals via plug-in of nonparametric density estimates. They show that asymptotic linearity/root-n consistency ofβ is possible under the optimal bandwidth for density estimation when the number of derivatives of the influence function is more than half the dimension of observation o; requiring a higher order kernel. (See Ichimura and Newey [2015] for additional discussion).
Other density estimation approaches. Although we have discussed estimation of the components ofP by nonparametric Nadaraya-Watson/kernel-density estimation for simplicity, other estimates of probability densities or generative models can be similarly used; we describe a few below. These approaches may work well empirically; we leave a thorough evaluation for future work.
Example 4 (Bayesian Dirichlet process mixture models). Hannah et al. [2011] develop Dirichlet process mixtures of generalized linear models with covariate density given by a mixture of exponential-family distributions, and conditional response density modeled via a generalized linear model (with a set of parameters for each mixture component). This specification of the probability densities implies a regression estimate.
Example 5 (Generative modeling with Gaussian mixture models). Jesson et al. [2021] use generative modeling with Gaussian mixture networks, and sampling from the generative model, to approximate the integral in order to scale-up Nadaraya-Watson regression to high-dimensional settings in sensitivity analysis. Similar approaches may be used to compute the density estimates in our setting.

E Details on simulation
We consider a piecewise linear outcome mean specification which is hard for nonparametric estimation. The propensity score is drawn as A ∼ Bern(sin(20X) + 0.5).
Details on computation Experiments were conducted on a MacBook Pro with 16g RAM.